Critical frequencies and parameters for linear delay systems: A Lyapunov matrix approach

Systems & Control Letters 62 (2013) 781–790

Contents lists available at SciVerse ScienceDirect

Systems & Control Letters journal homepage: www.elsevier.com/locate/sysconle

Critical frequencies and parameters for linear delay systems: A Lyapunov matrix approach Gilberto Ochoa a , Vladimir L. Kharitonov b , Sabine Mondié a,∗ a

Department of Automatic Control, CINVESTAV-IPN, México D.F., Mexico

b

Faculty of Applied Mathematics and Control Processes, St.-Petersburg State University, St.-Petersburg, Russia

article

info

Article history: Received 27 June 2012 Received in revised form 24 January 2013 Accepted 26 May 2013 Available online 5 July 2013 Keywords: Time delay systems Distributed delays Neutral systems Stability Critical delays

abstract A new method, based on recent results concerning the construction of the Lyapunov matrix of delay systems, reveals conditions under which the characteristic function has a root s0 such that −s0 is also a root. As the pure imaginary roots, which are crucial in the stability analysis of time delay systems, are of that type, this method gives rise to a novel approach for the delay independent and delay dependent stability analyses of systems with delay that are multiple of a basic delay, distributed delays and of a class of neutral type time delay systems. A number of examples are given to illustrate the approach and to show its strength. © 2013 Elsevier B.V. All rights reserved.

1. Introduction It is well known that a linear time invariant delay system of the form dx(t )

= A0 x(t ) + A1 x(t − h), (1) dt where A0 , A1 ∈ Rn×n , and h is a positive delay, is exponentially stable if and only if all roots of its characteristic quasipolynomial g (s, h) = det(sI − A0 − A1 e−hs ),

(2)

have negative real parts [1]. In the case when the system matrices depend continuously on some parameters and/or the delay h is not fixed one may try to detect the values of the parameters for which the system may loose its stability. For such critical values the following should occur: at least one of the roots of g (s, h) comes across the imaginary axis of the complex plane. There is a great variety of algorithms for the computation of the critical parameter values. A thorough study of the algorithms reveals that they have a common core: all of them are based on some conditions under which the characteristic quasipolynomial (2) has a root s0 such that −s0 is also a root of g (s, h). Once these roots are detected one has to select among them those with the real part equal to zero. The pure imaginary roots constitute the set

of critical frequencies, and the corresponding values of the system parameters form the set of critical values of the parameters. To illustrate this point let us give a brief tour of some of the most popular algorithms. The first set of algorithms is based on the observation that the spectrum of the Kronecker sum of n × n matrices A and B, A⊕B = A⊗I +I ⊗B, where I is the identity matrix, consists of all sums of the form λ + µ, where λ belongs to the spectrum of matrix A, and µ belongs to the spectrum of matrix B; see [2]. Now, let us assume that s0 and −s0 are roots of g (s, h); then s0 belongs to the spectrum of matrix (A0 + zA1 ), where z = e−hs0 , while −s0 belongs to the spectrum of matrix (A0 + z −1 A1 ). As a result, the determinant of the Kronecker sum (A0 + zA1 ) ⊕ (A0 + z −1 A1 ) is equal to zero. A rewriting of the Kronecker sum in the matrix pencil form reduces the problem to the existence on the unit circle of the complex plane of generalized eigenvalues of the corresponding matrix pencil. The resulting algorithms are known as constant matrix tests [3], or matrix pencil techniques [4–6]. The origin of another set of algorithms for the computation of the critical values of the system parameters lies in the next observation: it is not difficult to find m ≥ 0 such that the following two functions are bivariate polynomials g1 (s, z ) = det(sI − A0 − zA1 ) and g2 (s, z ) = z m det(sI + A0 + z −1 A1 ).

∗

Corresponding author. Tel.: +52 555061 3738. E-mail addresses: [email protected] (G. Ochoa), [email protected] (V.L. Kharitonov), [email protected] (S. Mondié). 0167-6911/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sysconle.2013.05.010

Once again, assume that s0 and −s0 are roots of g (s, h); then, by fixing the first variable s = s0 , and considering the resulting polynomials in the variable z, one can conclude that they have a

782

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

common root z0 = e−hs0 . It is well known that two polynomials have a common root if and only if the resultant matrix of the polynomials is singular; see [2]. Hence, the critical values of the system parameters are those for which the determinant of the resultant matrix, as a function of s, admits roots on the imaginary axis of the complex plane. The set of such roots constitutes the set of critical frequencies of the time delay system. There are several computational schemes based on this observation; see [3,7]. The Ts Rekasius substitution, z = 11− , where T is a real parameter, +Ts provides an interesting variation of this approach. In this case, after multiplication by an appropriate power of (1 + Ts), the function

 det sI − A0 −

1 − Ts 1 + Ts

 A1

,

turns into a polynomial g (s, T ) in the variable s with coefficients depending on the parameter T . It is just a matter of simple calculations to check that a pure imaginary root, s0 = iω0 , of the quasipolynomial g (s, h) is also a root of g (s, T ) for a special choice of T . As in this case the real polynomial g (s, T ) has a root s0 such that −s0 is also a root of the polynomial, and it admits an ‘‘even–odd’’ decomposition g (s, T ) = α(s2 , T ) + sβ(s2 , T ), where α(s2 , T ) and β(s2 , T ) have a common root s20 . To detect the system parameters for which α(s2 , T ) and β(s2 , T ) admit a common root one can employ, once again, the resultant matrix of the polynomials. It is worthy of mention here that up to the order of the rows this new resultant matrix coincides with the Hurwitz matrix of the polynomial g (s, T ); see for details [8,9]. A different algorithm for the computation of the critical values of the system parameters has been presented in [10,11]. It is based on the fact that any pure imaginary root of the characteristic function g (s, h) is also a root of the polynomial p(s) = det (sI + AT0 ) ⊗ (sI − A0 ) − AT1 ⊗ A1







