- Email: [email protected]
Computation of critical parameters for neutral type time delay systems
Computation of critical parameters for neutral type time delay systems G. Ochoa ∗ S. Mondi´ e ∗ V. L. Kharitonov ∗∗ ∗
Automatic Control Department, CINVESTAV-IPN, Av. IPN, C.P. 07360 D.F., M´exico, (e-mail: gochoa;[email protected]). ∗∗ Applied Mathematics and Control Process Department, St.-Petersburg State University, St.-Petersburg, Russia, (e-mail: [email protected])
Abstract: In this paper, a procedure for the computation of the Lyapunov matrix of neutral type systems, with delays multiple of a basic one, is recalled: It consists in solving a boundary value problem for a delay free system of matrix equations. The important property that the elements of the spectrum of the delay system that are symmetric with respect to the imaginary axis belong to the spectrum of the delay free system is also established. This property is exploited for the determination of the critical values of the time delay system. Keywords: Critical parameters, Lyapunov matrix, Neutral type systems 1. INTRODUCTION A number of contributions on the theoretical analysis of neutral type linear time delays systems have been contributed in the past few years: The form of the Lyapunov functionals with prescribed derivative parameterized in the Lyapunov matrix was presented in Rodriguez et al. (2004) for the single delay case. This matrix has been shown to be the solution of a delay free system subject to boundary conditions Louisell (2001), Kharitonov (2005) and its existence and uniqueness were established in Kharitonov (2009). The extension of some of these results to the case of delays multiple of a basic one was recently addressed in Ochoa et al. (2009). In this paper, we present for the case of neutral type systems commensurate delays, a methodology that reveals the values of the parameters (called critical parameters) for which the system may cross stability/instability boundaries. This methodology is based on the same concepts as those presented for distributed and point wise delay systems Ochoa et al. (2009). It exploits the interesting connection between the system spectrum and that of the delay free system used in the computation of the Lyapunov matrix: If s0 is a root of the delay system such that −s0 is also a root, then s0 belongs also to the spectrum of the delay free system. This methodology differs from other approaches for the computation of critical values, see Ivanescu et al. (2003), Jarlebring (2007), Louisell (1998), Olgac and Sipahi (2004), Walton and Marshall (1987) and the references therein, by the fact that it is strongly connected to time domain properties. The structure of the contribution is as follows: In Section 2 some basic results and the concept of Lyapunov matrix are introduced. Section 3 is devoted to recall the analytic This work was supported by CONACyT, project No. 61076
procedure for the computation of the Lyapunov matrix that is based in the fact that the Lyapunov matrix satisfies a delay free system of matrix equations. The relation between the spectrum of the delay free system and those of the neutral type time delay system is established in Section 4. Based on this relation a methodology for the computation of critical parameters is proposed in Section 5 and some examples are analyzed in Section 6. The contribution ends with some concluding remarks. 2. PRELIMINARY RESULTS We introduce next some definitions and useful results. 2.1 System description In this paper we consider the class of neutral type systems studied in Ochoa et al. (2009): m m Bk x˙ (t − hk) = Ak x (t − hk) , t ≥ 0. (1) k=0
k=0
Here Ak , Bk ∈ Rn×n are given matrices, h > 0 is the basic delay. Without any loss of generality we assume that B0 = I.
Given an initial condition ϕ ∈ C([−mh, 0], R), there exists a unique solution x(t, ϕ) of system (1). We denote by xt (ϕ) the segment of the solution xt (ϕ) → x(t + θ, ϕ), θ ∈ [−mh, 0]. Definition 1. Bellman and Cooke (1963) System (1) is exponentially stable if there exist constants σ > 0 and γ ≥ 1 such that for every solution x(t, ϕ) of the system the following estimates holds x(t, ϕ) ≤ γe−σt ϕmh , for t ≥ 0. Here, ϕmh is defined as ϕmh = max ϕ (θ) . θ∈[−mh,0]
• The algebraic property
2.2 Lyapunov-Krasovskii functionals
m m Bj U (h (j − k)) Ak +
The form of the functional with prescribed time derivative associated to system (1) is obtained in Ochoa et al. (2009), Rodriguez et al. (2004). This results reads as follows: Let us assume that system (1) is exponentially stable. Given a functional w (xt ) of the form w (xt ) =
m
x (t − hk) Wk x (t − hk) +
k=0
+
m 0
x (t + θ) Wm+k x (t + θ) dθ,
k=1−hk
where Wk , k = 0, 1, ..., 2m are positive definite matrices, there exists a quadratic functional v (xt ) , , such that d v (xt ) = −w (xt ) dt along the solutions of system (1). The functional v (xt ) is of the form m m v (xt ) = x (t) Bj U (h (j − k)) Bk x (t) + +2
x
(t) Bj
j=0 k=1
0
+ Ak U (h (j − k)) Bj ] = −W.
3. COMPUTATION OF THE LYAPUNOV MATRIX In Ochoa et al. (2009) a methodology for the computation of the Lyapunov matrix for system (1) was introduced. This procedure consists in showing that the Lyapunov matrix satisfies a delay free system of matrix equations where the initial conditions are not known and the problem is reduced to find the initial conditions via a two point boundary value problem. Due to its importance in the present contribution, we recall this semi-analytic procedure below. First, the following set of 2m auxiliary matrices is introduced Xj (τ ) = U (τ + jh) , j = −m, ..., 0, ..., m − 1. (5) Let us consider j ≥ 0 and compute the derivative of matrices Xj (τ ) with respect to τ , that is
j=0 k=0
m m
d d Xj (τ ) = U (τ + jh) , j = 0, ..., m − 1. dτ dτ Since the matrices defined in (5) satisfy the properties (2)-(4), by using the dynamic property (2), the above expression can be rewritten as
U (h (j − k) − θ) ×
−hk
× [Ak x (t + θ) − Bk x˙ (t + θ)] dθ+ m m 0 0 [x (t + θ1 ) Ak − +
Xj (τ ) = −
j=1 k=1−hj −hk
θ2 ) Bk ] U
−x˙ (t + (h (j − k) + θ1 − θ2 ) × × [Ak x (t + θ) − Bk x˙ (t + θ)] dθ2 dθ1 + m 0 x˙ (t + θ) [Wk + (hk + θ) Wm+k ] × +
+ equivalently, m
× x (t + θ) dθ. U (τ ) =
∞
m
k=1 m
U (τ + h (j − k)) Bk + U (τ + h (j − k)) Ak ,
k=0
k=1−hk
Here
Xj−k
(τ ) Bk =
k=0
K (t) W K (t + τ ) dτ,
0
is called the Lyapunov matrix of system (1) associated to m the matrix W = W0 + [Wk + hkWm+k ] and matrix
m
then
m
K (t) is the fundamental matrix of system (1).
Bk Xj+k (τ ) = −
k=0
2.3 Properties of the Lyapunov matrix The Lyapunov matrix satisfies a set of properties that are the basis for a methodology for its computation, Kharitonov (2005). These properties are: • The dynamic property m [U (τ − hk) Bk − U (τ − hk) Ak ] = 0, τ > 0, k=0
(2) (3)
Xj−k (τ ) Ak .
(6)
k=0
Now for j = −m, ..., −1 the derivative of the auxiliary matrices is Xj (τ ) = − [U (−τ − jh)]
k=1
• The symmetry property U (τ ) = U (−τ ) , τ ≥ 0,
(4)
j=0 k=0
m
Ak Xj+k (τ ) .
(7)
k=0
Collecting the matrix equations (6) and (7) we arrive to the following set of delay free matrix equations: m m Xj−k (τ ) Bk = Xj−k (τ ) Ak , j ≥ 0, k=0 m
k=0
Bk Xj+k
k=0 m
(τ ) = −
(8)
Ak Xj+k
(τ ) , j < 0.
k=0
The vectorization of the above system can be written as M1 z (τ ) = M2 z (τ ) , or z (τ ) = M z (τ )
with M = M1−1 M2 , for details about the existence of M −1 see Gohberg and Lerer (1976). Here M1 , M2 ∈ 2 2 R2mn ×2mn ⎡are defined as follows ⎤ ¯m 0 · · · ¯1 · · · B ¯0 B 0 0 B ¯m−1 B ¯m · · · ¯0 · · · B ⎢0 B 0 0 ⎥ ⎢ ⎥ .. .. . . .. .. ⎥ ⎢ .. .. . . ⎢ . . ⎥ . . . . . . ⎢ ¯m−1 B ¯m ⎥ ⎢ 0 0 ··· ⎥ 0 0 · · · B ⎥ M1 = ⎢ ˜m 0 · · · ˜0 B ˜1 · · · B ⎢B ⎥ 0 0 ⎢ ⎥ ˜0 · · · B ˜m−1 B ˜m · · · ⎢0 B ⎥ 0 0 ⎢ ⎥ ⎢ .. .. . . .. .. . . .. .. ⎥ ⎣ . . . . . . . . ⎦ ˜ ˜ 0 0 ··· 0 0 · · · Bm−1 Bm and ⎡ ⎤ 0 ··· 0 0 A¯0 A¯1 · · · A¯m ⎢ 0 0 0 ⎥ A¯0 · · · A¯m−1 A¯m · · · ⎢ ⎥ .. . . .. .. .. .. ⎥ ⎢ .. .. ⎢ . . . . . . . . ⎥ ⎢ ⎥ ¯ ¯ ⎢ 0 0 ··· 0 0 · · · Am−1 Am ⎥ ⎥. M2 = ⎢ ⎢−A˜0 −A˜1 · · · −A˜m 0 ··· 0 0 ⎥ ⎢ ⎥ ⎢ 0 −A˜0 · · · −A˜m−1 −A˜m · · · 0 0 ⎥ ⎢ ⎥ ⎢ .. .. . . .. .. .. .. ⎥ .. ⎣ . . . . . . . . ⎦ 0 0 ··· 0 0 · · · −A˜m−1 −A˜m
¯j = B ⊗ I where A¯j = Aj ⊗ I , A˜j = I ⊗ Aj , B j
˜j = I ⊗ B , j = 0, 1, ..., m. The initial conditions and B j of the above system are not known, but we can find them by solving the following boundary value problem: a first set of boundary conditions follows from the definition of the auxiliary matrices introduced in (5): Xj+1 (0) = Xj (h) , j = −m, ..., 0, ..., m − 2, and a second could be found by expressing the algebraic property (4) in terms of the auxiliary matrices (5): m −W = [Bk X0 (0) Ak + Ak X0 (0) Bk ] k=0
+
m m
Bk Xk−j−1 (h) Aj + Bj Xk−j−1 (h) Ak
j=0 k=j+1
+
m m j=0 k=j+1
Aj Xk−j−1 (h) Bk + Ak Xk−j−1 (h) Bj .
4. FREQUENCY DOMAIN STABILITY ANALYSIS In this section the key relationship between the spectrum of the delay free system (8) and the roots of the characteristic function of system (1) is established, Kharitonov (2006), Louisell (2001). Theorem 1. Let s0 be an eigenvalue of the neutral delay system (1), such that −s0 is also an eigenvalue of system (1). Then s0 belongs to the spectrum of the delay free system (8). Proof. The characteristic equation of system (1) is given by m G (s, h) = sBk e−shk − Ak e−shk .
γ G(s0 , h) = 0 and G (−s0 , h) μ = 0. (9) If a complex number s belongs to the spectrum of the delay free system (8) then there exist a set of 2m non trivial matrices Xj , j = −m, ..., 0, ...m − 1 such that the following holds: m [sXj−k Bk − Xj−k Ak ] =0, j ≥ 0, (10) k=0
m
[sBk Xj+k + Ak Xj+k ] =0,
j < 0.
k=0
Computing e−s0 hj μ [γ G(s0 , h)] = 0, j = 0, ..., m one arrives to the expression m s0 Bk es0 h(j−k) − Ak es0 h(j−k) = 0 μγ
(11)
k=0
for j = 0, ..., m.
Now if we calculate [G (−s0 , h)μ] γ es0 hj = 0, j = −m, ..., −1 it follows that m s0 Bk es0 h(j+k) + Ak es0 h(j+k) μγ = 0.
(12)
k=0
If we define the matrices ¯ l = es0 hl μγ , l = −m, ..., 0, ..., m − 1, X the matrix equations (11)-(12) take the form m m ¯ j−k Bk = ¯ j−k Ak , j ≥ 0, s0 X X k=0
m
¯ j+k = − s0 Bk X
k=0
k=0 m
¯ j+k , j < 0. Ak X
k=0
¯ l are non trivial, then It is clear that the matrices X the value s0 belongs to the delay free system (8). This procedure also holds for −s0 . Proposition 2. The spectrum of the delay free system (8) is symmetrical with respect to the imaginary axis. Proof. The fact that the spectrum of the delay free system of matrix equations (8) is symmetrical with respect of the imaginary axis follows from the idea that if for a value s there exist non trivial matrices Xi i = −m, ..., 0, ...m − 1 that satisfy system (10). By transposing and multiplying by (−1) system (10) it follows that m −sBk Xj−k =0, + Ak Xj−k k=0 m k=0
−sXj+k Bk − Xj+k Ak =0,
and by setting the matrices ˜j = X X −j−1 , j = −m, ..., 0, ..., m − 1 we conclude that the system is satisfied (10) for −s. 5. PROPOSED METHODOLOGY
k=0
Since the values s0 and −s0 belong to the spectrum of the neutral delay system (1), there exist non zero vectors μ and γ, such that
In this section two different methodologies are proposed: The first one applies to systems with an unknown parameter in the matrices of the system. The second methodology
6. EXAMPLES
applies to systems where the only unknown parameter is the delay h. 5.1 Methodology A Here we consider the class of systems that are of the form m m A˜k x (t − hk) . Bk x˙ (t − hk) = k=0
k=0
where A˜k = kμAk , μ is an unknown parameter that satisfies μ > 0. The methodology for the computation of the critical parameters is as follows: (1) Calculate the corresponding delay free system (2) Compute the characteristic polynomial p(s, ρ) of the delay free system (8). The value ρ represent the delay h or the unknown parameter μ (or both) (3) As the spectrum of the system is symmetrical with respect to the imaginary axis, all the powers of the polynomial are even, then we can set λ = s2 and analyze the polynomial p1 (λ, ρ) (4) Determining the purely imaginary roots of p(s, ρ) reduces to find the negative real roots of p1 (λ, ρ) (5) Using Sturm’s Theorem Uspensky (1948), determine the negative real roots - Compute the Sturm sequence of the polynomial p1 (λ, ρ) - Determine with the help of Sturm’s Theorem if p1 (λ, ρ) has negative real roots (6) If there exist a negative real roots compute the corresponding purely imaginary roots and form a set of candidate roots for the delay system (7) For each of these candidate roots, check if there exists a value of the parameter ρ that annihilate the characteristic function. If this is the case then the pair (ωi∗ , ρ∗i ) , receives the name of critical frequency and critical parameter (8) Form a table of ρ∗i in ascending order (consider ρ0 = 0) (9) For ρ∗i , ρ∗i+1 and ωi∗ determine the number of unstable roots of the characteristic function using the argu- ment principle to determine if the interval ρ∗i , ρ∗i+1 is a stable or unstable region. 5.2 Methodology B For this methodology we consider the class of neutral type time delay systems of the form (3), the methodology consists of the following steps: a. The steps (1)-(3) of Methodology A remains valid b. Compute the negative real roots and the corresponding purely imaginary roots c. Make the change of variable z = e−sh d. Evaluate the characteristic function f (s, h) for s = iω obtained from step (b) e. Compute the roots of the polynomial f (z) and select those who satisfies |z| = 1 f. Compute the critical delays associated to the critical frequency ω ∗ , that is −Arg (z) h∗ = ω∗ g. Continue with Steps (8)-(9) of Methodology A for the parameter h.
In this section we illustrate the use of the proposed methodologies for the computation of critical parameters. Example 3. Let us consider the following system studied in Ivanescu et al. (2003), Rodriguez et al. (2004) x˙ (t) + B1 x˙ (t − h) = A0 x (t) + A1 x (t − h) where −2 0 −1 0 , A1 = , A0 = 0 −0.9 −1 −1 −0.1 0 B1 = 0 −0.1
Since ρ (B1 ) = 0.1 < 1 the difference operator is stable. The corresponding delay free system is of the form (I ⊗ I) (B1 ⊗ I) z (τ ) = (I ⊗ B1 ) (I ⊗ I) (A0 ⊗ I) (A1 ⊗ I) z (τ ) . −(I ⊗ A1 ) −(I ⊗ A0 )
In order to find the critical delays, we first compute the characteristic polynomial of the delay free system p (s) = s8 − 5. 689 1s6 + 8. 1628s4 − 0.1955s2 − 0.3797.
By setting λ = s2 , and analyzing the auxiliary polynomial p1 (λ) we can directly check that the only negative real root 19 is λ = − 99 which corresponds to the purely imaginary root √ √ 11 19 s = i 33 . Now, the characteristic function of the delay system is f (s, h) = det(sI + sB1 e−sh − A0 − A1 e−sh ). Evaluating f (s, h) for s = iω and the change of variable z = e−sh ; f (s, h) take the following form f (z) = z 2 (0.998 0 − 0.0087i) + z (2. 938 3 + 0.749 1i) + (1. 608 0 + 1. 270 4i) with roots
z1 = −0.879 1 − 0.4765i, z2 = −1. 977 0 − 0.524 6 i.
Now we are looking for roots z1,2 that satisfies |z| = 1, the only root that satisfies this condition is z1 . The critical delay is computed as −Arg (z1 ) h∗ = = 6.0371. ω∗ Finally with the help of the argument principle we can conclude that the system is stable for h ∈ [0, 6.0371) and is unstable for h ∈ (6.0371, ∞). Example 4. In this example we consider a system with an unknown parameter μ that satisfies μ > 0. The system is of the form x˙ (t) + B1 x˙ (t − 1) = μA0 x (t) + A1 x (t − 1) where A0 , A1 and B1 are the same as in the previous example. The characteristic polynomial of the corresponding delay free system is 1 p(s, μ) = 4 q (s, μ) r (s, μ) v (s, μ) 99 where
Re f(ω1 (μ ), μ ) Im f(ω1 (μ ), μ )
3 2 1 0 −1 −2 −3
0
0.1
0.2
0.3
μ
0.4
0.5
Fig. 1. Real and Imaginary parts of f (iω1 , μ) 8 Re f(ω2 (μ ), μ ) Im f(ω2 (μ ), μ )
6
−5.5779 5.0852 3.5478 A0 = −1.7151 2.2988 0.8782 , 0.4407 3.2284 −1.3174 0 −0.2 0.4 0.3 0.1 0 0 −0.2 0 B1 = 0.5 −0.3 0 , B2 = , −0.2 −0.7 0 −0.1 0 −0.4 The matrices of the delay free system are ⎤ ⎡ 0 I ⊗ I B1 ⊗ I B2 ⊗ I I ⊗ I B1 ⊗ I B2 ⊗ I ⎥ ⎢ 0 M1 = ⎣ I ⊗ B2 I ⊗ B1 I ⊗ I 0 ⎦ 0 I ⊗ B2 I ⊗ B1 I ⊗ I and ⎡ ⎤ 0 0 0 A0 ⊗ I A0 ⊗ I 0 0 ⎢ 0 ⎥ M2 = ⎣ ⎦ 0 0 −I ⊗ A0 0 0 0 0 −I ⊗ A0 The expression of the characteristic polynomial of this delay free system, which is of degree 32, is omitted due to space limitations. With the change of variable λ = s2 the polynomial p1 (λ) is obtained. We can verify that this polynomial has a pair of negative real roots
4
4
λ1 = −0.8807 λ2 = −1.3665 × 10−4
2
with corresponding purely imaginary roots
0
s1 = 0.9384i
−2
s2 = 0.0116i. −4
0
0.2
0.4
μ
0.6
0.8
1
Fig. 2. Real and Imaginary parts of f (iω2 , μ) q (s, μ) = 99s − 81μ + 100 2
2
r (s, μ) = 99s2 − 400μ2 + 100
v (s, μ) = 9801s4 + s2 −47 740μ2 + 19 800
+ 32 400μ4 − 36 000μ2 + 10000
With the change of variable λ = s2 we can directly check that the auxiliary polynomial q1 (λ, μ) has one negative
, while the auxiliary real root for μ ∈ 0, 10 9
polynomial r1 (λ, μ) has one negative root for μ ∈ 0, 12 . With the help of the Sturm array we can conclude that the auxiliary polynomial v1 (λ, μ) does not have negative real roots ∀μ > 0. The characteristic function of the system is of the form
f (s, μ) = det sI + sB1 e−s − μA0 − A1 e−s .
We must check if there exist pairs {iωk (μ) , μ} k = 1, 2 that annihilate the characteristic function of the system. By direct calculations we can prove that for μ1 = 0.4615 and ω1∗ = 0.3867, see Figure 1, and for ω2∗ = 0.3866 and μ2 = 1.0256, see Figure 2, the real and imaginary parts vanish simultaneously Example 5. Now let us consider the system introduced in Michiels and Vyhl´ıdal (2005), Jarlebring (2007), x˙ (t) + B1 x˙ (t − h) + B2 x˙ (t − 2h) = A0 x (t) , where
Now the characteristic function of the system is of the form
f (s, h) = det sI + sB1 e−sh + sB2 e−2sh − A0 .
With the change of variable z = e−sh and considering s1 = 0.9384i, f (s, h) takes the form f (z) = − (0.0198i) z 6 − (0.0396i) z 5 − (0.6969 − 0.1388i)z 4 + (0.6452 + 0.119i) z 3 + (2.9265 − 0.9274 i) z 2 + (0.2915 +6.5017i) z − (3.994 6 + 4.748 9i) The only root zi that satisfies |z| = 1 is z1 = 0.9839 − 0.1785i, then the critical delay is Arg(z1 ) h∗1 = − = 0.1912. ω1∗ If we compute f (s, z) for ω2 = 0.0116 we obtain the polynomial
f (z) = − 7.667 × 10−8 i z 6 − 1.533 × 10−7 i z 5
− 2.162 × 10−4 − 5.367 × 10−7 i z 4
+ 0.0536 + 2.3 × 10−7 i z 3
+ 4.5408 × 10−4 − 0.12788i z 2
+ 4.523 × 10−5 + 0.0779 × 10−2 i z + (0.0523 − 0.09i)
that has roots z1 = −9.6132 + 7.6772i z3 = −39.8206 − 38.2932i z2 = 0.5582 + 0.8296i z4 = 36.811 7 + 38.296 0i z3 = 5.301 2 + 0.07.574i z5 = 4.762 5 − 8.585 4i
The only root that satisfies |z| = 1 is z2 then Arg(z2 ) 2π + ∗ = 453.7893. h∗2 = − ω2∗ ω2 The argument principle shows that the only interval where the system is stable is for h ∈ [0, 0.1912). 7. CONCLUSIONS A procedure for the computation of the critical parameters of time delay systems of neutral type is presented. The approach takes advantage of the key relation between the roots of the characteristic polynomial of a delay free system of matrix equations and the roots of the characteristic function of the neutral type time delay system. Some illustrative examples show the effectiveness of the method for determining the critical delays as well as other critical parameters. REFERENCES Bellman, R. and Cooke, K.L. (1963). DifferentialDifference Equations. Academic Press, New York. Gohberg, I.C. and Lerer, L.E. (1976). Resultants of matrix polynomials. Bulletin of the American Mathematical Society, 82, 565–567. Ivanescu, D., Niculescu, S.I., Dugard, L., Dion, J.M., and Verriest, E.I. (2003). On delay-dependent stability for linear neutral systems,. Automatica, 39, 255–261. Jarlebring, E. (2007). On critical delays for linear neutral delay systems. European Control Conference, Kos, Greece. Kharitonov, V.L. (2005). Lyapunov functionals and lyapunov matrices for neutral type time delay systems: a single delay case. International Journal of Control, 78, 783–800. Kharitonov, V.L. (2006). Lyapunov matrices for a class of time delay systems. Systems and Control Letters, 55, 610–617. Kharitonov, V.L. (2009). Lyapunov matrices: Existence and uniqueness issues. 8th Workshop on Time Delay Systems, Sinaia, Romania. Louisell, J. (1998). Numerics of the stability exponent and eigenvalue abscissas of a matrix delay system. In L. Dugard and E.I. Verriest (eds.), In Lecture Notes in Control and information sciences 228, volume 228. Springer. Louisell, J. (2001). A matrix method for determining the imaginary axis eigenvalues of a delay system. IEEE Transactions on Automatic Control, 46, 2008–2012. Michiels, W. and Vyhl´ıdal, T. (2005). An eigenvalue based approach for the stabilization of linear time-delay systems of neutral type. Automatica, 41, 991–998. Ochoa, G., Vel´ azquez, J.E.V., Kharitonov, V.L., and Mondi´e, S. (2009). Lyapunov matrices for neutral type time delay systems. In J. Loiseau, W. Michiels, S.I. Niculescu, and R. Sipahi (eds.), In Lecture Notes in Control and information sciences 388, volume 388. Springer. Olgac, N. and Sipahi, R. (2004). A practical method for analyzing the stability of neutral type lti-time delayed systems. Automatica, 40, 847–853. Rodriguez, S.A., Kharitonov, V.L., Dion, J.M., and Dugard, L. (2004). Robust stability of neutral systems:
a lyapunov krasovskii constructive approach. International Journal of Robust and Nonlinear Control, 14, 1345–1358. Uspensky, J.V. (1948). Theory of Equations. McGraw-Hill, New York. Walton, K.E. and Marshall, J.E. (1987). Direct method for tds stability analysis. IEEE proceedings, 134, 101–107.
Automatic Control Department, CINVESTAV-IPN, Av. IPN, C.P. 07360 D.F., M´exico, (e-mail: gochoa;[email protected]). ∗∗ Applied Mathematics and Control Process Department, St.-Petersburg State University, St.-Petersburg, Russia, (e-mail: [email protected])
Abstract: In this paper, a procedure for the computation of the Lyapunov matrix of neutral type systems, with delays multiple of a basic one, is recalled: It consists in solving a boundary value problem for a delay free system of matrix equations. The important property that the elements of the spectrum of the delay system that are symmetric with respect to the imaginary axis belong to the spectrum of the delay free system is also established. This property is exploited for the determination of the critical values of the time delay system. Keywords: Critical parameters, Lyapunov matrix, Neutral type systems 1. INTRODUCTION A number of contributions on the theoretical analysis of neutral type linear time delays systems have been contributed in the past few years: The form of the Lyapunov functionals with prescribed derivative parameterized in the Lyapunov matrix was presented in Rodriguez et al. (2004) for the single delay case. This matrix has been shown to be the solution of a delay free system subject to boundary conditions Louisell (2001), Kharitonov (2005) and its existence and uniqueness were established in Kharitonov (2009). The extension of some of these results to the case of delays multiple of a basic one was recently addressed in Ochoa et al. (2009). In this paper, we present for the case of neutral type systems commensurate delays, a methodology that reveals the values of the parameters (called critical parameters) for which the system may cross stability/instability boundaries. This methodology is based on the same concepts as those presented for distributed and point wise delay systems Ochoa et al. (2009). It exploits the interesting connection between the system spectrum and that of the delay free system used in the computation of the Lyapunov matrix: If s0 is a root of the delay system such that −s0 is also a root, then s0 belongs also to the spectrum of the delay free system. This methodology differs from other approaches for the computation of critical values, see Ivanescu et al. (2003), Jarlebring (2007), Louisell (1998), Olgac and Sipahi (2004), Walton and Marshall (1987) and the references therein, by the fact that it is strongly connected to time domain properties. The structure of the contribution is as follows: In Section 2 some basic results and the concept of Lyapunov matrix are introduced. Section 3 is devoted to recall the analytic This work was supported by CONACyT, project No. 61076
procedure for the computation of the Lyapunov matrix that is based in the fact that the Lyapunov matrix satisfies a delay free system of matrix equations. The relation between the spectrum of the delay free system and those of the neutral type time delay system is established in Section 4. Based on this relation a methodology for the computation of critical parameters is proposed in Section 5 and some examples are analyzed in Section 6. The contribution ends with some concluding remarks. 2. PRELIMINARY RESULTS We introduce next some definitions and useful results. 2.1 System description In this paper we consider the class of neutral type systems studied in Ochoa et al. (2009): m m Bk x˙ (t − hk) = Ak x (t − hk) , t ≥ 0. (1) k=0
k=0
Here Ak , Bk ∈ Rn×n are given matrices, h > 0 is the basic delay. Without any loss of generality we assume that B0 = I.
Given an initial condition ϕ ∈ C([−mh, 0], R), there exists a unique solution x(t, ϕ) of system (1). We denote by xt (ϕ) the segment of the solution xt (ϕ) → x(t + θ, ϕ), θ ∈ [−mh, 0]. Definition 1. Bellman and Cooke (1963) System (1) is exponentially stable if there exist constants σ > 0 and γ ≥ 1 such that for every solution x(t, ϕ) of the system the following estimates holds x(t, ϕ) ≤ γe−σt ϕmh , for t ≥ 0. Here, ϕmh is defined as ϕmh = max ϕ (θ) . θ∈[−mh,0]
• The algebraic property
2.2 Lyapunov-Krasovskii functionals
m m Bj U (h (j − k)) Ak +
The form of the functional with prescribed time derivative associated to system (1) is obtained in Ochoa et al. (2009), Rodriguez et al. (2004). This results reads as follows: Let us assume that system (1) is exponentially stable. Given a functional w (xt ) of the form w (xt ) =
m
x (t − hk) Wk x (t − hk) +
k=0
+
m 0
x (t + θ) Wm+k x (t + θ) dθ,
k=1−hk
where Wk , k = 0, 1, ..., 2m are positive definite matrices, there exists a quadratic functional v (xt ) , , such that d v (xt ) = −w (xt ) dt along the solutions of system (1). The functional v (xt ) is of the form m m v (xt ) = x (t) Bj U (h (j − k)) Bk x (t) + +2
x
(t) Bj
j=0 k=1
0
+ Ak U (h (j − k)) Bj ] = −W.
3. COMPUTATION OF THE LYAPUNOV MATRIX In Ochoa et al. (2009) a methodology for the computation of the Lyapunov matrix for system (1) was introduced. This procedure consists in showing that the Lyapunov matrix satisfies a delay free system of matrix equations where the initial conditions are not known and the problem is reduced to find the initial conditions via a two point boundary value problem. Due to its importance in the present contribution, we recall this semi-analytic procedure below. First, the following set of 2m auxiliary matrices is introduced Xj (τ ) = U (τ + jh) , j = −m, ..., 0, ..., m − 1. (5) Let us consider j ≥ 0 and compute the derivative of matrices Xj (τ ) with respect to τ , that is
j=0 k=0
m m
d d Xj (τ ) = U (τ + jh) , j = 0, ..., m − 1. dτ dτ Since the matrices defined in (5) satisfy the properties (2)-(4), by using the dynamic property (2), the above expression can be rewritten as
U (h (j − k) − θ) ×
−hk
× [Ak x (t + θ) − Bk x˙ (t + θ)] dθ+ m m 0 0 [x (t + θ1 ) Ak − +
Xj (τ ) = −
j=1 k=1−hj −hk
θ2 ) Bk ] U
−x˙ (t + (h (j − k) + θ1 − θ2 ) × × [Ak x (t + θ) − Bk x˙ (t + θ)] dθ2 dθ1 + m 0 x˙ (t + θ) [Wk + (hk + θ) Wm+k ] × +
+ equivalently, m
× x (t + θ) dθ. U (τ ) =
∞
m
k=1 m
U (τ + h (j − k)) Bk + U (τ + h (j − k)) Ak ,
k=0
k=1−hk
Here
Xj−k
(τ ) Bk =
k=0
K (t) W K (t + τ ) dτ,
0
is called the Lyapunov matrix of system (1) associated to m the matrix W = W0 + [Wk + hkWm+k ] and matrix
m
then
m
K (t) is the fundamental matrix of system (1).
Bk Xj+k (τ ) = −
k=0
2.3 Properties of the Lyapunov matrix The Lyapunov matrix satisfies a set of properties that are the basis for a methodology for its computation, Kharitonov (2005). These properties are: • The dynamic property m [U (τ − hk) Bk − U (τ − hk) Ak ] = 0, τ > 0, k=0
(2) (3)
Xj−k (τ ) Ak .
(6)
k=0
Now for j = −m, ..., −1 the derivative of the auxiliary matrices is Xj (τ ) = − [U (−τ − jh)]
k=1
• The symmetry property U (τ ) = U (−τ ) , τ ≥ 0,
(4)
j=0 k=0
m
Ak Xj+k (τ ) .
(7)
k=0
Collecting the matrix equations (6) and (7) we arrive to the following set of delay free matrix equations: m m Xj−k (τ ) Bk = Xj−k (τ ) Ak , j ≥ 0, k=0 m
k=0
Bk Xj+k
k=0 m
(τ ) = −
(8)
Ak Xj+k
(τ ) , j < 0.
k=0
The vectorization of the above system can be written as M1 z (τ ) = M2 z (τ ) , or z (τ ) = M z (τ )
with M = M1−1 M2 , for details about the existence of M −1 see Gohberg and Lerer (1976). Here M1 , M2 ∈ 2 2 R2mn ×2mn ⎡are defined as follows ⎤ ¯m 0 · · · ¯1 · · · B ¯0 B 0 0 B ¯m−1 B ¯m · · · ¯0 · · · B ⎢0 B 0 0 ⎥ ⎢ ⎥ .. .. . . .. .. ⎥ ⎢ .. .. . . ⎢ . . ⎥ . . . . . . ⎢ ¯m−1 B ¯m ⎥ ⎢ 0 0 ··· ⎥ 0 0 · · · B ⎥ M1 = ⎢ ˜m 0 · · · ˜0 B ˜1 · · · B ⎢B ⎥ 0 0 ⎢ ⎥ ˜0 · · · B ˜m−1 B ˜m · · · ⎢0 B ⎥ 0 0 ⎢ ⎥ ⎢ .. .. . . .. .. . . .. .. ⎥ ⎣ . . . . . . . . ⎦ ˜ ˜ 0 0 ··· 0 0 · · · Bm−1 Bm and ⎡ ⎤ 0 ··· 0 0 A¯0 A¯1 · · · A¯m ⎢ 0 0 0 ⎥ A¯0 · · · A¯m−1 A¯m · · · ⎢ ⎥ .. . . .. .. .. .. ⎥ ⎢ .. .. ⎢ . . . . . . . . ⎥ ⎢ ⎥ ¯ ¯ ⎢ 0 0 ··· 0 0 · · · Am−1 Am ⎥ ⎥. M2 = ⎢ ⎢−A˜0 −A˜1 · · · −A˜m 0 ··· 0 0 ⎥ ⎢ ⎥ ⎢ 0 −A˜0 · · · −A˜m−1 −A˜m · · · 0 0 ⎥ ⎢ ⎥ ⎢ .. .. . . .. .. .. .. ⎥ .. ⎣ . . . . . . . . ⎦ 0 0 ··· 0 0 · · · −A˜m−1 −A˜m
¯j = B ⊗ I where A¯j = Aj ⊗ I , A˜j = I ⊗ Aj , B j
˜j = I ⊗ B , j = 0, 1, ..., m. The initial conditions and B j of the above system are not known, but we can find them by solving the following boundary value problem: a first set of boundary conditions follows from the definition of the auxiliary matrices introduced in (5): Xj+1 (0) = Xj (h) , j = −m, ..., 0, ..., m − 2, and a second could be found by expressing the algebraic property (4) in terms of the auxiliary matrices (5): m −W = [Bk X0 (0) Ak + Ak X0 (0) Bk ] k=0
+
m m
Bk Xk−j−1 (h) Aj + Bj Xk−j−1 (h) Ak
j=0 k=j+1
+
m m j=0 k=j+1
Aj Xk−j−1 (h) Bk + Ak Xk−j−1 (h) Bj .
4. FREQUENCY DOMAIN STABILITY ANALYSIS In this section the key relationship between the spectrum of the delay free system (8) and the roots of the characteristic function of system (1) is established, Kharitonov (2006), Louisell (2001). Theorem 1. Let s0 be an eigenvalue of the neutral delay system (1), such that −s0 is also an eigenvalue of system (1). Then s0 belongs to the spectrum of the delay free system (8). Proof. The characteristic equation of system (1) is given by m G (s, h) = sBk e−shk − Ak e−shk .
γ G(s0 , h) = 0 and G (−s0 , h) μ = 0. (9) If a complex number s belongs to the spectrum of the delay free system (8) then there exist a set of 2m non trivial matrices Xj , j = −m, ..., 0, ...m − 1 such that the following holds: m [sXj−k Bk − Xj−k Ak ] =0, j ≥ 0, (10) k=0
m
[sBk Xj+k + Ak Xj+k ] =0,
j < 0.
k=0
Computing e−s0 hj μ [γ G(s0 , h)] = 0, j = 0, ..., m one arrives to the expression m s0 Bk es0 h(j−k) − Ak es0 h(j−k) = 0 μγ
(11)
k=0
for j = 0, ..., m.
Now if we calculate [G (−s0 , h)μ] γ es0 hj = 0, j = −m, ..., −1 it follows that m s0 Bk es0 h(j+k) + Ak es0 h(j+k) μγ = 0.
(12)
k=0
If we define the matrices ¯ l = es0 hl μγ , l = −m, ..., 0, ..., m − 1, X the matrix equations (11)-(12) take the form m m ¯ j−k Bk = ¯ j−k Ak , j ≥ 0, s0 X X k=0
m
¯ j+k = − s0 Bk X
k=0
k=0 m
¯ j+k , j < 0. Ak X
k=0
¯ l are non trivial, then It is clear that the matrices X the value s0 belongs to the delay free system (8). This procedure also holds for −s0 . Proposition 2. The spectrum of the delay free system (8) is symmetrical with respect to the imaginary axis. Proof. The fact that the spectrum of the delay free system of matrix equations (8) is symmetrical with respect of the imaginary axis follows from the idea that if for a value s there exist non trivial matrices Xi i = −m, ..., 0, ...m − 1 that satisfy system (10). By transposing and multiplying by (−1) system (10) it follows that m −sBk Xj−k =0, + Ak Xj−k k=0 m k=0
−sXj+k Bk − Xj+k Ak =0,
and by setting the matrices ˜j = X X −j−1 , j = −m, ..., 0, ..., m − 1 we conclude that the system is satisfied (10) for −s. 5. PROPOSED METHODOLOGY
k=0
Since the values s0 and −s0 belong to the spectrum of the neutral delay system (1), there exist non zero vectors μ and γ, such that
In this section two different methodologies are proposed: The first one applies to systems with an unknown parameter in the matrices of the system. The second methodology
6. EXAMPLES
applies to systems where the only unknown parameter is the delay h. 5.1 Methodology A Here we consider the class of systems that are of the form m m A˜k x (t − hk) . Bk x˙ (t − hk) = k=0
k=0
where A˜k = kμAk , μ is an unknown parameter that satisfies μ > 0. The methodology for the computation of the critical parameters is as follows: (1) Calculate the corresponding delay free system (2) Compute the characteristic polynomial p(s, ρ) of the delay free system (8). The value ρ represent the delay h or the unknown parameter μ (or both) (3) As the spectrum of the system is symmetrical with respect to the imaginary axis, all the powers of the polynomial are even, then we can set λ = s2 and analyze the polynomial p1 (λ, ρ) (4) Determining the purely imaginary roots of p(s, ρ) reduces to find the negative real roots of p1 (λ, ρ) (5) Using Sturm’s Theorem Uspensky (1948), determine the negative real roots - Compute the Sturm sequence of the polynomial p1 (λ, ρ) - Determine with the help of Sturm’s Theorem if p1 (λ, ρ) has negative real roots (6) If there exist a negative real roots compute the corresponding purely imaginary roots and form a set of candidate roots for the delay system (7) For each of these candidate roots, check if there exists a value of the parameter ρ that annihilate the characteristic function. If this is the case then the pair (ωi∗ , ρ∗i ) , receives the name of critical frequency and critical parameter (8) Form a table of ρ∗i in ascending order (consider ρ0 = 0) (9) For ρ∗i , ρ∗i+1 and ωi∗ determine the number of unstable roots of the characteristic function using the argu- ment principle to determine if the interval ρ∗i , ρ∗i+1 is a stable or unstable region. 5.2 Methodology B For this methodology we consider the class of neutral type time delay systems of the form (3), the methodology consists of the following steps: a. The steps (1)-(3) of Methodology A remains valid b. Compute the negative real roots and the corresponding purely imaginary roots c. Make the change of variable z = e−sh d. Evaluate the characteristic function f (s, h) for s = iω obtained from step (b) e. Compute the roots of the polynomial f (z) and select those who satisfies |z| = 1 f. Compute the critical delays associated to the critical frequency ω ∗ , that is −Arg (z) h∗ = ω∗ g. Continue with Steps (8)-(9) of Methodology A for the parameter h.
In this section we illustrate the use of the proposed methodologies for the computation of critical parameters. Example 3. Let us consider the following system studied in Ivanescu et al. (2003), Rodriguez et al. (2004) x˙ (t) + B1 x˙ (t − h) = A0 x (t) + A1 x (t − h) where −2 0 −1 0 , A1 = , A0 = 0 −0.9 −1 −1 −0.1 0 B1 = 0 −0.1
Since ρ (B1 ) = 0.1 < 1 the difference operator is stable. The corresponding delay free system is of the form (I ⊗ I) (B1 ⊗ I) z (τ ) = (I ⊗ B1 ) (I ⊗ I) (A0 ⊗ I) (A1 ⊗ I) z (τ ) . −(I ⊗ A1 ) −(I ⊗ A0 )
In order to find the critical delays, we first compute the characteristic polynomial of the delay free system p (s) = s8 − 5. 689 1s6 + 8. 1628s4 − 0.1955s2 − 0.3797.
By setting λ = s2 , and analyzing the auxiliary polynomial p1 (λ) we can directly check that the only negative real root 19 is λ = − 99 which corresponds to the purely imaginary root √ √ 11 19 s = i 33 . Now, the characteristic function of the delay system is f (s, h) = det(sI + sB1 e−sh − A0 − A1 e−sh ). Evaluating f (s, h) for s = iω and the change of variable z = e−sh ; f (s, h) take the following form f (z) = z 2 (0.998 0 − 0.0087i) + z (2. 938 3 + 0.749 1i) + (1. 608 0 + 1. 270 4i) with roots
z1 = −0.879 1 − 0.4765i, z2 = −1. 977 0 − 0.524 6 i.
Now we are looking for roots z1,2 that satisfies |z| = 1, the only root that satisfies this condition is z1 . The critical delay is computed as −Arg (z1 ) h∗ = = 6.0371. ω∗ Finally with the help of the argument principle we can conclude that the system is stable for h ∈ [0, 6.0371) and is unstable for h ∈ (6.0371, ∞). Example 4. In this example we consider a system with an unknown parameter μ that satisfies μ > 0. The system is of the form x˙ (t) + B1 x˙ (t − 1) = μA0 x (t) + A1 x (t − 1) where A0 , A1 and B1 are the same as in the previous example. The characteristic polynomial of the corresponding delay free system is 1 p(s, μ) = 4 q (s, μ) r (s, μ) v (s, μ) 99 where
Re f(ω1 (μ ), μ ) Im f(ω1 (μ ), μ )
3 2 1 0 −1 −2 −3
0
0.1
0.2
0.3
μ
0.4
0.5
Fig. 1. Real and Imaginary parts of f (iω1 , μ) 8 Re f(ω2 (μ ), μ ) Im f(ω2 (μ ), μ )
6
−5.5779 5.0852 3.5478 A0 = −1.7151 2.2988 0.8782 , 0.4407 3.2284 −1.3174 0 −0.2 0.4 0.3 0.1 0 0 −0.2 0 B1 = 0.5 −0.3 0 , B2 = , −0.2 −0.7 0 −0.1 0 −0.4 The matrices of the delay free system are ⎤ ⎡ 0 I ⊗ I B1 ⊗ I B2 ⊗ I I ⊗ I B1 ⊗ I B2 ⊗ I ⎥ ⎢ 0 M1 = ⎣ I ⊗ B2 I ⊗ B1 I ⊗ I 0 ⎦ 0 I ⊗ B2 I ⊗ B1 I ⊗ I and ⎡ ⎤ 0 0 0 A0 ⊗ I A0 ⊗ I 0 0 ⎢ 0 ⎥ M2 = ⎣ ⎦ 0 0 −I ⊗ A0 0 0 0 0 −I ⊗ A0 The expression of the characteristic polynomial of this delay free system, which is of degree 32, is omitted due to space limitations. With the change of variable λ = s2 the polynomial p1 (λ) is obtained. We can verify that this polynomial has a pair of negative real roots
4
4
λ1 = −0.8807 λ2 = −1.3665 × 10−4
2
with corresponding purely imaginary roots
0
s1 = 0.9384i
−2
s2 = 0.0116i. −4
0
0.2
0.4
μ
0.6
0.8
1
Fig. 2. Real and Imaginary parts of f (iω2 , μ) q (s, μ) = 99s − 81μ + 100 2
2
r (s, μ) = 99s2 − 400μ2 + 100
v (s, μ) = 9801s4 + s2 −47 740μ2 + 19 800
+ 32 400μ4 − 36 000μ2 + 10000
With the change of variable λ = s2 we can directly check that the auxiliary polynomial q1 (λ, μ) has one negative
, while the auxiliary real root for μ ∈ 0, 10 9
polynomial r1 (λ, μ) has one negative root for μ ∈ 0, 12 . With the help of the Sturm array we can conclude that the auxiliary polynomial v1 (λ, μ) does not have negative real roots ∀μ > 0. The characteristic function of the system is of the form
f (s, μ) = det sI + sB1 e−s − μA0 − A1 e−s .
We must check if there exist pairs {iωk (μ) , μ} k = 1, 2 that annihilate the characteristic function of the system. By direct calculations we can prove that for μ1 = 0.4615 and ω1∗ = 0.3867, see Figure 1, and for ω2∗ = 0.3866 and μ2 = 1.0256, see Figure 2, the real and imaginary parts vanish simultaneously Example 5. Now let us consider the system introduced in Michiels and Vyhl´ıdal (2005), Jarlebring (2007), x˙ (t) + B1 x˙ (t − h) + B2 x˙ (t − 2h) = A0 x (t) , where
Now the characteristic function of the system is of the form
f (s, h) = det sI + sB1 e−sh + sB2 e−2sh − A0 .
With the change of variable z = e−sh and considering s1 = 0.9384i, f (s, h) takes the form f (z) = − (0.0198i) z 6 − (0.0396i) z 5 − (0.6969 − 0.1388i)z 4 + (0.6452 + 0.119i) z 3 + (2.9265 − 0.9274 i) z 2 + (0.2915 +6.5017i) z − (3.994 6 + 4.748 9i) The only root zi that satisfies |z| = 1 is z1 = 0.9839 − 0.1785i, then the critical delay is Arg(z1 ) h∗1 = − = 0.1912. ω1∗ If we compute f (s, z) for ω2 = 0.0116 we obtain the polynomial
f (z) = − 7.667 × 10−8 i z 6 − 1.533 × 10−7 i z 5
− 2.162 × 10−4 − 5.367 × 10−7 i z 4
+ 0.0536 + 2.3 × 10−7 i z 3
+ 4.5408 × 10−4 − 0.12788i z 2
+ 4.523 × 10−5 + 0.0779 × 10−2 i z + (0.0523 − 0.09i)
that has roots z1 = −9.6132 + 7.6772i z3 = −39.8206 − 38.2932i z2 = 0.5582 + 0.8296i z4 = 36.811 7 + 38.296 0i z3 = 5.301 2 + 0.07.574i z5 = 4.762 5 − 8.585 4i
The only root that satisfies |z| = 1 is z2 then Arg(z2 ) 2π + ∗ = 453.7893. h∗2 = − ω2∗ ω2 The argument principle shows that the only interval where the system is stable is for h ∈ [0, 0.1912). 7. CONCLUSIONS A procedure for the computation of the critical parameters of time delay systems of neutral type is presented. The approach takes advantage of the key relation between the roots of the characteristic polynomial of a delay free system of matrix equations and the roots of the characteristic function of the neutral type time delay system. Some illustrative examples show the effectiveness of the method for determining the critical delays as well as other critical parameters. REFERENCES Bellman, R. and Cooke, K.L. (1963). DifferentialDifference Equations. Academic Press, New York. Gohberg, I.C. and Lerer, L.E. (1976). Resultants of matrix polynomials. Bulletin of the American Mathematical Society, 82, 565–567. Ivanescu, D., Niculescu, S.I., Dugard, L., Dion, J.M., and Verriest, E.I. (2003). On delay-dependent stability for linear neutral systems,. Automatica, 39, 255–261. Jarlebring, E. (2007). On critical delays for linear neutral delay systems. European Control Conference, Kos, Greece. Kharitonov, V.L. (2005). Lyapunov functionals and lyapunov matrices for neutral type time delay systems: a single delay case. International Journal of Control, 78, 783–800. Kharitonov, V.L. (2006). Lyapunov matrices for a class of time delay systems. Systems and Control Letters, 55, 610–617. Kharitonov, V.L. (2009). Lyapunov matrices: Existence and uniqueness issues. 8th Workshop on Time Delay Systems, Sinaia, Romania. Louisell, J. (1998). Numerics of the stability exponent and eigenvalue abscissas of a matrix delay system. In L. Dugard and E.I. Verriest (eds.), In Lecture Notes in Control and information sciences 228, volume 228. Springer. Louisell, J. (2001). A matrix method for determining the imaginary axis eigenvalues of a delay system. IEEE Transactions on Automatic Control, 46, 2008–2012. Michiels, W. and Vyhl´ıdal, T. (2005). An eigenvalue based approach for the stabilization of linear time-delay systems of neutral type. Automatica, 41, 991–998. Ochoa, G., Vel´ azquez, J.E.V., Kharitonov, V.L., and Mondi´e, S. (2009). Lyapunov matrices for neutral type time delay systems. In J. Loiseau, W. Michiels, S.I. Niculescu, and R. Sipahi (eds.), In Lecture Notes in Control and information sciences 388, volume 388. Springer. Olgac, N. and Sipahi, R. (2004). A practical method for analyzing the stability of neutral type lti-time delayed systems. Automatica, 40, 847–853. Rodriguez, S.A., Kharitonov, V.L., Dion, J.M., and Dugard, L. (2004). Robust stability of neutral systems:
a lyapunov krasovskii constructive approach. International Journal of Robust and Nonlinear Control, 14, 1345–1358. Uspensky, J.V. (1948). Theory of Equations. McGraw-Hill, New York. Walton, K.E. and Marshall, J.E. (1987). Direct method for tds stability analysis. IEEE proceedings, 134, 101–107.