HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 47(1) pp. 71–77 (2019) hjic.mk.uni-pannon.hu DOI: 10.33927/hjic-2019-11 COMPARISON OF THE APPROXIMATION METHODS FOR TIME-DELAY SYSTEMS: APPLICATION TO MULTI-AGENT SYSTEMS ÁRON FEHÉR *1,2 AND LŐRINC MÁRTON1 1Department of Electrical Engineering, Sapientia Hungarian University of Transylvania, C.1., Târgu Mureş, 547367, ROMANIA 2Department of Mathematics, University of Pannonia, Egyetem u. 10, Veszprém, 8200, HUNGARY This paper presents a review of dominant pole and model-approximation algorithms for delayed systems that can be applied to multi-agent systems. A novel algorithm is proposed to determine an approximation method for multi-agent systems in the platoon configuration with a communication delay. Simulations are presented to show the applicability of the proposed algorithm. Keywords: time delay, differential equations, asymptotic properties, multi-agent systems 1. Introduction A Multi-Agent System (MAS) consists of multiple ac- tive agents (e.g. vehicles), passive agents (e.g. obstacles), in addition to cognitive agents and their environment [1]. Every active agent is at least partially autonomous and uses a distributed control algorithm [2]. The agents can communicate with each other. This communication struc- ture is defined by a graph, where each vertex corresponds to one agent and each edge to a communication direction. A platoon is a special MAS configuration, with a linear communication graph. The leading (first) agent implements a reference tracking algorithm. Every other agent implements a consensus with the adjacent agents [3]. With the increase in the number of agents and the physical distance between the agents, the communication delay cannot be neglected. While the stability of such sys- tems is ensured by the consensus protocol [4], the delay will influence the transient behavior of the MAS [5]. The goal of this paper is to compare the existing ap- proximation methods for the transient behavior analysis of MAS with communication delays and present a novel analysis method which can be applied to any MAS sys- tem which satisfies a smallness delay condition. 2. Modelling of MAS A MAS is considered with agents that exhibit single- integrator dynamics [6]. The state-space model of an agent becomes ẋi(t) = ui(t), where xi ∈ R denotes *Correspondence: fehera@ms.sapientia.ro the state of the ith agent and ui ∈ R represents the input, i = 1, 2, . . . ,n. A MAS has an underlying communication graph, in which the vertex is an agent and the edge a communica- tion path [7], so the agent ith in the system is the vertex vi. Let Ni be the set of neighbors of vi, so that Ni con- tains all vertices that are connected to vi. Consensus algorithm The consensus problem of a MAS is the procedure of gathering every state from the initial condition to a common steady-state. If the commu- nication graph is connected, the consensus with regard to an agent can be reached with the input (consensus proto- col) ui(t) = ∑ j∈Ni (xj(t) −xi(t)) . (1) The adjacency matrix A = (aij) ∈ Rn×n of a graph with n nodes is defined as aij := { 1, if i 6= j and vi are adjacent to vj 0, otherwise . (2) The degree matrix D = (dij) ∈ Rn×n of a graph with n nodes shows the number of neighbors for each vertex and can be defined as dij := { deg(vi), if i = j 0, otherwise , (3) where deg(vi) denotes the degree or the number of edges incident to vertex i. mailto:fehera@ms.sapientia.ro 72 FEHÉR AND MÁRTON Figure 1: The communication topology of a vehicle pla- toon system consisting of n vehicles with time delay τ in the communication graph. The dashed line symbolizes more nodes in between, while the dotted line represents the reference input of the nodes. With this notation the dynamics of the MAS with the consensus protocol is given by ẋ(t) = −Lx(t), x(0) = x0, (4) where L denotes the Laplacian matrix [8] which is con- structed as L = D − A, D stands for the degree ma- trix, A represents the adjacency matrix of the graph, x = ( x1 x2 . . . xn )> ∈ Rn denotes the state vec- tor consisting of the n states of the MAS, and x0 ∈ Rn is a constant vector of the initial states. According to [9] the eigenvalues of a MAS consisting of n agents with a connected communication graph can be ordered as 0 = λ1 > λ2 ≥ ···≥ λn. (5) The steady states (equilibria) of the MAS xss are the el- ements of the null space of L. Given, according to the definition of the Laplacian matrix, ∑ j∈Ni lij = 0 [6] for every solution x of (Eq. 4): lim t→∞ x(t) = xss = 1 n ∑n i=1 xi(0)1, where 1 = ( 1 1 . . . 1 )> ∈ Rn. MAS with delayed communication In the case of de- layed communication, the consensus protocol is of the form ui(t) = ∑ j∈Ni (xj(t− τ) −xi(t)) , (6) where τ ≥ 0 denotes the constant communication delay which is present among adjacent agents. According to Refs. [4] and [10], the MAS with com- munication delay (also referred to as Multi-Agent System with Delays (DMAS)) is ẋ(t) = −Dx(t) + Ax(t− τ) + R, (7) x(θ) = x0 ∈ R n, θ ∈ [−τ, 0], 3. Platoon of vehicles A platoon of vehicles is a special class of MAS with a communication topology shown in Fig. 1. The kinematic model of a vehicle can be written as ẋi(t) = ui(t), (8) where xi(t) denotes the position of the vehicle and ui(t) represents the control signal in the form{ u1(t)=kp1 (r−x1(t)) ui(t)=ξi(xi−1(t−τ)−d−xi(t)),∀i = 2, ...,n (9) where i ∈ (1,n), n ∈ N,n > 1, d denotes the prescribed inter-vehicle distance, ξi > 0, and kpi > 0 are constant control gains. r is the position reference of the first vehi- cle. Eqs. 8 and 9 can be written as a system of differential equations with delay: ẋ(t) = −Φx(t) + Γx(t− τ) + f 1 r −f 2 d, (10) where Φ = diag{ ( kp1 ξ2 . . . ξn−1 ξn ) } (11) denotes the degree matrix with the state feedback, Γ =   0 0 0 . . . 0 0 0 ξ2 0 0 . . . 0 0 0 ... ... ... . . . ... ... ... 0 0 0 . . . ξn−1 0 0 0 0 0 . . . 0 ξn 0   (12) represents the adjacency matrix, and f 1 = ( kp1 0 . . . 0 )> , (13) f 2 = ( 0 ξ2 . . . (n− 1)ξn )> . (14) The homogeneous part of the relation in Eq. 10 is of the same form as the relation in Eq. 7, and the term f 1 r−f 2 d represents the reference as well as inter-vehicle distance induced inflows. 4. Existing methods for the approximation of time-delay systems In this section, the various current approximation meth- ods for time-delay systems are reviewed. The form of the studied systems is shown by the following relation: ẋ(t) = −Φx(t) + Γx(t− τ),x(h) = x0, (15) which is the homogeneous part of Eq. 10, where θ(h) denotes the initial condition with h ∈ [−τ, 0], and has a quasi-polynomial characteristic equation: λIn + Φ − Γe−τλ = 0, (16) where the delay component induces an exponential term. The roots of Eq. 16 determine the transient behaviour of the system in Eq. 15. As Eq. 15 possesses an infinite number of solutions, it is important to develop such an equivalent system which has a finite number of eigenval- ues and is a good approximation of the original system in Eq. 15. If a good delay-free approximation is available, it Hungarian Journal of Industry and Chemistry APPROXIMATION OF MULTI-AGENT SYSTEMS WITH DELAYS 73 could be applied to the transient behaviour analysis and control design for delayed systems. Let a special solution be denoted by x̃, which is uniquely determined by the value x̃(0), and independent of θ(h), h ∈ [−τ, 0), thus forming an n parameter fam- ily. In a linear autonomous system, this corresponds to the eigensolution generated by exactly n characteristic roots (multiplicities included) which lies in the half-plane Re λ > −1/τ, see Ref. [11]. Theorem 1. [11] Consider a Delay Differential Equation (DDE) system of the form shown in relation Eq. 15 with a Lipschitz criterion of: (‖Φ‖ + ‖Γ‖)τe < 1, (17) further noted as smallness condition, for every solution x of Eq. 15 a globally defined solution x̃ : R → Rn exists that satisfies the growth condition supt≤0 ‖x̃(t)‖et/τ < ∞ such that ‖x(t) − x̃(t)‖→ 0 exponentially as t →∞. For further related discussions, see Refs. [12] and [13]. The aforementioned theorem yields n dominant eigenval- ues which can accurately represent a DDE system. The system shown in the relation in Eq. 15 uses the consensus protocol, which ensures that the eigenvalues are located in the left half-plane in the complex region. The final lo- cation of the dominant eigenvalues is created as a semi- circle with origin 0, radius 1/τ, and a negative real part. 4.1 The modified chain approximation The modified chain approximation method creates an ap- proximating system directly from the state-space repre- sentation of the delayed system [14]. For a DDE in the form of Eq. 15, the modified chain approximation method yields an approximating linear system in the form of ẏ 0 (t) = −Φy 0 (t) + m τ Inym(t) ẏ 1 (t) = Γy 0 (t) − m τ Iny1(t) (18) ... ẏ k (t) = m τ Inyk−1(t) − m τ Inyk(t), 2 ≤ k ≤ m The output vector z = y 0 represents the approximation of the solution x of the relation in Eq. 15. The approxi- mating system in the form of a matrix is shown in ẏ(t) = Gy(t) (19) with G =   −Φ 0n 0n · · · 0n mτ In Γ −m τ In 0n · · · 0n 0n 0n m τ In −mτ In · · · 0n 0n ... ... ... . . . ... ... 0n 0n 0n · · · mτ In − m τ In   (20) where y ∈ Rmn denotes the state vector, G ∈ R(m+1)n×(m+1)n represents the matrices of the system, 0n ∈ Rn×n stands for the zero matrix, n is the number of agents and m denotes the number of approximating equations. According to Ref. [14] the output of the system (Eq. 19) defined as z(t) = y 0 (t) linearly converges into the solution of the original DDE system such that c > 0 and supt≥0‖x(t) −z(t)‖≤ c m . 4.2 The Lambert W function The Lambert W function can be used to find the domi- nant eigenvalues of a quasi-polynomial equation (Eq. 16). Every W(s) function that satisfies W(s)eW(s) = s (21) by definition is referred to as a Lambert W function [15], where s is either a scalar or matrix complex number func- tion. The Lambert W function has multiple branches de- noted as Wk(s) with k = 0,±1,±2, . . . ,±∞ . Example 4.1. If a scalar DDE is present in the form of ẋ(t) = ax(t) + bx(t− τ) + cu(t), (22) with a,b,c,θ ∈ R, x(h) = θ(h) for h ∈ [−τ, 0], the quasi-polynomial of the homogeneous part can be written as (λ−a)eτλ = b. (23) If both sides are multiplied by τe−aτ , then (λ−a)τeτ(λ−a) = bτe−aτ, (24) which satisfies Eq. 21 with W(bτe−aτ ) = (λ−a)τ, and the eigenvalues can be calculated using the branches of the Lambert W function as λk = 1 τ Wk(bτe −aτ ) + a. (25) In our case, in terms of the relation in Eq. 15, the aforementioned solution is generalized as λk = 1 τ Wk(ΓτQk) − Φ, (26) 47(1) pp. 71–77 (2019) 74 FEHÉR AND MÁRTON where Qk can be calculated by solving the equation nu- merically Wk(ΓτQk)e Wk(ΓτQk)−Φτ = Γτ (27) for Qk [16]. Although many numerical solvers provide native sup- port for the solution of the scalar Lambert W function, the general case requires additional solver tools. There- fore, the LambertWDDE Toolbox [17] was created. The function find_Sk assumes τ, Γ and −Φ. The returned val- ues are the eigenvalues λk and the Qk parameters for a given k branch. The toolbox can create the approximat- ing solution for the given system as x̃(t) = m2∑ k=m1 eλktCIk, (28) where x̃ is the approximating solution of the original de- lay system and the parameter CIk can be computed with the help of find_CI for a given m1 < k < m2 branch. In the scalar case, CIk takes the form of CIk = x0 + be −λkτ ∫ τ 0 θ(t− τ)dt 1 + bτe−λkτ . (29) 4.3 The Quasi-polynomial root-finder algo- rithm The quasi-polynomial root-finder algorithm calculates the dominant eigenvalues of a system based directly on the quasi-polynomial equation in Eq. 16. The quasi-polynomial equation of the system in Eq. 16 can be written as: P(λ) = N∑ k=0 Qk(λ)e −αkτλ, (30) where Qk is a polynomial with real coefficients and αk ∈ R. The objective is to compute the spectrum in the region of the complex plane D ⊂ C with boundaries βmin < Re(D) < βmax and ωmin < Im(D) < ωmax. Let the surfaces defined by the real and imaginary parts of P(λ) be: Re(P(β,ω)) = 0 (31) Im(P(β,ω)) = 0 (32) The eigenvalues can be located at the points of in- tersection of the zero-level curves of the surfaces Re(P(β,ω)) = 0 and Im(P(β,ω)) = 0 as shown in Ref. [18]. The accuracy of the algorithm is increased by Newton’s method and by adapting the grid density of D as shown in Ref. [16]. The Quasi-polynomial root-finder algorithm (QPmR) [19] is implemented in MATLAB. The function expects the region of interest [βmin,βmax,ωmin,ωmax] to be in the complex plane of the polynomial coefficient matrix of the quasi-polynomial where one row corresponds to one polynomial multiplied by the same exponential term. The delay vector, computational accuracy and grid step are also required. In our case this translates into a region of interest [−1 τ , 0,−1 τ , 1 τ ] if the smallness condition (Eq. 17) is sat- isfied. The first row of the matrix of polynomial coeffi- cients contains the coefficients of the delay-free part so that the delay vector is of the form [0,τ, 2τ, ...]. The algorithm covers the given region with a mesh grid, then evaluates the quasi-polynomial at each point of the grid by splitting it into a real and an imaginary part. The zero-level curves are then mapped with the help of the contour plotting algorithm. The computational error is checked and if it is too large, the algorithm is restarted using a modified grid density as described in Ref. [16]. If the computational error is smaller than the given level of tolerance, the computed dominant eigenvalues are re- turned. 5. Explicit matrix approximation method An approximation method was devised where the con- vergence rate is exponential, the degree of the resulting system in the form of Eq. 4 matches exactly the degree of the delayed MAS given in Eq. 15, and the same properties are exhibited in specific cases. The Banach fixed-point theorem was used as dis- cussed in Ref. [20] to explicitly find a linear system of the form of Eq. 4 which approximates the homogeneous part of the system in Eq. 10. If an (X,f) metric space is present and T : B → B is a contraction with a bounded set B ⊂ X, and q < 1 such that f(T(x),T(y)) ≤ qf(x,y) (33) by definition T admits a unique fixed point x̃ such as T(x̃) = x̃, and this fixed point can be found by starting from an arbitrary element x0 ∈ B with the sequence xn = T(xn−1), (34) where xn → x̃. A DDE system is defined in Eq. 15 by the corre- sponding smallness condition of Eq. 16. Since all normed spaces are metric spaces, the metric space X = Rn×n is set with f as the induced matrix norm. The contraction T : B → B is present such that T(Λ) = −Φ + Γe−Λτ, (35) and B = {Λ ∈ Rn×n|‖Λ‖ ≤ (‖Φ‖ + ‖Γ‖)e}. The rela- tion in Eq. 33 holds true for (‖Φ‖ + ‖Γ‖)τe < 1. Let Λ1, Λ2 ∈ B such that ‖Λ1‖ > ‖Λ2‖. The left side of the inequality in Eq. 33 can be written as ‖T(Λ1) −T(Λ2)‖ = ‖Γ‖‖e−Λ1τ − e−Λ1τ‖. It is evident that ‖Γ‖ ≤ ‖Γ‖ + ‖Φ‖ and the maximum norm can be used as ‖e−Λ1τ − e−Λ1τ‖≤ τ‖Λ1 − Λ2‖eτ max{‖Λ1‖,‖Λ2|}. Hungarian Journal of Industry and Chemistry APPROXIMATION OF MULTI-AGENT SYSTEMS WITH DELAYS 75 It can be seen that ‖T(Λ1)−T(Λ2)‖≤ τ‖Λ1−Λ2‖(‖Γ‖+‖Φ‖)eτ(‖Γ‖+‖Φ‖)e and by applying the smallness condition ‖Γ‖ + ‖Φ‖ ≤ 1/(τe) (Eq. 17), ‖T(Λ1) −T(Λ2)‖≤‖Λ1 − Λ2‖ is obtained, which proves that T : B → B is a contrac- tion. This shows that by solving Λ = −Φ + Γe−Λτ (36) for the Λ matrix, a system is created dx̃ dt = Λx̃, (37) which approximates the DDE system of Eq. 15. As a re- sult of the proposed iterative method, the eigenvalues and eigenvectors of Eq. 37 approximate, with a given degree of precision, the dominant eigenvalues of Eq. 14. 5.1 Comparison with the existing methods • Chain approximation: + The result is an approximating system with known system matrices (both eigenvalues and eigenvectors are known). - The resulting system is of a higher degree than the delay system. - The convergence rate of the algorithm is linear. • Lambert W function: + The result is a trajectory approximation. + The convergence rate of the algorithm is expo- nential. - The algorithm requires numerical solvers for an exponential matrix equation. - Multiple branches of the Lambert W function must be used to create an accurate approxima- tion, and the number of eigenvalues in a branch cannot be predetermined generally. • QPmR algorithm: + The algorithm determines the exact number of eigenvalues in the given complex domain, with the given computational error. + The convergence rate of the algorithm is expo- nential. - The algorithm uses the quasi-polynomial equa- tion, thus, will not contain any information on the eigenvectors. - The algorithm uses numerical solvers to com- pute the zeros of the zero-level curves created from the quasi-polynomial equation. • Approximation of an explicit matrix: + The result is an approximating system. + The resulting system is of the same degree as the approximating system. - The algorithm calculates a matrix exponen- tial numerically, which is a compute-intensive task. 6. Simulations and results Let us consider two cases: a MAS consisting of five and twenty-five agents, respectively, with first-order dynam- ics in a platoon configuration as shown in the relation of Eq. 10, with the constant initial condition θ(h) = x0. For the comparisons, a 6th order chain approxima- tion was used. In the case of the Lambert W function, the initial matrix Q0 = ( 1 1 1 1 ) and the branches k = −2,−1, 0, 1 were used. For the quasi-polynomial root- finder algorithm (QPmR), a symbolic calculation to find the characteristic quasi-polynomial equation of the sys- tem was used, and the plane of the search was set to [−τ, 0]× [−τj,τj]. The algorithm for explicit matrix ap- proximation was used with the initial matrix Λ0 = 0n×n. The error threshold 1e−7 was used in every iterative al- gorithm. The dominant eigenvalues of the system consisting of five agents, with minor differences, was identified by ev- ery approximation method. In the case of the larger sys- tem, the chain and explicit matrix approximations could generate a result, while the quasi-polynomial root-finder algorithm and the Lambert W function were determined by numerical calculations. As such, the smaller system was chosen as a point of comparison for the algorithms. Tables 1 and 2 contain a comparative summary of the four aforementioned algorithms: the number of iterations, the overall computation time and the dominant eigenval- ues identified for a platoon consisting of five agents as shown in Fig. 2. Fig. 3 shows that the linear system is generated by the explicit matrix approximation in just twelve steps for the DMAS that consists of twenty-five agents. Fig. 4 shows that the resultant approximating system exhibits the same steady-state and transient behavior as the original delayed system using the same initial condi- tions. Since the initial position falls within the range of Table 1: The number of iterations and the overall compu- tation time for the algorithms. Name of algorithm No. of cycles Computation time Chain approximation 1 0.02 s Lambert W function 35 7.5 min QPmR algorithm 12 17.18 s Explicit algorithm 4 0.03 s 47(1) pp. 71–77 (2019) 76 FEHÉR AND MÁRTON Table 2: The eigenvalues returned by the studied algo- rithms. Name of algorithm Eigenvalues Chain approximation -0.1080, -0.9471, -1, -2.4736, -4.2511 Lambert W function -0.1081, -0.9482, -1, -2.4658, -4.1603 QPmR algorithm -0.1081, -0.9481, -1, -2.4661, -4.1603 Explicit algorithm -0.1080, -0.9481, -1, -2.4661, -4.1603 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 Re( ) -1 -0.5 0 0.5 1 Im ( ) Figure 2: Dominant poles of the DMAS consisting of five agents. 0 20 40 60 80 100 120 Number of iterations (1) -60 -50 -40 -30 -20 -10 0 10 20 30 tr ( n ) tr( n ) |tr(| n )-tr( n-1 )|=1e-6 Figure 3: The explicit matrix approximation returns a valid system after 12 iterations for 25 agents. 0 to 100 m, the error of the approximating system is 8 cm in the transient domain and 5 mm under steady-state conditions, as shown in Fig. 5. 7. Conclusions An iterative algorithm was proposed and tested based on which a degree-preserving approximation model can be created for a class of MAS in a platoon formation consist- ing of agents that exhibit first-order dynamics in the pres- ence of a small communication delay. The algorithm uses methods of numerical computation. The obtained MAS is of the same degree and steady state as the delayed MAS, moreover, it correctly approximates the transient behav- ior. The algorithm was compared with existing methods for approximating eigenvalues and systems. Simulations show that the presented algorithm is suitable for the anal- ysis of complex MAS. 0 50 100 150 200 250 300 350 400 450 500 Time (s) 0 10 20 30 40 50 60 70 80 90 D is ta n c e ( m ) Figure 4: The trajectories of the platoon of the delayed MAS compared to the platoon of the approximated MAS created by the explicit matrix approximation. 0 10 20 30 40 50 60 70 80 90 100 Time (s) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 D is ta n c e ( m ) Figure 5: The trajectory of the approximation error. Acknowledgement The authors would like to express their gratitude to Dr. Mihály Pituk (University of Pannonia, Hungary) for his useful comments. REFERENCES [1] Kubera, Y.; Mathieu, P.; Picault, S.: Everything can be Agent!, in Proceedings of the ninth Interna- tional Joint Conference on Autonomous Agents and Multi-Agent Systems (AAMAS’2010) (W. van der Hoek, G.A. Kaminka, Y. Lespérance, M. Luck, S. Sen, eds.) (International Foundation for Au- tonomous Agents and Multiagent Systems, Toronto, Ontario, Canada), 1547–1548 [2] Panait, L.; Luke, S.: Cooperative multi-agent learn- ing: The state of the art, Auton. Agents Multi-Agent Syst., 2005 11(3), 387–434, DOI: 10.1007/s10458-005- 2631-2, ISBN: 978-981-10-2491-7 [3] Zabat, M.; Stabile, N.; Farascaroli, S.; Browand, F.: The aerodynamic performance of platoons: A final report, Institute of Transportation Studies, Research Reports, Working Papers, Proceedings, Institute of Transportation Studies, UC Berkeley, 1995 [4] Cheng-Lin, L.; Fei, L.: Consensus Problem of Delayed Linear Multi-agent Systems, Springer- Briefs in Electrical and Computer Engineering (Springer Singapore), 2017, DOI: 10.1007/978-981-10- 2492-4, ISBN:978-981-10-2491-7 Hungarian Journal of Industry and Chemistry https://doi.org/10.1007/s10458-005-2631-2 https://doi.org/10.1007/s10458-005-2631-2 https://doi.org/10.1007/978-981-10-2492-4 https://doi.org/10.1007/978-981-10-2492-4 APPROXIMATION OF MULTI-AGENT SYSTEMS WITH DELAYS 77 [5] Wim, M.; Silviu-Iulian, N.: Stability and stabiliza- tion of time-delay systems: An eigenvalue-based approach (SIAM), 2007, DOI: 10.1137/1.9780898718645 [6] Lewis, F.L.; Zhang, H.; Hengster-Movric, K.; Das, A.: Cooperative control of multi-agent systems: Op- timal and adaptive design approaches (Springer), 2014, DOI: 10.1007/978-1-4471-5574-4 [7] Trudeau, R.J.: Introduction to Graph Theory, Dover Books on Mathematics (Dover Publications), 2nd edn., 1994, ISBN: 978-0486678702 [8] Chaiken, S.; Kleitman, D.J.: Matrix tree theorems, J. Comb. Theory A, 1978 24(3), 377 – 381, DOI: 10.1016/0097-3165(78)90067-5 [9] Mesbahi, M.; Egerstedt, M.: Graph theoretic meth- ods in multiagent networks, Princeton Series in Applied Mathematics (Princeton University Press), 2010, DOI: 10.1515/9781400835355 [10] Cepeda-Gomez, R.; Olgac, N.: An exact method for the stability analysis of linear consensus proto- cols with time delay, IEEE T. Autom. Control, 2011 56(7), 1734–1740, DOI: 10.1109/TAC.2011.2152510 [11] Arino, O.; Pituk, M.: More on linear differ- ential systems with small delays, J. Differen- tial Equations, 2001 170(2), 381 – 407, DOI: 10.1006/jdeq.2000.3824 [12] Győri, I.; Pituk, M.: Asymptotically ordinary de- lay differential equations, Functional Differential Equations, 2005 12, 187–208 [13] Győri, I.; Pituk, M.: Asymptotic formulas for a scalar linear delay differential equation, Electron. J. Qual. Theory Differ. Equ., 2016 1(72), 1–14, DOI: 10.14232/ejqtde.2016.1.72 [14] Krasznai, B.; Győri, I.; Pituk, M.: The modi- fied chain method for a class of delay differen- tial equations arising in neural networks, Math. Comput. Model., 2010 51(5), 452 – 460, DOI: 10.1016/j.mcm.2009.12.001 [15] Corless, R.M.; Gonnet, G.H.; Hare, D.E.G.; Jef- frey, D.J.; Knuth, D.E.: On the Lambert W func- tion, Adv. Comput. Math., 1996 5(1), 329–359, DOI: 10.1007/BF02124750 [16] Vyhlídal, T.; Lafay, J.F.; Sipahi, R.: Delay sys- tems: From theory to numerics and applications, Advances in Delays and Dynamics (Springer In- ternational Publishing), 2013, DOI: 10.1007/978-3-319- 01695-5 [17] Yi, S.; Duan, S.; Nelson, P.W.; Ulsoy, A.G.: The Lambert W Function Approach to Time Delay Systems and the LambertW_DDE Toolbox, IFAC Proceedings Volumes, 2012 45(14), 114–119, DOI: 10.3182/20120622-3-US-4021.00008 [18] Vyhlídal, T.; Zítek, P.: Quasipolynomial mapping based rootfinder for analysis of Time delay systems, in Proc. of IFAC Workshop on Time-Delay Sys- tems, TDS 2003, DOI: 10.1016/S1474-6670(17)33330-X [19] Vyhlídal, T.; Zítek, P.: QPmR – Quasipolynomial rootfinder: Algorithm update and examples, Ad- vances in Delays and Dynamics, vol. 1 (Springer International Publishing, Cham), 2013 299–312, DOI: 10.1007/978-3-319-01695-5_22, http://www.cak.fs. cvut.cz/algorithms/qpmr [20] Latif, A.: Banach contraction principle and its gen- eralizations, Topics in Fixed Point Theory (Springer International Publishing, Cham), 2014 33–64, DOI: 10.1007/978-3-319-01586-6_2 47(1) pp. 71–77 (2019) https://doi.org/10.1137/1.9780898718645 https://doi.org/10.1007/978-1-4471-5574-4 https://doi.org/10.1016/0097-3165(78)90067-5 https://doi.org/10.1016/0097-3165(78)90067-5 https://doi.org/10.1515/9781400835355 https://doi.org/10.1109/TAC.2011.2152510 https://doi.org/10.1006/jdeq.2000.3824 https://doi.org/10.1006/jdeq.2000.3824 https://doi.org/10.14232/ejqtde.2016.1.72 https://doi.org/10.14232/ejqtde.2016.1.72 https://doi.org/10.1016/j.mcm.2009.12.001 https://doi.org/10.1016/j.mcm.2009.12.001 https://doi.org/10.1007/BF02124750 https://doi.org/10.1007/BF02124750 https://doi.org/10.1007/978-3-319-01695-5 https://doi.org/10.1007/978-3-319-01695-5 https://doi.org/10.3182/20120622-3-US-4021.00008 https://doi.org/10.3182/20120622-3-US-4021.00008 https://doi.org/10.1016/S1474-6670(17)33330-X https://doi.org/10.1007/978-3-319-01695-5_22 http://www.cak.fs.cvut.cz/algorithms/qpmr http://www.cak.fs.cvut.cz/algorithms/qpmr https://doi.org/10.1007/978-3-319-01586-6_2 https://doi.org/10.1007/978-3-319-01586-6_2 Introduction Modelling of mas Platoon of vehicles Existing methods for the approximation of time-delay systems The modified chain approximation The Lambert W function The Quasi-polynomial root-finder algorithm Explicit matrix approximation method Comparison with the existing methods Simulations and results Conclusions