,

which, in fact, is the characteristic polynomial of the following delay free system of two matrix equations [12]



X0′ (τ ) = X0 (τ )A0 + X−1 (τ )A1 ′ T T X− 1 (τ ) = −A1 X0 (τ ) − A0 X−1 (τ ),

(3)

that is used for the computation of the Lyapunov matrices of system (1). This explicit relation between the spectrum of the original time delay system (1) and that of the delay free system (3), which establishes a bridge between the frequency domain stability analysis of time delay systems and that in the time domain, has not been thoroughly studied, and no attempts have been made to propose similar algorithms for the computation of the critical values of the system parameters for more general classes of time delay systems. We first illustrate with the simple scalar time delay equation the potential advantages of this approach. Example 1. Consider the scalar delay equation x˙ (t ) + ax(t ) + bx(t − h) = 0, with characteristic quasipolynomial g (s, h) = s + a + be−sh . The delay free system (3) is now of the form



x0 (τ ) = −ax0 (τ ) − bx−1 (τ )

√

has a pair of pure imaginary roots ±iω = ±i b2 − a2 . Substituting the roots into the characteristic quasipolynomial and solving the resulting equation for h one arrives at the following critical values of the delay: h=

arccos(−a/b) + 2kπ

√

b 2 − a2

= bx0 (τ ) + ax−1 (τ ),

and its characteristic polynomial p(s) = s2 + b2 − a2 . In the case b2 − a2 < 0 the polynomial has no imaginary roots, so in this case the quasipolynomial g (s, h) has no zeros on the imaginary axis for any h ≥ 0. If, additionally, a + b > 0, then the delay equation is delay independent stable. In the case b2 − a2 > 0 the polynomial

k = 0, ±1, ±2, . . . .

This approach can be extended with no difficulty, to the case of vector systems (1). The aim of this contribution is to exploit similarly the results on Lyapunov matrices for multiple delays systems, distributed delays systems [13–15] and neutral type time delay systems [16,17]. Furthermore, as the obtained delay free systems are of counterflow type [18], their characteristic polynomial p(s) only has even powers of s. As a consequence, the original problem of finding pure imaginary roots of p(s), reduces to determine the nonpositive real roots of polynomial p1 (λ), with λ = s2 , a problem that becomes non trivial as the dimension of the system and number of delays increases. For tackling it, we choose to use Sturm’s Theorem [19], which is based on the number of sign changes of Sturm’s sequence in λ ∈ (−∞, 0). We would like to stress out that an other advantage of this choice is that it allows the analysis of systems that depend on some unknown parameter. Motivated by the above discussion and example, we start Section 2 with an extension of the connection between the spectrum of the characteristic equation of systems with delays that are multiple of a basic one, that also present a distributed delay, and that of a special delay free system of matrix equations arising in the computation of Lyapunov matrices for time delay systems. In Section 3, an extension of these results to a class of neutral type time delay systems delay is proposed. Then a general methodology of how to use these results for the computation of the critical values of system parameters is explained in Section 4. Several illustrative examples that validate the proposed scheme, and show the strength of this approach are given in Section 5. The contribution ends with some concluding remarks. Notation: Throughout this paper, Rn×n is the set of all n×n matrices, the superscript ‘‘T ’’ stands for the transpose matrix. For a given z ∈ C the module is denoted by |z | and its argument by Arg(z ). We denote by U ′ (τ ) the derivative of U (τ ) with respect of τ . 2. Systems with concentrated and distributed delays In this section, we consider time delay systems of the form dx(t ) dt

=

m 

Ak x(t − hk) +



0

H (θ )x(t + θ )dθ ,

(4)

−hm

k=0

where Ak ∈ Rn×n , h > 0 is the basic delay, H (θ ) is a continuous real-valued matrix function that can be expressed as H (θ ) =

¯ m 

ηk (θ ) Hk .

(5)

k=0

Here Hk ∈ Rn×n and the scalar functions ηj (θ ) satisfy d dθ

ηj (θ ) =

¯ m 

γjk ηk (θ) ,

(6)

k=0

and

′

x′−1

,

ηj (−θ − hm) =

¯ m 

αj,k ηk (θ ) .

(7)

k=0

2.1. Delay free system A key step in this new approach to the stability analysis of linear time-delay systems, is the computation of a delay free system

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

of matrix differential equations, that is closely related with the Lyapunov matrix. To this end, let us first introduce the Lyapunov matrix associated to system (4). Definition 2. Given a matrix W = W T . Then, the solution of the equation U ′ (τ ) =

m 

U (τ − hk) Ak +

0



U (τ + θ) H (θ ) dθ ,

(8)

−hm

k=0

U T (τ ) = U (−τ ) ,

(9)

and

− W = AT0 U (0) + U (0)A0 +

m  

ATk U (hk) + U T (hk)Ak



Xj′ (τ ) =

m 

0

H T (θ) U T (θ ) dθ +



Xj′ (τ ) = −

U (θ ) H (θ) dθ

(10)

k = −m, . . . , 0, . . . , m − 1

(11)

