Solution of Systems of Linear Delay Differential Equations via Laplace Transformation Sun Yi, A. Galip Ulsoy Senior Member, IEEE and Patrick W. Nelson Abstract— An approach using the Lambert W function for the analytical solution, free and forced, to systems of delay differential equations with a single delay has been developed in [8], [9]. The solution is expressed in the form of an infinite series of modes written in terms of the matrix Lambert W function. In this paper, we utilize the analytical solution to present a solution in the Laplace domain, present validation examples, and emphasize the analogy of the solution method to systems of ordinary differential equations.
I. INTRODUCTION Time-delay systems (TDS) are systems in which a significant time delay exists between the applications of input to the system and their resulting effect. Such systems arise from an inherent time delay in the components of the system or a deliberate introduction of time delay into the system for control purposes. TDS can be represented by delay differential equations (DDEs) which belong to the class of functional differential equations (FDEs), and have been extensively studied over the past decades [1]. DDEs, also known as difference differential equations, were initially introduced in the 18th century by Laplace and Condorcet. The basic theory concerning stability of systems described by equations of this type was developed by Pontryagin in 1942. Important works have been written by Bellman and Cooke in 1963, Smith in 1957, Pinney in 1958, Halanay in 1966, El’sgol’ts and Norkin in 1971, Myshkis in 1972, Yanushevski in 1978, Marshal in 1979, and Hale in 1977. The reader is referred to the detailed review in [2]. The principal difficulty in studying DDEs lies in its special transcendental character. One of the well-known approximation methods is the Pad´e approximation, which results in a shortened repeating fraction for the approximation of the characteristic equation of the delay [3]. Another approach to the delay problem is to look at the entire delay spectrum. DDEs are often solved using numerical methods, asymptotic solutions, and graphical tools. A related study on analytic solutions of linear DDEs can be found in Falbo, 1995 [4]. A Fourier-like analysis of the existence of the solution and its properties for the nonlinear DDEs is studied by Wright, 1946 [5]. Similar approaches to linear and nonlinear DDEs are also reported by Bellman in 1963 [6]. The uniqueness of the solution and its properties for the linear DDEs are also studied by Wright This research was supported by the National Science Foundation (NSF Grant No. 0555765). Sun Yi and A. G. Ulsoy are with the Mechanical Engineering Department, University of Michigan, Ann Arbor, MI 48109-2125 USA
[email protected]. Patrick W. Nelson is with the Department of Mathematics, University of Michigan, Ann Arbor, MI 48109.
in 1948 [7]. An analytic approach to obtain the complete solution of systems of DDEs based on the concept of the matrix Lambert W function was developed by Asl and Ulsoy in 2003 [8]. And their analytical approach has been extended to general systems of DDEs and non-homogeneous DDEs [9]. The advantage of this approach lies in the fact that the form of the solution obtained is analogous to the general solution form of ordinary differential equations (ODEs), thus the concept of the state transition matrix in ODEs can be generalized to DDEs (see Table I). The Laplace transform is known to be useful in obtaining solutions of linear DDEs with constant coefficients [6]. In this paper, the Laplace transform approach for solving linear systems of DDEs is developed using the matrix Lambert W function method. This approach makes it possible to solve DDEs analytically, contributes a clear understanding of the properties of the delayed systems, and is expected to enable the analysis of systems of DDEs for controllability, pole placement, etc. II. S OLUTION TO S YSTEMS OF DDE S U SING L AMBERT W FUNCTION
THE
In this section we briefly summarize some key results from [8], [9]. First, for the first-order scalar DDE, x(t) ˙ + ax(t) + ad x(t − h) = bu(t), t > 0 x(t) = g(t), t ∈ [−h, 0) x(t) = x0 , t=0
(1)
the free solution in terms of the Lambert W function, Wk , is [8] ∞ X
1 Wk (−ad heah ) − a h k=−∞ (2) where the coefficient CkI is determined from the preshape function, g(t) and the initial condition, x0 . Every function W (h), such that W (h)eW (h) = h, is called a Lambert W function. The Lambert function, W (h), is complex valued, with a complex argument h, and has an infinite number of branches Wk (h), where k = −∞, −1, 0, 1, . . . , ∞ [12]. The forced solution is Z t X ∞ x(t) = eSk (t−ξ) CkN bu(ξ)dξ (3) x(t) =
eSk t CkI , where Sk =
0 k=−∞
where CkN in (3) is a constant and the method for its evaluation is presented in [9]. Hence, the total solution to
TABLE I C OMPARISON OF THE S OLUTIONS TO ODE S AND DDE S ODEs
DDEs
Scalar Case x(t) ˙ + ax(t) = bu(t), t > 0 x(t) = x0 , t = 0 Z t x(t) = e−at x0 + e−a(t−ξ) bu(ξ)dξ
x(t) ˙ + ax(t) + ad x(t − h) = bu(t), t > 0 x(t) = g(t), t ∈ [−h, 0); x(t) = x0 , t = 0 Z t X ∞ ∞ X x(t) = eSk t CkI + eSk (t−ξ) CkN bu(ξ)dξ
0
0 k=−∞
k=−∞
where, Sk = Matrix-Vector Case x˙ (t) + Ax(t) = Bu(t), t > 0 x(t) = x0 , t = 0 Z t x(t) = e−At x0 + e−A(t−ξ) Bu(ξ)dξ 0
˙ x(t) + Ax(t) + Ad x(t − h) = Bu(t), t > 0 x(t) = g(t), t ∈ [−h, 0); x(t) = x0 , t = 0 Z t X ∞ ∞ X x(t) = eSk t CIk + eSk (t−ξ) CN k Bu(ξ)dξ 0 k=−∞
k=−∞
where, Sk =
∞ X
x(t) =
eSk t CkI + {z
|
k=−∞
|
}
free
t
∞ X
0 k=−∞
eSk (t−ξ) CkN bu(ξ)dξ {z
forced
}
(4) The matrix Lambert W function approach can be applied to solve systems of DDEs in matrix-vector form [9] ˙ + Ax(t) + Ad x(t − h) = Bu(t) t > 0 x(t) x(t) = g(t) t ∈ [−h, 0) x(t) = x0 t=0
1 Wk (−Ad hQk ) − A h
ˆ ˆ ˆ Hk = Zk Jk Z−1 k . Jk = diag(Jk1 (λ1 ), Jk2 (λ2 ), . . . , Jkp (λp )), ˆ where Jki (λi ) is m × m Jordan block and m is multiplicity ˆ i . Then, the matrix Lambert W function of the eigenvalue λ can be computed as [10]
the system in (1) becomes Z
1 Wk (−ad heah ) − a h
(5)
where A and Ad are n × n matrices, x(t) is an n × 1 state vector, B is an n × r matrix, and u(t), an r × 1 vector, is a function representing the external excitation. Assuming a free solution form as
Wk (Hk ) = ˆ 1 )), . . . , Wk (Jkp (λ ˆ p ))}Z−1 Zk {diagWk (Jk1 (λ k
(11)
where ˆi) = Wk (Jki (λ ˆ i ) W ′ (λ ˆ Wk (λ k i) · · · ˆi ) · · · 0 Wk (λ .. .. .. . . . 0 0 ···
(m−1) ˆ 1 (λi ) (m−1)! Wk
0 .. . ˆi) Wk (λ
(12)
(7)
CIk is determined from the given initial conditions, x0 , and the preshape function, g(t), as in the scalar case. The forced solution is Z t X ∞ x(t) = eSk (t−ξ) CN (13) k Bu(ξ)dξ
1 Wk (−Ad hQk ) − A (8) h And the following condition is used to solve for the unknown matrix Qk .
where CN k is n×n coefficient matrices, and the total solution is Z t X ∞ ∞ X Sk t I x(t) = e Ck + eSk (t−ξ) CN k Bu(ξ)dξ
x(t) = eSt x0
(6)
the free solution to (5) is derived as x(t) =
∞ X
eSk t CIk
k=−∞
where
Sk =
Wk (−Ad hQk )e Wk denotes the k which satisfies
th
Wk (−Ad hQk )−Ah
= −Ad h
(9)
branch of the matrix Lambert W function
Wk (Hk )eWk (Hk ) = Hk
0 k=−∞
0 k=−∞
k=−∞
|
{z
free
}
|
{z
forced
}
(14) The solution to DDEs in terms of the matrix Lambert W function, and its analogy to that of ODEs are summarized in Table I.
(10)
Corresponding to each branch, k, of the Lambert function, Wk , there is a solution Qk from (9), and for Hk = Ad hQk , we compute the Jordan canonical form Jk from
III. L APLACE T RANSFORM
OF
S CALAR DDE S
Consider the scalar free DDE in (1) without external force u(t). Because x(t) generally has a non-zero preshape
function, it can be transformed as, Z ∞ e−st x(t − h)dt = 0 Z h Z ∞ = e−st x(t − h)dt + e−st x(t − h)dt 0 h Z h Z ∞ −st = e g(t)dt + e−s(t+h) x(t)dt 0
According to (21), nk (s) has the interesting and useful property,
(15)
0
= G(s) + e−sh X(s) Then, the transform of the free equation is sX(s) − x0 + ad e−sh X(s) + ad G(s) + aX(s) = (s + ad e−sh + a)X(s) − x0 + ad G(s) = 0
(16)
Therefore, X(s) =
x0 − ad G(s) s + ad e−sh + a
(17)
On the other hand, the solution obtained by the approach using the Lambert W function in (2) can be transformed as I C−1 C0I C1I + + + ··· s − S−1 s − S0 s − S1 ∞ X CkI = s − Sk k=−∞ (18) Sk is defined in (2). Writing X(s) in (18) over a common denominator yields
X(s) = · · · +
∞ X CkI nk (s) X(s) = d(s)
nk (s = Sl ) = 0, when k 6= l = · · · (Sl − Sk−2 )(Sl − Sk−1 )(Sl − Sk+1 )(Sl − Sk+2 ) · · · when k = l (24) With this property, we can compute the residue CkI from (23) by substituting Sk for s, that is, .. . J(S ){x − ad G(S0 )} 0 0 C0I = n0 (S0 ) J(S0 ){x0 − ad G(S0 )} = · · · (S0 − S−2 )(S0 − S−1 )(S0 − S1 ) · · · J(S1 ){x0 − ad G(S1 )} C1I = n1 (S1 ) J(S1 ){x0 − ad G(S1 )} = · · · (S1 − S−1 )(S1 − S0 )(S1 − S2 ) · · · .. . where J(s) is calculated from (22) as ∞ Y
(s − Sk ) = J(s)(s + ad e−sh + a)
k=−∞
⇐⇒ J(s) = (19)
where, d(s) and nk (s) are polynomials in the Laplace variables, s ∞ Y
J(s) = (s − Sk )
k=−∞
lim
s→S0
(20)
= · · · (s − S−1 )(s − S0 )(s − S1 ) · · · d(s) nk (s) ≡ (s − Sk ) = · · · (s − Sk−2 )(s − Sk−1 )(s − Sk+1 )(s − Sk+2 ) · · · (21) Then, comparing denominators of both sides in (19) and (17) yields d(s) =
∞ Y
(s − Sk ) = J(s)(s + ad e−sh + a)
(s − Sk )
(26)
k=−∞
(s + ad e−sh + a)
=
lim
∞ Y
(s − Sk )
k=−∞
(s + ad e−sh + a) ∞ ∂ Y (s − Sk ) ∂s k=−∞
∂ (s + ad e−sh + a) ∂s · · · (S0 − S−2 )(S0 − S−1 )(S0 − S1 ) · · · 1 − ad he−S0 h
s→S0
(27)
If we substitute the above result into (25), we obtain CkI =
x0 − ad G(Sk ) 1 − ad he−Sk h
(28)
With the above results, let us next consider the scalar nonhomogeneous case in (1). The Laplace transform of (1) is,
and comparing numerators yields CkI nk (s) = J(s){x0 − ad G(s)}
=
(22)
k=−∞
∞ X
∞ Y
and using L’Hopital’s rule,
k=−∞
d(s) ≡
(25)
(23)
sX(s) − x0 + ad e−sh X(s) + ad G(s) + aX(s) = (s + ad e−sh + a)X(s) − x0 + ad G(s) = bU (s)
(29)
k=−∞
where J(s) is an undetermined polynomial in s, introduced to equate the numerators and the denominators, which are of infinite order, on both sides of (19) and (17). Even though J(s) is undetermined here, it eventually cancels out later.
Therefore, the transform of the response, x(t) is X(s) =
x0 − ad G(s) bU (s) + s + ad e−sh + a s + ad e−sh + a | {z } | {z } free
forced
(30)
Solving for the unknown X(s) yields 1.2
X(s) = (sI + A + Ad e−sh )−1 {x0 − Ad G(s) + BU(s)} (35) Therefore, x(t) = L−1 (sI + A + Ad e−sh )−1 {x0 − Ad G(s)} | {z } free + L−1 (sI + A + Ad e−sh )−1 {BU(s)} | {z }
1
0.8
Responses, x(t)
0.6
0.4
0.2
0
with 20 branches (k=-9,..., 10) with 4 branches (k=-2,..., 1)
-0.2
forced
with 2 branches (k=-1, 0)
-0.4
-0.6
-0.8 -1
0
1
2
3
4
5
Time, t
Fig. 1. Solution obtained using the Laplace transform combined with the matrix Lambert W function method of 2, 4, 20 branches (straight). Compared to those obtained using the numerical method (dashed), ‘dde 23’ in Matlab, it shows good agreement.
Using the transform of the convolution property of the Laplace transformation, the inverse of the forced term in (30) is bU (s) 1 = bU (s) s + ad e−sh + a s + ad e−sh + a Z t X ∞ ⇐⇒ eSk (t−ξ) CkN bu(ξ)dξ
(31)
0 k=−∞
where, CkN =
1 1 − ad he−Sk h
Note that CkI is dependent on the initial conditions, x0 and the preshape function, g(t), but CkN is not, and (33)
k=−∞
Example 1: Consider the case where a = ad = b = h = 1 in (1) and g(t) = 1, x0 = 1, then, Sk is computed as defined in (2), CkI is computed from (28), and CkN is computed from (32). Table II shows the values, with the branches k = −1, 0, 1 of the Lambert W function, for the total solution in (4) to the system in (1) with the coefficients. Applying these values to (4), we can obtain the solution to (1). The results obtained using 2, 4, and 20 branches are shown in Figure 1 and compared to that obtained using a numerical integration method (‘dde23’ in Matlab). As seen in the figure, with more branches, the results show better agreement. When 20 branches are used, the agreement is excellent. IV. G ENERALIZATION TO S YSTEMS
OF
= · · · + (sI − S−1 )−1 CI−1 + (sI − S0 )−1 CI0 + · · ·
Comparing (37) with the free solution part in (36) yields (sI + A + Ad e−sh )−1 {x0 − Ad G(s)} ∞ X = (sI − Sk )−1 CIk
(38)
k=−∞
= · · · + (sI − S−1 )−1 CI−1 + (sI − S0 )−1 CI0 + · · ·
(38) provides the condition for calculating CIk . Here we provide a 2 × 2 example. If the coefficients are a1 a2 ad1 ad2 A= , Ad = (39) a3 a4 ad3 ad4 we can write the term in (38) as
(32)
∞ X 1 = eSk t CkN s + ad e−sh + a
(36) On the other hand, we have the free solution to (5) as (7). The solution in (7) can be transformed as ∞ X X(s) = (sI − Sk )−1 CIk (37) k=−∞
DDE S
For the system of DDEs in (5), if we take the Laplace transform of both sides we find, sX(s) − x0 + Ad e−sh X(s) + Ad G(s) + AX(s) = BU(s) (34)
(sI + A + Ade−sh ) s + a1 + ad1 e−sh = a3 + ad3 e−sh
a2 + ad2 e−sh s + a4 + ad4 e−sh
(40)
And the inverse of the matrix term in (40) is −1 (sI + A + Ad e−sh ) s + a4 + ad4 e−sh 1 = ∆(s) −(a3 + ad3 e−sh )
−(a2 + ad2 e−sh ) s + a1 + ad1 e−sh (41) where ∆(s) is same as the characteristic equation of (5), that is ∆(s) = s2 + {a1 + a4 + (ad1 + ad4 )e−sh }s +(a1 a4 + a2 a3 ) + (a1 ad4 + ad1 a4 − a2 ad3 + ad2 a3 )e−sh +(ad1 ad4 − ad2 ad3 )e−sh (42) Note that the right side of (38), 1 0 pk1 pk2 sI − Sk = s − 0 1 pk3 pk4 1 0 λk1 0 = s − Vk V−1 , k 0 1 0 λk2 pk1 pk2 where Sk = pk3 pk4 (43) And the inverse is 1 s − pk4 pk2 (sI − Sk )−1 = pk3 s − pk1 (s − λk1 )(s − λk2 ) (44)
TABLE II I NTERMEDIATE RESULTS FOR COMPUTING THE SOLUTION FOR EXAMPLE 1 VIA THE MATRIX L AMBERT W FUNCTION branch
···
k = −1
k=0
k=1
···
Sk
···
−0.6050 − 1.7882i
−0.6050 + 1.7882i
−0.0528 + 7.7184i
···
CkI
···
0.4410 − 0.1541i
0.4410 + 0.1541i
0.0313 − 0.0086i
···
CkN
···
0.2712 − 0.3477i
0.2712 + 0.3477i
−0.0009 − 0.1296i
···
Applying (41) and (44) in (38), we can find the coefficients CIk in (14). For example, to obtain the coefficient of the principal branch, CI0 from (38), 1 s + a4 + ad4 e−sh −(a2 + ad2 e−sh ) ∆(s) −(a3 + ad3 e−sh ) s + a1 + ad1 e−sh × {x(0) − Ad G(s)} 1 s − p04 p02 = p03 s − p01 (s − λ01 )(s − λ02 ) +(sI − S−1 )−1 CI−1 + (sI − S1 )−1 CI1 + · · ·
(46)
∞ X
eSk t CN k
(51)
k=−∞
−1 N C0 + · · · = · · · + (sI − S−1 )−1 CN −1 + (sI − S0 )
∂ (s − λk1 )(s − λk2 ) lim ∂s × s→λ ∆(s) ki −sh λki + a4 + ad4 e −(a2 + ad2 e−sh ) −(a3 + ad3 e−sh ) λki + a1 + ad1 e−sh =
λki − pk4 pk3
pk2 λki − pk1
CN k ,
(53)
fori = 1, 2
(48)
Solving the two equations in (53) simultaneously, we can compute CN k . The above approach can readily be generalized to the case of higher order systems of DDEs.
(49)
Example 2: Consider the case, from [13], where h = 1 and the two initial conditions, g(t) = [1 0]T and x0 = [1 0]T , where T indicates transpose, 1 3 −1.66 0.697 A= , Ad = (54) −2 5 −0.93 0.330
we can write that ∂ (s − λ01 )(s − λ02 ) lim ∂s × s→λ ∆(s) 0i −sh λ0i + a4 + ad4 e −(a2 + ad2 e−sh ) −(a3 + ad3 e−sh ) λ0i + a1 + ad1 e−sh × {x0 − Ad G(λ01 )} λ0i − p04 p02 = CI0 , fori = 1, 2 p03 λ0i − p01
(sI + A + Ad e−sh )−1 ⇐⇒
(47)
This problem can be resolved by application of L’Hopital’s rule as in the scalar case. And assuming where m 6= n; i = 1, 2
where
Here we again consider the 2 × 2 case as an example. Using the same notation as (49) for the general branch, without the terms related to initial condition and preshape function, we should make s approach λki , then,
Then, substitution of λ01 for s in (46) makes the other terms on the right hand side zero except the first term. However, as in the scalar case, because ∆(s) is the characteristic equation, we again encounter a problem:
λmi 6= λni ,
0 k=−∞
The coefficients CN k are computed from the following equation ∞ X (sI + A + Ad e−sh )−1 = (sI − Sk )−1 CN k (52) k=−∞
+(s − λ01 )(s − λ02 )(sI − S−1 )−1 CI−1 +(s − λ01 )(s − λ02 )(sI − S1 )−1 CI1 + · · ·
(s − λ01 )(s − λ02 ) 0 = ∆(s) 0
k=−∞
(45)
If we multiply the both sides by (s − λ01 )(s − λ02 ) (s − λ01 )(s − λ02 ) × ∆(s) −sh s + a4 + ad4 e −(a2 + ad2 e−sh ) −(a3 + ad3 e−sh ) s + a1 + ad1 e−sh × {x(0) − Ad G(s)} s − p04 p02 = CI0 p03 s − p01
where CI0 is computed by solving the two equations in (49) simultaneously. With the result in the previous section, continued from (36) and using the convolution property of the Laplace transform, we can write Z t X ∞ ∞ X Sk t I eSk (t−ξ) CN x(t) = e Ck + k Bu(ξ)dξ (50)
with an external forcing term of cos(t) Bu(t) = sin(t)
(55)
TABLE III I NTERMEDIATE RESULTS FOR COMPUTING THE SOLUTION FOR EXAMPLE 2 IN (54) VIA THE MATRIX L AMBERT W FUNCTION k = −1
branch
Sk
−0.3499 − 4.9801i 2.4174 + 0.1309i
−1.6252 + 0.1459i −5.1048 − 4.5592i
CIk
CN k
k=0
0.0173 − 0.0010i −0.0584 − 0.0585i
0.0767 + 0.1876i 0.0239 + 0.0013i
0.0142 − 0.0419i 0.0196 + 0.1925i
0.3055 2.1317
−1.4150 −3.3015
1
x1
Responses, x(t)
0.6
0.4
0.2
x2
0
-0.2
-0.4
-0.6
-0.8 -1
0
1
2
3
4
0.9847 0.2162
0.3424 −0.0563
1.2
0.8
k=1
5
Time,t
−0.0789 0.4855
−0.3499 + 4.9801i 2.4174 − 0.1309i
−1.6252 − 0.1459i −5.1048 + 4.5592i
0.0173 + 0.0010i −0.0584 + 0.0585i
0.0767 − 0.1876i 0.0239 − 0.0013i
0.0142 + 0.0419i 0.0196 − 0.1925i
of the state transition matrix, can potentially be extended to systems of DDEs. For example, the approach presented based on the matrix Lambert W function, may be useful in controller design via eigenvalue assignment for systems of DDEs. Similarly, concepts of observability, controllability, state estimator design may be tractable and are being studied by the authors. The analytical approach using the matrix Lambert function for ‘time-varying’ DDEs based on Floquet theory is also being currently investigated. These, and others, are all potential topics for future research, which can build upon the foundation presented in this paper. R EFERENCES
Fig. 2. Solution obtained using the Laplace transform combined with the matrix Lambert W function of 11 branches method (straight). Compared to those obtained using the numerical method (dashed), ‘dde 23’ in Matlab, they show good agreement.
Then, Sk is computed from (8) and (9), and CIk is computed by applying (41) and (44) in (38); CN k is computed from (53). Table III shows the values, with the branches k = −1, 0, 1 of the Lambert W function, for the total solution in (14) to the system in (5) with the coefficients in (54). Applying these values into (14), we can obtain the solution to (5). The result obtained using 11 branches is shown in Figure 2 and compared to that obtained using the numerical integration method (‘dde23’ in Matlab). As seen in the figure, the agreement is excellent. V. C ONCLUSIONS
AND
F UTURE W ORKS
In this paper, the Laplace transform approach to solve linear systems of DDEs is developed based upon the matrix Lambert W function method [8], [9]. To provide a closed form solution to systems of linear DDEs, in a form similar to systems of ordinary differential equations, is the essential advantage of the matrix Lambert W function method. With the Laplace transform approach presented here, the analogy can be reinforced. Table I summarizes the comparison between the ODE and DDE cases. The concept of the state transition matrix in ODEs can be generalized to DDEs using the matrix Lambert W function. This suggests that some analyses used in systems of ODEs, based on the concept
[1] J. P. Richard, “Time-delay systems: an overview of some recent advances and open problems,” Automatica, vol. 39, 2003, pp. 16671694. [2] H. Gorecki, S. Fuksa, P. Grabowski, and A. Korytowski, Analysis and Synthesis of Time Delay Systems, John Wiley and Sons, PWN-Polish Scientific Publishers Warszawa, 1989. [3] J. Lam, “Model Reduction of Delay Systems Using Pade Approximants,” Int. J. Control, Vol. 57, No. 2, pp. 377-391, 1993 [4] C. E. Falbo, “Analytic and Numerical Solutions to the Delay Differential Equations,” Joint Meeting of the Northern and Southern California Sections of the MAA, San Luis Obispo, CA, 1995 [5] E. M. Wright, “The Non-Linear Difference-Differential Equation,” Q. J. Math., Vol. 17, pp. 245-252, 1946 [6] R. E. Bellman and K. L. Cooke, Differential-Difference Equations, New York: Academic Press, 1963. [7] E. M. Wright, “The Linear Difference-Differential Equation With Asymptotically Constant Coefficients,” Am. J. Math., Vol. 70, No. 2, pp. 221-238, 1948 [8] F. M. Asl, and A. G. Ulsoy, “Analysis of a System of Linear Delay Differential Equations,” ASME J. Dyn. Syst. Meas. Control, Vol. 125, No. 2, pp. 215-223, 2003 [9] S. Yi and A. G. Ulsoy, “Solution of a system of linear delay differential equations using the matrix Lambert function,” Proc. 25th American Control Conference, Minneapolis, MN, Jun. 2006, pp. 2433-2438. [10] M. C. Pease, Methods of Matrix Algebra. New York: Academic Press, 1965, pp. 406. [11] H. Shinozaki, Robust Stability Analysis of Linear Time-Delay Systems by Lambert W Function. Master thesis, Dept. Electro. Inform. Sci., Kyoto Institute of Technology, Kyoto, Japan, 2003 [12] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J.Jeffrey, and D. E. Knuth, “On Lambert’s W function,” Adv. Comput. Math., Vol. 5, No.4, 1996, pp. 329-359. [13] T. N. Lee and S. Dianat, “Stability of Time-Delay Systems,” IEEE Trans. Aut. Control, Vol. AC-26, No. 4, pp. 951-953, 1981 [14] M. Malek-zavarei, M. Jamshidi, Time-delay systems, North Holland, New York, 1987