¯ + 1) matrices • m (m   0 j = −m, . . . , −1 Zj,k (τ ) = ηk (θ) Xj (τ − θ ) dθ , . . . ¯ . (13) k = 0, . . . , m −hm Computing the derivative with respect to τ of the auxiliary matrices (11)–(13), and considering that they satisfy the Lyapunov matrix properties (8)–(10), we arrive to the following delay free system of matrix equations Xj−k (τ ) Ak +

¯ m 

Yj,k (τ ) Hk

{j = 0, . . . , m − 1

ATk Xj+k (τ ) −

k=0



HkT Zj,k (τ )

k=0

{j = −m, . . . , −1 Yj′,k (τ ) = ηk (0)Xj (τ ) − ηk (−hm) Xj−m (τ )

−

¯ m 

γk,l Yj,l (τ )



l =0

(14)

j = 0, . . . , m − 1 ¯ k = 0, . . . , m

Zj′,k (τ ) = −ηk (0)Xj (τ ) + ηk (−hm) Xj+m (τ )

+

¯ m  l=0

HkT Zj,k (τ )

k=0

(16)

Φ (τ ) = X¯ (τ ) (η(0) ⊗ I ) − Xˆ (τ ) (η (−hm) ⊗ I )   − Φ (τ ) Γ T ⊗ I Ψ ′ (τ ) = X¯ (τ ) (η (−hm) ⊗ I )   − Xˆ (τ ) (η(0) ⊗ I ) + Ψ (τ ) Γ T ⊗ I .

Yj,k (τ ) =

¯ m 

αk,l Zj−m,l (τ )

and

Zj,k (τ ) =

¯ m 

αk,l Yj+m,l (τ ) .

γk,l Zj,l (τ )

Once these relations are established, one may reduce the delay free system (16) to the system Xj′ (τ ) =

m 

Xj−k (τ ) Ak +

k =0

Xj′ (τ ) = −

¯ m 

Yj,k (τ ) Hk

{j = 0, . . . , m − 1

k =0

m 

ATk Xj+k (τ ) −

¯ m  m 

k =0

αk,l HkT Yj+m,l (τ )

k=0 l=0

(17)

{j = −m, . . . , −1 Φ (τ ) = X¯ (τ ) (η(0) ⊗ I ) − Xˆ (η (−hm) ⊗ I )   − Φ (τ ) Γ T ⊗ I . ′

The following result links the spectrum of the delay system (4) with that of the delay free system (17). This relation can be exploited in order to compute the critical delays, frequencies or parameters values.

¯ m



¯ m 

2.2. Frequency stability analysis

k=0 m

Xj′ (τ ) = −

ATk Xj+k (τ ) −

l =0

¯ + 1) matrices of the form • m (m   0 j = 0, . . . , m − 1 Yj,k (τ ) = ηk (θ) Xj (τ + θ ) dθ , . . . (12) ¯ k = 0, . . . , m −hm

k=0

m  k =0

• 2m auxiliary matrices

Xj (τ ) =

{j = 0, . . . , m − 1

l=0

Now, let us introduce the following sets of auxiliary matrices

m 

Yj,k (τ ) Hk

k=0

{j = −m, . . . , −1

is called the Lyapunov matrix of system (4) associated with W .

′

¯ m 

k =0

0

−hm

−hm

Xk (τ ) = U (τ + kh) ,

Xj−k (τ ) Ak +

Now, in view of (7) matrices Yi,k (τ ) and Zj,l (τ ) satisfy the relations

k=1

+

the delay free system (14) can be expressed as

′

that satisfies properties



783

Theorem 3. Let s0 be an eigenvalue of the time delay system (4) such that −s0 is also an eigenvalue of the system. Then the value s0 belongs to the spectrum of the delay free system of matrix equations (17). Furthermore, the spectrum of the delay free system is symmetrical with respect to the imaginary axis. Proof. The characteristic matrix of system (4) is

j = −m, . . . , −1 ¯. k = 0, . . . , m



G (s, h) = sI −

If we define

 m¯ ,m¯  m¯ ,m¯  m¯ ,m¯ Γ = γi,j i,j=0 , ∆ = αi,j i,j=0 , Φ (τ ) = Yi,j (τ ) i,j=0 ,  m¯ ,m¯ Ψ (τ ) = Zi,j (τ ) i,j=0 ,  T X¯ (τ ) = X0T (τ )X1T (τ ) · · · XmT −1 (τ ) ,  T T T T Xˆ (τ ) = X− , m (τ )X−(m−1) (τ ) · · · X−1 (τ )   η (θ ) = η0 (θ) η1 (θ) · · · ηm¯ (θ )

m 

Ak e−shk −

k=0

¯ m 

f (k) (s)Hk ,

k=0

where (15)

f (k) (s) =



0

ηk (θ ) esθ dθ .

−hm

If s0 and −s0 are both eigenvalues of system (4) then, there exist nonzero vectors µ and δ such that

δ T G (s0 , h) = 0,

GT (−s0 , h) µ = 0.

784

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

Premultiplying the first of the above equalities by the factor

µes0 hj j = 0, . . . , m − 1, and postmultiplying the second one by δ T es0 hj j = −m, . . . , −1 we obtain the equations m 

s0 es0 hj µδ T − µδ T

− µδ T es0 hj

¯ m 

s0 f (k) (s0 ) +

γk,l f (l) (s0 ) = ηk (0) − ηk (−hm) e−s0 hm .

l =0

es0 h(j−k) Ak

Premultiplying the above equality by es0 hj µδ T

k=0

¯ m 

Comparing the right hand side of both equalities and evaluating at s = s0 , yield

s0 f (k) (s0 ) es0 hj µδ T = ηk (0)es0 hj µδ T − ηk (−hm) es0 h(j−m) µδ T

f (k) (s0 ) Hk = 0

(18)

−

k=0

−s0 µδ T es0 hj −

m 

¯ m 

γk,l f (l) (s0 ) es0 hj µδ T .

l =0

ATk µδ T es0 h(j+k)

Then, it follows from the matrix definitions (21) and (15) that

k=0

−

¯ m 

f

(k)

(− )

s0 HkT

µδ e

T s0 hj

= 0.

(19)

k=0

s0 Yj,k = ηk (0)Xj − ηk (−hm) Xj−m −



Since ηj (θ ) satisfy conditions (3) and (7) the following relation

f (j) (−s) = eshm

¯ m 

αj,k f (k) (s),

k=0

m 

ATk XjT−k −

k=0

− s0 µδ T es0 hj −

m 

−s0 XjT =

ATk µδ T es0 h(j+k)

¯  ¯ m m 

α

T k,l Hk

T s0 h(j+m) (l)

µδ e

f

(s0 ) = 0.

(20)

k=0 l=0

Yj,k = µδ T es0 hj f (k) (s0 ) ,

(21)

Xj−k Ak +

k=0

¯ m 

s0 X j = −

 k=0

Yj,k Hk

0

−hm

d dθ

0

−hm

d dθ

and

¯ m 

αk.l Y˘−(j+1)+m,l = YjT,k ,

l =0

and using the fact that ΛΓ ηT (θ ) = −Γ ΛηT (θ ).



αk,l HkT Yj+m,l .

k=0 l=0

ηk (θ) e dθ = s



0



¯ m 

ηk (θ ) esθ dθ

l =0

γk,l



Remark 4. It is worthy of mention that for retarded systems of the form

dt

−hm

+

(23)

dx(t )

ηk (θ) esθ dθ = ηk (0) − ηk (−hm) e−shm ,

sθ

αk,l Yj+m,l Hk

k=0 l=0

¯ m

and on the other hand



XjT+k Ak +

  −s0 (Λ ⊗ I ) Φ T = − (Λ ⊗ I ) ηT (0) ⊗ I X¯ T   + (Λ ⊗ I ) ηT (−hm) ⊗ I Xˆ T

(22)

Now, on the one hand we have that



k=0

¯  ¯ m m 

−s0 Φ T = − (η(0) ⊗ I )T X¯ T + (η (−hm) ⊗ I )T Xˆ T  T + Γ T ⊗ I ΦT .

X˘ −(j+1) = XjT

¯ m

ATk Xj+k −

HkT YjT,k

Finally, the result follows from defining the matrices

k=0 m

¯ m 

+ (Λ ⊗ I ) (Γ ⊗ I ) Φ T .

Eqs. (18), (20) can be expressed as m 

(24)

By premultiplying the last equation by the factor (Λ ⊗ I ) from the left, it follows that

By defining the matrices Xj = µδ T es0 hj ,

m  k=0

k=0



Collecting Eqs. (22)–(24), as the matrices Xj , Yj,k are non trivial, we conclude that the value s0 belongs to the spectrum of the delay free system (17). The previous analysis and conclusion remain valid for the value −s0 . The fact that the spectrum of the delay free system is symmetrical with respect to the imaginary axis, can be proved transposing and negating system (22)–(24), that is

−s0 XjT = −

therefore, (19) takes the form

s0 X j =

γk,l Yj,l

l=0

s0 Φ = X¯ (η(0) ⊗ I ) − Xˆ (η (−hm) ⊗ I ) − Φ Γ T ⊗ I .

holds

−

¯ m 

0

−hm

esθ ηl (θ) dθ .

=

m 

Ak x (t − hk) ,

(25)

k=0

neither the delay free system (17) nor its spectrum depends on the delay value h. On the other hand, the spectrum of the system (25) depends on the delay. Remark 5. Because of the symmetry with respect to the imaginary axis of the spectrum of the delay free system, the characteristic function p (s, h) can be written as p1 (λ, h) where λ = s2 . A numerically important consequence of this fact is that the problem of determining the purely imaginary roots of p (s, h) reduces to finding the nonpositive real roots of the real polynomial p1 (λ, h). This in turn implies that the maximum number of critical frequencies ¯ + 1). is 2n2 m (m

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

3. Neutral type time delay systems In this section we analyze the class of neutral type systems of the form m 

Bk

dx (t − hk) dt

k=0

m 

=

Ak x (t − hk) ,

t ≥ 0.

(26)

• The dynamic property m  

U ′ (τ − hk) Bk − U (τ − hk) Ak = 0,



τ ≥ 0.

• The symmetry property (28)

• The algebraic property (29)

j=0 k=0

Let us recall the definition of the auxiliary matrices (11) j = −m, . . . , 0, . . . , m − 1.

Computing the derivatives of matrices Xi (τ ) and using the Lyapunov matrix properties (27)–(29), lead to the delay free system of matrix equations

k=0

m 



Xj−k (τ ) Ak ,

k=0



m  

s0 BTk es0 h(j+k) + ATk es0 h(j+k) µγ T = 0.



(34)

k=0

If we define the matrices X¯ l = es0 hl µγ | ,

l = −m, . . . , 0, . . . , m − 1,

m 

s0 X¯ j−k Bk =

m 

m 

X¯ j−k Ak ,

j ≥ 0,

k=0

s0 BTk X¯ j+k = −

m 

ATk X¯ j+k ,

j < 0.

k=0

It is clear that the matrices X¯ l are non trivial, then the value s0 belongs to the spectrum of the delay free system (30). This procedure also holds for −s0 . The fact that the spectrum of the delay free system of matrix equations (30) is symmetrical with respect to 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 (32). Transposing and multiplying by (−1) system (32) implies that

k=0

(30)

m

BTk Xj′+k (τ ) = −

(33)

m    −sBTk XjT−k + ATk XjT−k = 0,

j = 0, . . . , m − 1

k=0

m



k=0

BTj U (h (j − k)) Ak + ATk U T (h (j − k)) Bj = −W .

Xj′−k (τ ) Bk =

s0 Bk es0 h(j−k) − Ak es0 h(j−k) = 0,

k=0

τ ≥ 0.

Xj (τ ) = U (τ + jh) ,

m  

the matrix equations (33)–(34) take the form (27)

k=0

m 

µγ T

k =0

Here Ak , Bk ∈ Rn×n are given matrices and h > 0 is the basic delay. B0 is equal to the identity matrix and the difference equation,  m k=0 Bk x (t − hk) = 0 is assumed to be strongly stable [20]. In this case, for a given matrix W = W T , the Lyapunov matrix [21] satisfies the following properties.

m  m 

expressions

k=0

U (τ ) = U T (−τ ) ,

785

 If we compute GT (−s0 , h)µ γ T es0 hj = 0, for j = −m, . . . , −1   and es0 hj µ γ T G(s0 , h) = 0, for j = 0, . . . , m one arrives to the 

ATk Xj+k (τ ) ,

j = −m, . . . , −1.

m    −sXjT+k Bk − XjT+k Ak = 0, k=0

k=0

and the result follows by setting 3.1. Frequency domain stability analysis

T X˜ j = X− j −1 ,

In this section a key relationship between the spectrum of the delay free system (30) and the roots of the characteristic function of system (26) is established. Theorem 6. Let s0 be an eigenvalue of the neutral delay system (26), such that −s0 is also an eigenvalue of system (26). Then, s0 belongs to the spectrum of the delay free system (30). Moreover, the spectrum of the delay free system (30) is symmetrical with respect to the imaginary axis. Proof. The characteristic equation of system (26) is given by G (s, h) =

m  



k=0

γ G(s0 , h) = 0 and G (−s0 , h) µ = 0. T

T

sXj−k Bk − Xj−k Ak = 0,

j ≥ 0,

sBTk Xj+k + ATk Xj+k = 0,

j < 0.



k=0

k=0





0

F (θ ) x (t + θ ) dθ ,

(35)

−h

can also be associated to the spectrum of a delay free system. Indeed, differentiation of (35) leads to a neutral delay system with distributed delay of the form

−A

dx (t − h) dt

= F (0) x(t ) − F (−h) x (t − h) 

0

F˙ (θ ) x (t + θ ) dθ .

−

(36)

−h

(31)

If a complex number s belongs to the spectrum of the delay free system (30) then there exists a set of 2m non trivial matrices Xj , j = −m, . . . , 0, . . . , m − 1 such that the following holds:

m  

x(t ) = Ax (t − h) +

dt

Since s0 and −s0 belong to the spectrum of the neutral delay system (26), there exist nonzero vectors µ and γ , such that

m  

Remark 7. The spectrum of continuous time difference equations of the form

dx(t )

sBk e−shk − Ak e−shk .

j = −m, . . . , 0, . . . , m − 1. 

(32)

The stability of systems (35) and (36) is not equivalent because of the zero dynamics introduced by the differentiation [22,23]. Nevertheless it is clear that the computation of critical frequencies and parameters of (35) can be developed by using the delay free system associated to system (36) reported in [21]. For the particular cases when A = 0 or F (θ) = F0 , differentiation of (35) leads to systems of the form (4) or (25), respectively. Remark 8. The considerations of Remark 5 remain valid. In this case, the maximum number of critical frequencies is 2mn2 .

786

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

4. Methodology In the present section a methodology for the computation of critical delays and critical parameters is introduced. Here, we consider the class of systems of the form dx(t ) dt

=

m 

A˜ k x(t − hk) +

0



H (θ )x(t + θ )dθ , −hm

k=0

or

0 < δ1 < δ2 < · · · < δm . Its characteristic quasipolynomial is 2 f (s) = q(s) + β e−sh , where q(s) = (s2 + δ12 )(s2 + δ22 ) . . . (s2 + δm ). The equation can be written in the form dx(t ) dt

= A0 x(t ) + A1 x(t − h),

where A0 and A1 are (2m)×(2m) matrices. The characteristic polynomial of the corresponding delay free system for the Lyapunov matrix is p(s) = [q(s)]4m−2 q2 (s) − β 2 .



m 

B˜ k

k=0

dx (t − hk) dt

=

m 

A˜ k x (t − hk) ,

k=0

where A˜ k = µk Ak and B˜ k = µk Bk . The steps of the methodology are as follows. 1. Calculate the corresponding delay free system. 2. Compute the characteristic polynomial p (s, δ) of the delay free system, the value δ represent the delay h or the unknown parameter µ, which could be a control parameter. 3. Set λ = s2 and analyze the auxiliary polynomial p1 (λ, δ) of Remarks 5 and 8, looking for nonpositive real roots instead of pure imaginary ones. 4. Sturm’s Theorem (see Appendix), is instrumental in finding those negative real roots, as explained in the three steps described below. (a) Compute Sturm’s sequence of the auxiliary polynomial p1 (λ, δ). (b) Determine if the auxiliary polynomial has negative real roots by applying Sturm’s Theorem. (c) Compute the corresponding purely imaginary roots and form a set of candidate roots for the delay system. 5. Check, for each one of the candidate roots, if there exists a pair (s, δ) that annihilate the characteristic function. If the answer is positive, the value s∗ = iω∗ receives the name of critical frequency and δ ∗ critical parameter. Remark 9. When the only unknown parameter is the delay h and the delays are concentrated (retarded or neutral type system) it is possible to compute analytically the critical delays and values of frequencies. Steps 1–3 of the methodology remain valid, while steps 4 and 5 are replaced by the following. Make the change of variable z = e−sh and evaluate the characteristic function g (s, z ) at each of the candidate roots s = iω, now the polynomial depends only on the variable z. Compute the roots of the polynomial g (z ), and select those that satisfy |z | = 1, then s∗ = iω∗ is a critical frequency and h∗ = −Arg(z )/ω∗ is the critical delay. 5. Examples

Example 10 ([24,25]). Let us consider an oscillatory system with delayed feedback with gain β , of the form d

2

dt 2

+ δ12



Observe that for the case β ̸= 0 the pure imaginary roots of the first factor of the characteristic polynomial p(s) are not roots of the quasipolynomial f (s). So, only the pure imaginary roots of the second factor may participate in the computation of the critical values of the delay h. Let us consider a simple oscillator d2 dt 2

z (t ) + z (t ) + β z (t − h) = 0,

with characteristic quasipolynomial g (s, h) = s2 + 1 − β e−sh .

(37)

The characteristic polynomial of the corresponding delay free system (3) is p(s, β) = q(s)2 r (s), where q(s) = s2 + 1, and r (s) = q(s)2 − β 2 . Notice that the roots of r (s) may produce critical delay values. The critical frequencies for β ∈ (−1, 1) are s∗1,2 =

√ √ ±iω1 = ±i β + 1 and s∗3,4 = ±iω2 = ±i 1 − β . The corre(2k+1)π π sponding critical delays are h∗1,k = , and h∗2,k = 2k , for ω1 ω2 k = 0, ±1, ±2, . . . . The stability regions for those critical values

are depicted in Fig. 1. Regions I, IV, VI, X and XIV correspond to stability regions, while the remaining correspond to unstable ones. Let us consider the case choice of β the critical β = 0.1. For this frequencies are s∗1,2 = ±i

11 10

and s∗3,4 = ±i

9 . 10

The critical delay

values associated to ω1 are 14.9769 and for ω2 are h∗2,1 = can determine the stability of each interval of interest by testing the stability of a point in this interval, using for instance the Argument Principle. From this analysis we can conclude that for h = 0 the system is unstable and for h ∈ (0, 2.9953) the system is stable. For h ∈ (2.9953, 6.6230) the system is once again unstable, and for h ∈ (6.6230, 8.9861) the system enters to a stable region, then for h ∈ (8.9861, 13.2461) the system enters to an unstable one. While, for (13.2461, 14.9769) the system is stable. Notice that these results coincide with those given in [24,25].

= 2.9953, h∗1,1 = 8.9861, h∗1,2 = 6.6230, h∗2,2 = 13.2461. Finally, one

h∗1,0

Next, we revisit the numerical case study proposed in [9] where different methods for determining the imaginary axis crossing frequencies are compared. Example 11. Let us consider the delay system

In this section a number of examples are analyzed to show the strength and versatility of our approach. Although the present paper is focused on critical frequency and parameters, for completeness, the stability properties of some of the examples are described. After appropriate ordering of the critical delays or parameter, the stability of each interval is obtained by testing the stability of a point, with the help of the argument principle.





d

2

dt 2

+ β z (t − h) = 0,



+ δ22 . . .



2

d

dt 2

+ δm2



z (t )

dx(t ) dt

= A0 x(t ) + A1 x (t − h) ,

with

 −1 A 0 = −3 −2

13.5 −1 −1

 −1 −2 , −4

 A1 =

−5.9 2 2

7.1 −1 0

−70.3 5 6

 .

The corresponding delay free system (3) has the characteristic polynomial p(s) = s18 + 5.81s16 − 18 682s14 − 2.4869 × 105 s12

+ 35 185s10 + 4.1197 × 109 s8 + 9.4110 × 1010 s6 + 7.1137 × 1011 s4 + 1.8931 × 1012 s2 + 1.0144 × 1012 .

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

787

Fig. 2. Stability regions.

It follows that the delay independent stability region is the region I on Fig. 3. For a2√ + 2b < 0 the second factor has pure imaginary roots ±iω = ±i −a2 − 2b. The substitution of these roots into g (s, h) gives the condition

Fig. 1. Stability regions. Table 1 (a) Critical frequencies and critical parameters. (b) Stability of the intervals. (a)

(b)

ω

h

h

Stability

15.5032 3.0351 2.9123 2.1109 0.8404

0.22198 0.1624 0.1859 0.8725 7.208

[0.0.1624) (0.1624, 0.1859) (0.1859, 0.2219) (0.2219, 0.6272) (0.6272, ∞)

Stable Unstable Stable Unstable Unstable

Taking λ = s2 , the negative real roots of p1 (λ) are −240.3497, −9.2114, −8.4820, −4.4562, −0.70635. They produce the pure imaginary roots of p(s): ±i15.5032, ±i3.0351, ±i2.9123, ±i2.1109, ±i0.8404. Next, one can proceed as in [8] and compute for each of these roots the delays h for which the quasipolynomial (2) vanishes. In Table 1(a) the set of smallest positive critical delays generated by each pure imaginary root is listed and in Table 1(b) the stability property of each interval is presented. Remark 12. According to the conclusions made by the authors of [9] on the different methods for the determination of imaginary axis crossing roots, the best one appears to be the one based on the Rekasius substitution, which reduces to the determination of the negative real roots of a polynomial with real coefficients of degree n2 . Notice that our approach reduces to a problem with same numerical difficulty. In the next example we consider a simple scalar equation with a distributed delay term. Example 13. The scalar equation dx(t ) dt

= ax(t ) + b



a < 0,

has a characteristic function g (s, h) = s − a − b

1 − e−hs s



Example 14. Consider the system dx(t )

= A0 x(t ) + µA1 x(t − 1) +

dt



0

[H0 + θ H1 ] x(t + θ )dθ ,

−1

where µ is a real parameter of the system and

 A0 =

−1

 H0 =



0

0 , −1

0.3 0

0 , 0.3

 A1 =



0 −1

 H1 =

0 −0.3



1 , 0 0.3 . 0



Its characteristic function is g (s, µ) = det sI − A0 − µA1 e−s − f (0) (s)H0 − f (1) (s)H1 ,





where 1 − e −s s

,

f (1) (s) =

e−s + he−s − 1 s2

.

q(s) = s8 + s6 −2µ2 − 3.2 + s4 µ4 + 3.2µ2 − 1.2µ + 2.47

,

and the characteristic polynomial of the delay free system (17) is p(s) = s(s2 − a2 − 2b). The substitution of the root s = 0 of the first factor of p(s) into g (s, h) leads to the condition

− a − bh = 0,

Now, a two dimensional system with distributed delay and an unknown parameter is considered.

The characteristic polynomial of the corresponding delay free system (17) is now p(s) = q(s)r (s) where

−h



(39)

Collecting the above, we see that for a < 0, the boundaries of the delay dependent stability region, I ∪ II on Fig. 2, are described by Eqs. (38) and (39).

f (0) (s) =

0

x(t + θ )dθ ,

  √  2 − 2b − a a h −a2 − 2b = arctan . −a 2 − b

(38)

which reduces to b = 0 when h → ∞. Next, observe that for a2 + 2b > 0 the second factor of p(s) has no pure imaginary roots.

      + s2 1.2µ3 − 0.18µ2 + 1.92µ − 0.288   + 0.36µ2 − 0.108µ + 0.1845     r (s) = s8 + s6 2µ2 − 3.2 + s4 µ4 − 3.2µ2 + 1.2µ + 2.38   + s2 1.2µ3 + 0.18µ2 − 1.2µ + 0.648   + 0.36µ2 − 0.108µ + 0.1845 . Let us set λ = s2 . The analysis of Sturm’s array of the polynomial q1 (λ) shows that for any real µ, it has non negative real roots. Similarly, Sturm’s array of the auxiliary polynomial r1 (λ) gives the

788

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

following root configuration:

µ

Number of negative real roots of q1 (λ) 2 0 2 4

[0, 0.3164) [0.3164, 1.2212) [1.2212, 2.3735) [2.3735, ∞)

For µ ∈ [0, 10], the roots are sketched in Fig. 3. Next, we must check if there exists a pair (iωk (µ), µ) , k = 1, 2, for which g (s, µ) vanishes. Direct calculations show that no pair of values (ω2 (µ), µ) and (ω4 (µ), µ), annihilates this characteristic function. Fig. 4 shows that the function g (iω1 (µ), µ) vanishes for the values (ω1 (µ) = 0.702, µ = 1.224) and (ω1 (µ) = 6.4293, µ = 6.5061). We observe in Fig. 5 that the real and imaginary parts of g (iω3 (µ), µ) vanish simultaneously, when (ω3 (µ) = 3.3945, µ = 3.7624) and (ω3 (µ) = 9.5252, µ = 9.6474). Fig. 3. Roots ω1 (µ) , ω2 (µ) , ω3 (µ) and ω4 (µ).

Next, we consider an integral system. Example 15. Let us now analyze the system x(t ) =

0



F (θ ) x (t + θ ) dθ ,

(40)

−h

where



F (θ ) = eAθ

A=

and



−1

0 . 1

2

According to Remark 7, computing the derivative of system (40) leads to the delay system with a distributed delay of the form: dx(t )

= A0 x(t ) + A1 (h)x (t − h) +

dt



0

H (θ ) x (t + θ) dθ , −h

where A0 = I , A1 = A sinh(h) − I cosh(h), and H (θ) = − [I sinh (θ) + A cosh (θ )]. The characteristic function is of the form

 g ( s , h) =

s − 2 + e−h(s−1)

s + e−h(s+1)



s2 − 1

 .

Fig. 4. Real and imaginary parts of g (iω1 (µ) , µ).

The characteristic polynomial of the corresponding delay free system (17) is

12 

p (s, h) = s8 s2 − 1



s2 + (cosh h + sinh h)2 − 4



  × s2 + (cosh h − sinh h)2 . Clearly the only term of the polynomial that does not have roots on

12



the imaginary axis is s2 − 1 . Now, we must check if there exists a pair (iω, h) that annihilate the characteristic function g (s, h). Direct computations lead us to determine that the only pair that satisfies g (s∗ , h∗ ) = 0 is s∗ = 0 and h∗ = 0.6931; see Fig. 6. Finally, we consider a neutral type time delay system with an unknown parameter. Example 16. Let us consider the following system dx(t ) dt

+ B1

dx (t − 1) dt

= µA0 x(t ) + A1 x (t − 1) ,

where

 A0 = B1 =

−2 0

 −0.1 0



0 , −0.9



0 , −0.1

 A1 =

−1 −1



0 , −1

Fig. 5. Real and imaginary parts of g (iω3 (µ) , µ).

and µ > 0 is not known. The spectral radius ρ (B1 ) = 0.1 < 1, thus the difference operator is stable. The characteristic polynomial of

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

789

0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1

0.2

0

0.4

0.6

0.8

1 Fig. 7. Real and imaginary parts of g iω1∗ , µ .



Fig. 6. Real and imaginary parts of g (s∗ , h).



the corresponding delay free system is p(s, µ) =

1 994

q (s, µ) r (s, µ) v (s, µ) ,

where q (s, µ) = 99s2 − 81µ2 + 100





r (s, µ) = 99s2 − 400µ2 + 100





v (s, µ) = (9801s4 + s2 (−47 740µ2 + 19 800) + 32 400µ4 − 36 000µ2 + 10 000). With the change of variable λ = s2 we can directly check that the auxiliary  10  polynomial q1 (λ, µ) has one negative real root for µ ∈ 0, 9 , while the auxiliary polynomial r1 (λ, µ) has one negative

root for µ ∈ 0, 12 . With the help of Sturm’s 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





g (s, µ) = det sI + sB1 e



−s

− µA0 − A1 e

−s



Fig. 8. Real and imaginary parts of g iω2∗ , µ .



.

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 Fig. 7, and for ω2∗ = 0.3866 and µ2 = 1.0256, see Fig. 8, the real and imaginary parts vanish simultaneously. 6. Concluding remarks A novel approach for determining the critical values of parameters of time delay systems for which the characteristic function of the system has pure imaginary roots is presented. It covers standard classes of linear delay systems with concentrated delays and distributed delays, as well as neutral type systems. The proposed methodology based on Sturm’s Theorem is validated by a selection of non trivial examples. We would like to emphasize the fact that, unlike the classical algebraic tools such as the Kronecker sum and resultant matrix, this approach establishes an interesting connection between the frequency domain stability analysis of time delay systems, and numerical schemes used for the computation of the Lyapunov matrices of the systems. Appendix. Sturm’s Theorem For the reader’s convenience, Sturm’s Theorem is recalled below.



Definition 17 ([19]). Let p(s) be a polynomial with real coefficients. Then Sturm’s sequence (P0 (s), P1 (s), . . . , Pn (s)) is defined as: P0 (s) = p(s) P1 (s) = p′ (s) Pn (s) = −rem (Pn−2 (s), Pn−1 (s))

n≥2

where rem (Pn−2 , Pn−1 ) denotes the reminder of the division of Pn−2 (s) and Pn−1 (s). The sequence stops when Pk (s) is a nonzero constant. Theorem 18 ([19]). Let p(s) be a polynomial with real coefficients and (P0 (s), P1 (s), . . . , Pn (s)) its Sturm’s sequence. Let a, b ∈ R, with a < b, such that p(a) ̸= 0 and p(b) ̸= 0. Then the number of distinct real roots of the polynomial p (s) in (a, b) is v(a)−v(b), where v(c ) denote the number of variations in the sign of Sturm’s sequence (P0 (s), P1 (s), . . . , Pn (s)) evaluated in s = c. References [1] R. Bellman, K.L. Cooke, Differential-Difference Equations, Academic Press, New York, 1963. [2] P. Lancaster, M. Tismenetsky, The Theory of Matrices, second ed., in: Comp. Sciences Appl. Math. Series, Academic Press, Orlando, 1985.

790

G. Ochoa et al. / Systems & Control Letters 62 (2013) 781–790

[3] J. Chen, G. Gu, C. Nett, A new method for computing delay margins for stability of linear delay systems, Systems Control Lett. 26 (2) (1995) 107–117. [4] A.F. Ergenc, N. Olgac, H. Fazelina, Extended kronecker summation for cluster treatment of LTI systems with multiple delays, SIAM J. Control Optim. 46 (1) (2007) 143–155. [5] E. Jarlebring, Critical delays and polynomial eigenvalue problems, J. Comput. Appl. Math. 224 (1) (2009) 296–306. [6] S.-I. Niculescu, Stability and hiperbolicity of linear systems with delayed states: a matrix pencil approach, IMA J. Math. Control Inform. 15 (4) (1998) 331–347. [7] K.E. Walton, J.E. Marshall, Direct method for TDS stability analysis, IEE Proc. 134 (D) (1987) 101–107. [8] N. Olgac, R. Sipahi, An exact method for the stability analysis of time-delayed linear time invariant (LTI) systems, IEEE Trans. Automat. Control 47 (5) (2002) 793–797. [9] R. Sipahi, N. Olgac, A comparative survey in determining the imaginary characteristic roots of LTI time delayed systems, in: Proc. of the 16th IFAC World Congress, 2005. [10] J. Louisell, A matrix method for determining the imaginary axis eigenvalues of a delay system, IEEE Trans. Automat. Control 46 (12) (2001) 2008–2012. [11] J. Louisell, Numerics of the stability exponent and eigenvalue abscissas of a matrix delay system, in: L. Dugard, E.I. Verriest (Eds.), in: Lecture Notes in Control and Information Sciences, vol. 228, Springer, 1998. [12] E.F. Infante, W.B. Castelan, A Lyapunov functional from a matrix differencedifferential equation, J. Differential Equations 29 (3) (1978) 439–451. [13] H. Garcia-Lozano, V.L. Kharitonov, Lyapunov matrices for time delay systems with commensurate delays, in: 2-nd IFAC Symposium on Systems, Structure and Control, Oaxaca, Mexico, 2004, pp. 102–106. [14] V.L. Kharitonov, Lyapunov matrices for a class of time delay systems, Systems Control Lett. 55 (7) (2006) 610–617.

[15] V.L. Kharitonov, S. Mondié, G. Ochoa, Frequency stability analysis of linear systems with general distributed delays, in: J.J. Loiseau, W. Michiels, S.I. Niculsecu (Eds.), in: Lecture Notes in Control and Information Sciences, vol. 388, Springer, 2009. [16] G. Ochoa, J.E. Velázquez Velázquez, V.L. Kharitonov, S. Mondié, Lyapunov matrices for neutral type time delay systems, in: J.J. Loiseau, W. Michiels, S.-I. Niculsecu (Eds.), in: Lecture Notes in Control and Information Sciences, vol. 388, Springer, 2012. [17] G. Ochoa, S. Mondié, V.L. Kharitonov, Computation of imaginary axis eigenvalues and critical parameters for neutral time delay systems, in: R. Sipahi, T. Vyhlídal, S.-I. Niculescu, P. Pepe (Eds.), in: Lecture Notes in Control and Information Sciences, vol. 423, Springer, 2009. [18] J.E. Marshall, H. Gorecki, A. Korytowski, K. Walton, Time-Delay Systems: Stability and Performance Criteria with Applications, Ellis Horwood, New York, 1992. [19] J.V. Uspensky, Theory of Equations, McGraw-Hill, New York, 1948. [20] J.K. Hale, S.M. Verduyn Lunel, Strong stabilization of neutral functional differential equations, IMA J. Math. Control Inform. 19 (1–2) (2002) 5–23. [21] V.L. Kharitonov, Lyapunov matrices for a class of neutral type time delay systems, in: 6th Workshop on Time Delay Systems, L’Aquila, Italy, 2006. [22] K. Gu, S.-I. Niculescu, Additional dynamics in transformed time-delay systems, IEEE Trans. Automat. Control 45 (3) (2000) 572–575. [23] V. Kharitonov, D. Melchor-Aguilar, Additional dynamics for general class of time-delay systems, IEEE Trans. Automat. Control 48 (6) (2003) 1060–1064. [24] C. Abdallah, P. Dorato, J. Benites-Read, R. Byrne, Delayed positive feedback can stabilize oscillatory systems, in: Proc. American Contr. Conf., 1993, pp. 3106–3107. [25] V.L. Kharitonov, S.-I. Niculescu, J. Moreno, W. Michiels, Static output feedback stabilization: necessary conditions for multiple delay controllers, IEEE Trans. Automat. Control 50 (1) (2005) 82–86.