Fazia Bedouhene
email: fazia.bedouhene@ummto.dz

Islam Boussaada
email: islam.boussaada@l2s.centralesupelec.fr

Silviu-Iulian Niculescu
email: silviu.niculescu@l2s.centralesupelec.fr

On Pole Placement and Spectral Abscissa Characterization for Time-delay Systems

Keywords: Time-delay systems, Stability, Spectral abscissa, Control design, Pole assignment, Non oscillation

This paper presents a systematic frequency domain approach to analyse the stability of reduced-order linear systems with single delay. More precisely, we address the problem of the spectral abscissa characterization and the coexistence of non oscillating modes for such functional differential equations. The dominancy of such non oscillating modes is analytically shown for the considered reduced order Time-delay systems.

INTRODUCTION

Investigation of dynamical systems with time-delay is an active research area that connects a wide range of scientific disciplines including mathematics, physics, engineering, biology, economics etc. The present paper focuses on stability and stabilizing-controllers design for linear time-invariant retarded time-delay systems. The study of conditions on the equation parameters that guarantees the exponential stability of solutions is a question of ongoing interest and remains an open problem especially when the systems are of high order or having multiple and/or distributed delays. In particular, in frequency-domain, the problem reduces to the analysis of the distribution of the roots of the corresponding characteristic equation, see for instance [START_REF] Bellman | Differential-difference equations[END_REF][START_REF] Cooke | On zeroes of some transcendental equations[END_REF][START_REF] Walton | Direct method for tds stability analysis[END_REF][START_REF] Stépán | Retarded Dynamical Systems: Stability and Characteristic Functions[END_REF][START_REF] Hale | Introduction to functional differential equations[END_REF][START_REF] Michiels | Stability and stabilization of time-delay systems[END_REF][START_REF] Sipahi | Stability and stabilization of systems with time delay[END_REF].

The starting point of the present work is an interesting property, discussed in recent studies, called Multiplicity-Induced-Dominancy denoted in the sequel (MID). As a matter of fact, it is shown that multiple spectral values for Time-delay systems can be characterized using a Birkhoff/Vandermonde-based approach; see for instance (Boussaada andNiculescu, 2016b,a, 2014;Boussaada et al., 2016). More precisely, in the previous works, it is emphasized that the admissible multiplicity of the zero spectral value is bounded by the generic Polya and Szegö bound denoted P S B , which is nothing but the degree of the corresponding quasipolynomial (i.e the number of the involved polynomials plus their degree minus one), see for instance [START_REF] Pólya | Problems and Theorems in Analysis[END_REF]. The multiplicity of a root itself is not important as such but its connection with the dominancy of this root is a meaningful tool for control synthesis. To the best of our knowledge, the first time an analytical proof of the dominancy of a spectral value for the scalar equation with a single delay was presented in [START_REF] Hayes | Roots of the transcendental equation associated with a certain difference-differential equation[END_REF]. The dominancy property is further explored and analytically shown in the case of second-order systems and a rightmost root assignment based design using delayed state-feedback is proposed in [START_REF] Boussaada | On the coalescence of spectral values and its effect on the stability of time-delay systems: Application to active vibration control[END_REF](Boussaada et al., , 2018b) ) where its applicability in damping active vibrations for a piezo-actuated beam is proved. See also (Boussaada et al., 2018a;Boussaada and Niculescu, 2018) which exhibit an analytical proof for the dominancy of the spectral value with maximal multiplicity for second-order systems controlled via a delayed proportional-derivative controller.

By this paper, we would like to extend such an analytical characterization of the spectral abscissa for retarded timedelay system with real spectral values which are not necessarily multiple. The effect of the coexistence of such non oscillatory modes on the asymptotic stability of the trivial solution will be explored. In particular, the coexistence of P S B real spectral values makes them rightmost-roots of the corresponding quasipolynomial. Furthermore, if they are negative, this guarantees the asymptotic stability of the trivial solution.

The remaining paper is organized as follows; in Section 2 we summarize some important facts on the coexistence non oscillatory modes for the first and second order time-delay differential equation with a single delay. Section 3 contains the main result of this paper. We address the problem of the spectral abscissa characterization and the coexistence of non oscillating modes for the third-order differential equation with a single delay. Section 4 is dedicated to illustrate the main result of the paper on Mach number regulation in a wind tunnel by the assignment of four dominant equidistributed real roots. In Section 5, some concluding remarks end the paper.

PRELIMINARIES

Let consider the generic n-order system with a single time delay: ẋ(t) = A 0 x(t) + A 1 x(tτ ).

(1) Here τ is a positive constant delay and the matrices A j ∈ M n (R) for j = 0 . . . 1. It is well known that the asymptotic behavior of the solutions of (1) is determined from the spectrum designating the set of the roots of the associated characteristic function (denoted in the sequel ∆(s, τ )). Namely, the characteristic function corresponding to system (1) is a quasipolynomial ∆ :

C × R + → C of the form: ∆(s, τ ) = det(sI -A 0 -A 1 e -τ s ).
(2) Asymptotic stability of the trivial solution and oscillatory behavior of (1) are known. In particular, the zero solution of this equation is asymptotically stable if and only if all roots of (3) lie in the left half plane, and a given solution of ( 1) is said to be non oscillatory if it corresponds to a real root of (3). This section summarizes the main findings reported in [START_REF] Amrane | On qualitative properties of low-degree quasipolynomials: Further remarks on the spectral abscissa and rightmost-roots assignment[END_REF], concerning the question of coexistence of P S B negative real roots for (3), when in particular ∆(s, τ ) = P 0 (s) + P 1 (s)e -τ s .

(3) where deg(P 0 ) = n, deg(P 1 ) = 0 for n = 1 and n = 2.

Let's consider the following systems

ẋ(t) + ax(t) + bx(t -τ ) = 0, (4) ẍ(t) + a ẋ(t) + bx(t) + αx(t -τ ) = 0. ( 5
)
The characteristic equation associated to (4) and ( 5) are respectively as follows:

∆ 1 (s, τ ) := s + a + b exp(-sτ ) = 0, (6) ∆ 2 (s, τ ) = s 2 + as + b + α exp(-sτ ) = 0. ( 7
)
Theorem 1. For a given delay τ > 0, the system (4) admits two distinct real spectral values at s = s 1 and s = s 2 , with s 2 < s 1 , if and only if

     a = a(s 1 , s 2 , τ ) := s 2 exp(-s 1 τ ) -s 1 exp(-s 2 τ ) exp(-s 2 τ ) -exp(-s 1 τ ) ; b = b(s 1 , s 2 , τ ) := s 1 -s 2 exp(-s 2 τ ) -exp(-s 1 τ ) . (8)
• Moreover, both spectral values s 1 and s 2 of (4) are negative, if and only if equation a(s 1 , s 2 , τ ) = 0 admits a positive solution in τ . Furthermore, in such a case the zero solution of ( 4) is asymptotically stable.

• The spectral value s 1 is nothing but the spectral abscissa corresponding to (4). Theorem 2. The system (5) admits three distinct real spectral values s 3 , s 2 and s 1 with s 3 < s 2 < s 1 if and only if the parameters a, b and α satisfy

                       a(τ ) := 1 Q(τ ) i,j,k∈Λ i<j,i =j =k (-1) i+j (s 2 i -s 2 j ) exp(-s k τ ) b(τ ) = - 1 Q(τ ) i,j,k∈Λ i<j,i =j =k (-1) i+j s i s j (s i -s j )e -s k τ α(τ ) = - 1 Q(τ ) i,j∈Λ i<j (s i -s j ) (9) where Q(τ ) = i,j,k∈Λ i<j,k =i,j (-1) i+j (s i -s j ) exp(-s k τ ).
In this case, α is necessarily negative.

• The spectral value s 1 is negative if and only if there exists τ 0 > 0 such that a(τ 0 ) + s 2 = 0. This guarantees the asymptotic stability of the system.

• The root s 1 is the spectral abscissa of (5).

Application to a control problem

Let us focus on the second order system ẍ(t) + a ẋ(t) + bx(t) = u(t), (10) where u is the unknown control and a and b are known parameters. Assume that the system (10) is unstable in the uncontrolled case, namely when u(t) = 0. This arises for instance if a < 0. Our aim is to design a control u under the form: u(t) = -αx(t)βx(tτ ), (11) that stabilizes the closed loop system: ẍ(t) + a ẋ(t) + (b + α)x(t) + βx(tτ ) = 0, (12) by pole assignment method. It is therefore a question to assign three negative spectral values, s 1 , s 2 and s 3 , to the characteristic equation associated to (12). This pole assignment is then interpreted by the conditions (9), in which case b := b + α, allowing the computation of the parameters α, β and τ > 0 of the control (11). By choosing s 1 , s 2 and s 3 equi-distributed (s 1s 2 = s 2s 3 = d > 0), we obtain, in view of (12), the following values of the parameters τ , α and β:

                   b = 1/2 a 2 + as 1 -ad -2 s 1 d + s 1 2 + 1/2 d 2 , β = -1/2 -8 s 1 d + 3 d 2 -4 ad + 4 s 1 2 + 4 as 1 + a 2 × -3 d + 2 s 1 + a 2 s 1 -d + a - -s 1 +d d , τ = ln -3 d + 2 s 1 + a 2 s 1 -d + a d -1
The asymptotic stability of the nontrivial solution of the closed-loop system (12) is guaranteed by the dominance property of s 1 established in Theorem 2, and the positivity of the delay τ .

2.2 About the geometric structure of the envelope curve associated to (6) and (7)

(1) It is well known that the classical envelope curve has a connected geometric structure (see for instance [START_REF] Silviu-Iulian | Delay Effects on Output Feedback Control of Dynamical Systems[END_REF]). Interestingly, when considering two real distinct spectral values, the connected structure of the envelope may be lost, see Fig. 1. (2) Such a geometry is encountered in [START_REF] Qiao | A linear time-invariant consensus dynamics with homogeneous delays: Analytical study and synthesis of rightmost eigenvalues[END_REF], where an analytical study and synthesis of rightmost eigenvalues of ẋ(t) = Ax(tτ ) is considered.

(3) Likewise, the geometric structure of the envelope curve of the quasipolynomial (7), defined by

x 2 + y 2 -A 0 2 -A 1 2 e -τ x = 0, with A 0 = -a 1 -b 0 and A 1 = 0 0 -α 0
, may loose its connection as observed in the first order equation, depending on the distribution of the roots s 1 , s 2 , s 3 and the delay τ . More precisely, three cases can be observed according to the distance of the root s 3 with respect to the centered circle, of radius R = A 0 2 . Case s 1 = -2, s 2 = -3, s 3 = -13. Case s 1 = -2, s 2 = -3, and s 3 = -R (right), and 0 > s 3 > -R (left).

ON POLE-PLACEMENT FOR THIRD-ORDER EQUATIONS WITH A SINGLE DELAY

Four real poles assignment is possible

Let consider the generic 3-order equation with a single time delay

x (3) (t)+a 2 x (2) (t)+a 1 x (1) (t)+a 0 x(t)+αx(t-τ ) = 0. ( 13
)
The characteristic equation associated to ( 13) is given by: ∆ 3 (s, τ ) = s 3 + a 2 s 2 + a 1 s + a 0 + α exp(-sτ ) = 0. ( 14) To simplify some formulas, let us introduce the notation [s 1 , s 2 ] t := ts 1 + (1t) s 2 . Theorem 3. i) System (13) admits four distinct real spectral values s 4 , s 3 , s 2 and s 1 with s 4 < s 3 < s 2 < s 1 if and only if

Q(τ ) := 4 k=1 (-1) k exp(-s k τ ) 4 i,j=1 i<j,k =i,j (s i -s j ) = 0.
In this case, the coefficients a i , i = 1, . . . , 3, and α are uniquely determined as a continuous function with respect to the delay τ > 0. The parameter variable α(τ ) is necessarily positive for every τ > 0, and satisfies s2,[s3,[s4,s1] λ ] θ ] t dλdθdt ii) The spectral value s 1 is negative if and only if there exists τ 0 > 0 such that a 2 (τ 0 ) + s 2 + s 3 = 0. This guarantees the asymptotic stability of the system.

α (τ ) = 1 Q (τ ) 4 i<j,i,j=1 (s i -s j ) = 1 τ 3 [0,1] 3 (1 -t) 2 (1 -θ) e -τ [

Sketch of the proof

First, let us find conditions on the coefficients a i , for i = 0, • • • , 2, and α ensuring the coexistence of P S B real spectral values (here P S B = 4) for the quasipolynomial (14), recall that P S B is Polya and Szegö bound. For, assume that (14) admits four real spectral values s 1 > s 2 > s 3 > s 4 . This means that s 3 i +a 2 s 2 i +a 1 s i +a 0 +α exp(-s i τ ) = 0, for all i = 1, • • • , 4.

(15) From which we deduce the following values of the coefficients a i , for i = 0, • • • , 2, and α:

a 2 (τ ) = 1 Q(τ ) 4 k=1 (-1) k+1 e -s k τ 4 i,j=1 i<j k =i,j (s i -s j ) 4 i=1 i =k s i ; a 1 (τ ) = 1 Q(τ ) 4 k=1 (-1) k e -s k τ 4 i,j=1 i<j k =i,j (s i -s j ) 4 i,j=1 i<j k =i,j s i s j ; a 0 (τ ) = 1 Q(τ ) 4 k=1 (-1) k+1 e -s k τ 4 i,j=1 i<j k =i,j (s i -s j ) 4 i,j=1 i<j k =i,j s i s j ; α(τ ) = 1 Q(τ ) 4 i,j=1 i<j (s i -s j );
Clearly, a i , for i = 0, • • • , 2, and α are well-defined, for any value of the delay τ , if and only if the following condition Q (τ ) = 0, ∀τ > 0. ( 16) is fulfilled.

So, let us check the condition (16) . By rearranging the terms in Q, and using the mean value theorem, we get

Q (τ ) = τ 3 4 i<j i,j=1 (s i -s j ) Q(τ ) > 0 (17) with Q(τ ) = [0,1] 3 (1 -t) 2 (1 -θ) e -τ [s2,[s3,[s4,s1] λ ] θ ] t dλdθdt.
The existence and uniqueness of the coefficients a 2 (τ ), a 1 (τ ), a 0 (τ ), α (τ )) is then proved. The positiveness of α is provided by ( 17) and the sign of Q.

Since the mapping τ → a 2 (τ ) + s 2 + s 3 is continuous and increasing from -∞ to -s 1 when τ varies in R + * , this means that the mapping τ → a 2 (τ ) + s 2 + s 3 takes positive values if and only if s 1 < 0. Also , if and only if, there exists (a unique) root τ 0 > 0 to equation a 2 (τ ) + s 2 + s 3 = 0.

The dominancy of s 1 for (13)

To study the stability of the system (13), we need to study the dominancy of s 1 by using an adequate factorization of the quasi-polynomial ∆ 3 in ( 14). Theorem 4. The root s 1 is the spectral abscissa of (13).

Sketch of the proof

Rewrite the quasipolynomial ∆ 3 as:

∆ 3 (s, τ ) = (s -s 1 ) (s -s 2 ) (s -s 3 ) P (s, τ ) with P (s, τ ) = s 3 + a 2 s 2 + a 1 s + a 0 + α exp (-τ s) (s -s 1 ) (s -s 2 ) (s -s 3 ) . Define the quantities: b 0 := a 0 + s 1 s 2 s 3 ; b 1 := a 1 - 3 i,j=1
i =j s i s j ; and b 2 := a 2 + 3 i=1 s i . Some tedious algebraic manipulations allows to write P (s, τ ) P (s, τ ) = 1τ 3 α P (s, τ ) where s,[s2,[s1,s3] λ ] θ ] t dλdθdt.

P (s, τ ) = [0,1] 3 (1 -t) 2 (1 -θ) e -τ [
To prove dominancy property for s 1 , let us assume that there exists some s 0 = ζ + jη a root of ( 14) such that ζ > s 1 . This means that P (s 0 , τ ) = 0. Combining this fact and the positiveness of α, we get s,[s1,[s2,s3]

1 = τ 3 α P (s, τ ) ≤ τ 3 α [0,1] 3 (1 -θ) e -τ [
λ ] θ ] t dλdθdt (18) < τ 3 α [0,1] 3 (1 -t) 2 (1 -θ) e -τ [s1,[[s2,s3] λ ,s4] θ ] t dλdθdt = τ 3 αQ (τ ) = 1 which is inconsistent.

MACH NUMBER REGULATION IN A WIND TUNNEL: EQUIDISTRIBUTED DOMINANT-ROOTS ASSIGNMENT

Roughly speaking, the Mach number regulation in a wind tunnel is based on the Navier-Stokes equations for unsteady flow and contains control laws for temperature and pressure regulation.

As an illustrative example for the applicative potential of the proposed main result, let revisit the following simplified model of Mach number regulation proposed in [START_REF] Manitius | Feedback controllers for a wind tunnel model involving a delay: Analytical design and numerical simulation[END_REF] and consists of a system of three state equations with a delay in one of the state variables. It is stressed that in steady-state operating conditions, the dynamic response of the Mach number perturbations ξ 1 to small perturbations in the guide vane angle actuator ξ 2 are governed by:

     ξ1 (t) = -aξ 1 (t) + k a ξ 2 (t -τ ) ξ2 (t) = ξ 3 (t) ξ3 (t) = -ω 2 ξ 2 (t) -2ζωξ 3 (t) + ω 2 u(t) (19)
where a, ω, ζ, k and τ are parameters depending on the operating point and presumed constant when the perturbations ξ i are small. Moreover, following the experimental parameter values of the wind tunnel developed at NASA Langley Research Center, the parameters a, ω, ζ, τ are positive.

In [START_REF] Manitius | Feedback controllers for a wind tunnel model involving a delay: Analytical design and numerical simulation[END_REF], a feedback consisting of a linear combination of state variables and weighted integrals of some of the state variables over a period equal to the time delay, where the spectrum of the closed-loop system is finite (consists of three eigenvalues). However, our method does not render the closed-loop system finite dimensional but only involves controlling its rightmost root. In (Boussaada et al., 2018a) the control law u(t) = -α ω 2 ξ 2 (t) -β0 ω 2 ξ 2 (tτ ) -β1 ω 2 ξ 3 (tτ ) is proposed allowing to the closed-loop quasipolynomial function: ∆(s, τ ) = (s+a) (sβ 1 +β 0 )e -sτ +s 2 +2 s ζ ω+ω 2 +α . (20) Thanks to such a factorization and since a is a positive parameter, the aim in (Boussaada et al., 2018a) were to establish conditions on parameters such that the rightmost root of the second factor of (20) has a negative real part and the MID property for second order systems were exploited. Here, we propose the control law: u(t) = αξ 1 (t) + βξ 2 (t) + γξ (22) as well as the precise value of the spectral abscissa and the distance between two successive assigned roots:

s 1 = - aτ + ln (2) τ , d = ln (2) τ .
Fig. 4. The spectrum distribution of ( 21) satisfying ( 22) exhibiting the four-roots placement in two configurations. In blue, the spectrum corresponding to a = 5, τ = 2. In red, the spectrum corresponding to a = 3, τ = 1. In both cases, the dominancy of the four assigned equidistributed real rootsis underlined.

Fig. 1 .

 1 Fig. 1. Envelope curve of the characteristic equation (6).Case of co-existence of two simple real roots s 1 = -1, s 2 = -2. Here τ * ≈ 0.67288 and τ * * = ln(2).

Fig. 2 .

 2 Fig. 2. Envelope curve of the characteristic equation (7).Case s 1 = -2, s 2 = -3, s 3 = -13.

yFig. 3 .

 3 Fig. 3. Envelope curve of the characteristic equation (7).Case s 1 = -2, s 2 = -3, and s 3 = -R (right), and 0 > s 3 > -R (left).

  3 (t) which gives in closed-loop the following quasipolynomial function ∆(s, τ ) =s 3 + (2 ζωω 2 γ + a)s 2 + (1β -aΓ)ω 2 + 2aζω s + (1β) aω 2ω 2 α e -sτ ka(21) Using Theorem 3 and assuming that s l = s 1 -(l -1) d for l = 2, . . . , 4 one obtains the following values of the controller gains:

CONCLUSION

In this note, we extended some recent results by the authors on pole-placement for Time-delay systems. This new result emphasizes a new delayed controller-design based on the trivial solution's decay rate assignment. The potential applicability of the approach is illustrated through the regulation of Mach number in a wind tunnel. Further insights on the applicability of the presented method in damping active vibrations can be found in (Tliba et al., 2019).


Islam Boussaada

Silviu-Iulian Niculescu

Tracking the Algebraic Multiplicity of Crossing Imaginary Roots for Generic Quasipolynomials: A Vandermonde-Based Approach

A standard way in analyzing dynamical systems consists in identifying and understanding the eigenvalues bifurcations when crossing the imaginary axis. Efficient methods for crossing imaginary roots identification exist. However, to the best of the authors knowledge, the multiplicity of such roots was not deeply investigated. In recent papers by the authors [1],

. This paper extends such a constructive approach in investigating the multiplicity of crossing imaginary roots j ω where ω = 0 and establishes a link with a new class of functional confluent Vandermonde matrices. As a byproduct of the proposed approach, a symbolic algorithm for computing the LU-factorization for such matrices is provided. This allows computing a bound which is sharper than Polya-Szegö generic bound arising from the principle argument.

I. INTRODUCTION

Consider the following class of time-delay systems:

ẋ = N k=0 A k x(t -τ k ) (1)
where x = (x 1 , . . . , x n ) ∈ R n denotes the state-vector, under appropriate initial conditions belonging to the Banach space of continuous functions C([-τ N , 0], R n ). Here τ j , j = 1 . . . N are strictly increasing positive constant delays, τ 0 = 0, τ = (τ 1 , . . . , τ N ) and the matrices A j ∈ M n (R) for j = 0 . . . N . It is well known that the asymptotic behavior of the solutions of ( 1) is determined from its spectra ℵ designating the set of the roots of the associated characteristic function denoted in the sequel ∆(z, τ ), often called quasipolynomial, that is a transcendental polynomial in the Laplace variable z in which appears exponential terms induced by delays, see, [START_REF] Diekmann | Applied Mathematical Sciences, Functional, complex, and nonlinear analysis[END_REF] and [START_REF] Bellman | Differential-difference equations[END_REF]. The study of zeros of such a class of entire function [START_REF] Levin | Distribution of zeros of entire functions[END_REF] plays a crucial role especially in the analysis of asymptotic stability of the zero solution associated with System [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF]. Indeed, the zero solution is asymptotically stable if all the spectral values of [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF] are in the open left-half complex plane [START_REF] Michiels | Stability and stabilization of time-delay systems[END_REF]. Accordingly to this observation, the parameter space which is spanned by the coefficients of the polynomials P i , can be split into stability and instability domains (Nothing else that the socalled D-decomposition, see for instance [START_REF] Michiels | Stability and stabilization of time-delay systems[END_REF] and references therein). These two domains are separated by a boundary corresponding to the spectra consisting in roots with zero real parts called in the sequel Crossing Imaginary Roots (CIR). Identification of CIR for Time-delay systems may be done for instance by Rekasius Substitution [START_REF] Olgac | An exact method for the stability analysis of timedelayed lti systems[END_REF] or Matrix Pencil [START_REF] Michiels | Stability and stabilization of time-delay systems[END_REF] to mention only few efficient approaches. However, to the best of the authors knowledge, the multiplicity of such roots was not deeply investigated, even, the computation of the upper bound of the algebraic multiplicity of CIR still unsolved. The characterization of CIR i.e. identification of CIR as well as their associated multiplicities are useful to overcome several issues related to the qualitative description of solution of a given dynamical system with or without delays. Furthermore, the delays for which the characteristic equation has CIR are referred to as critical delays and the sensitivity of a given CIR, i.e., the derivative with respect to the delay, produce a complete stability picture by keeping track of the number of roots entering and leaving the right half-plane, see for instance [START_REF] Michiels | Stability and stabilization of time-delay systems[END_REF] and [START_REF] Jarlebring | Invariance properties in the root sensitivity of time-delay systems with double imaginary roots[END_REF]. Moreover, characterizing CIR has a twofold interest: in linear analysis as well as in nonlinear analysis perspectives. From a linear viewpoint, for instance, a complete characterization of CIR allows to estimate the number of the spectral values with positive real part (elements of ℵ + = {z ∈ C, ∆(z, τ ) = 0, (z) > 0}) for a given Time-delay system [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF]. Indeed, the main theorem from [START_REF] Hassard | Counting roots of the characteristic equation for linear delay-differential systems[END_REF] emphasizes the link between (ℵ + ) and ℵ 0 (the set of all CIR), both counted by multiplicity. Taking into account the quoted result and its potential consequences on the design of a new approach on the linear stability analysis of timedelay systems, a greater importance to the study of CIR become obvious. Such an issue is also emphasized from a nonlinear analysis viewpoint, which gives another motivation for the present investigation. When the unstable spectra is an empty set or equivalently the spectra is

ℵ = ℵ -∪ ℵ 0 (elements of ℵ -= {z ∈ C, ∆(z, τ ) = 0, (z) < 0}
), a complete knowledge of CIR as well as their multiplicities becomes crucial predominatingly when the center manifold and the normal forms theory are involved for giving an accurate local qualitative description of the studied dynamical system.

A stimulating idea from [START_REF] Niculescu | Stabilizing a chain of integrators using multiple delays[END_REF], which consists in exploiting the properties of Vandermonde matrices to control one chain of integrators by delay blocks, suggest to the authors in [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF] and [START_REF] Boussaada | Computing the codimension of the singularity at the origin for time-delay systems in the regular case: A vandermonde-based approach[END_REF] to extend the investigation of such structured matrices. A complete characterization of the root at the origin can be found in [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF], [START_REF] Boussaada | Computing the codimension of the singularity at the origin for time-delay systems in the regular case: A vandermonde-based approach[END_REF]. It is shown that the multiplicity of the the root at the origin for a generic system is bounded by the degree of the associated quasipolynomial, a bound which was already depicted in [START_REF] Polya | Problems and Theorems in Analysis[END_REF]. 1 Furthermore, in the sparse case, a sharper bound is established in [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF] showing that the optimal multiplicity is related to the number of the non vanishing parameters rather than the degree of the quasipolynomial. The present work is a natural continuation of the last cited papers, its aim is three fold. First, it provides a symbolic algorithm for the LU-decomposition for a class of functional confluent Vandermonde matrices with trigonometric inputs (section III). Secondly, it extends the Vandermonde constructive approach established in [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF] to characterize CIR on beyond the zero frequency (sections IV and V). Finally, in section V the bound of the algebraic multiplicity for a given CIR is provided.

II. PREREQUISITES

A. Problem formulation and used notations

A given Time-delay system (1) is characterized by a quasipolynomial function ∆ : C × R N + → C of the form:

∆(z, τ ) = det z I -A 0 - N k=1 A k e -τ k z (2)
which gives,

∆(z, τ ) = P 0 (z) + M k ∈S N,n P M k (z) e σ M k z = P 0 (z) + ÑN,n k=1 P M k (z) e σ k z (3)
where σ M k = -M k τ T , τ = (τ 1 , . . . , τ N ) is the delays vector, S N,n is the set of all the possible row vectors

M k = (M k 1 , . . . , M k N ) belonging to R N * such that 1 ≤ M k 1 +. . .+ M k
N ≤ n and ÑN,n = (S N,n ). For instance, S 3,2 = {(1, 0, 0), (0, 1, 0), (0, 0, 1), (2, 0, 0), (1, 1, 0), (1, 0, 1), (0, 2, 0), (0, 1, 1), (0, 0, 2)} , is ordered first by increasing sums ( N i=1 M k i ) then by lexicographical order. For instance, for S 3,2 , one has M 2 = (0, 1, 0), Ñ3,2 = 9 and the cross-talk delay -σ 5 = τ 1 + τ 2 , see [START_REF] Sipahi | Complete stability robustness of third-order {LTI} multiple time-delay systems[END_REF]. It is also assumed that a i,j stands for the coefficient of the monomial z j in P i , a 0 = (a 0,0 , . . . , a 0,n ) is the vector composed from the coefficients of the polynomial P 0 and p = a 1,0 , . . . , a 1,n1-1 , . . . , a Ñ ,0 , . . . , a Ñ ,n Ñ -1 .

Without any lost of generality and for the sake of simplicity, the quasipolynomial function can be written

∆(z, τ ) = P 0 (z) + Ñ k=1 P k (z) e σ k z , (4)
where all the σ k can be assumed to be distinct and they can be considered as auxiliary delays for the quasipolynomial. Indeed, when for some value of the delay vector τ there exists some k = k such that M k τ T = M k τ T then the number of the polynomials is reduced by considering a new family of polynomials P such that PM k = P M k + P M k . It follows from ( 2) that P 0 is a monic polynomial of degree n in z and the polynomials P M k are such that deg(

P M k ) = n - N s=1 M k s ≤ (n -1), ∀M k ∈ S N,n .
A given complex number z 0 is said to be a spectral value for (1) if ∆(z 0 , τ * ) = 0 for some delay vector τ * . Let denote by ∂ k z ∆(z, τ ) the k-th derivative of ∆(z, τ ) given by ( 3) with respect to the variable z. A spectral value z 0 is of

algebraic multiplicity m ≥ 1 if ∆(z 0 , τ * ) = ∂ z ∆(z 0 , τ * ) = . . . = ∂ m-1 z ∆(z 0 , τ * ) = 0 and ∂ m z ∆(z 0 , τ * ) = 0.
A generic bound for the multiplicity of a given spectral value is deduced from a complex analysis result by Pólya and Szegö in [13, p. 144]. Their proof is mainly based on Rouché Theorem and the principal argument. It gives a bound of the number of quasipolynomial roots in any horizontal strip β ≤ Im(z) ≤ α. Setting α = β = ω gives the number of spectral values z with Im(z) = ω, which is a natural bound for the multiplicity of a given CIR z = jω, denoted in the sequel P S . See also [START_REF] Wielonsky | A Rolle's theorem for real exponential polynomials in the complex domain[END_REF] for a modern formulation of the mentioned result.

B. Integrator+two delay blocks:

For the general linear scalar time-delay systems of arbitrary order with two delays, the geometric criteria from [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF] allows to identify all CIR. However, the multiplicity of such roots was not explored.

For each pair of positive and distinct delays τ 1 , τ 2 the equation ẋ(t) + a 1 x(t -τ 1 ) + a 2 x(t -τ 2 ) = 0 admits a triple spectral value at

z 0 = -(τ 1 + τ 2 )/(τ 1 τ 2 ) if and only if a 1 = τ 2 e -τ1z0 / τ 2 τ 1 -τ 2 1 , a 2 = τ 1 e -τ2z0 / τ 2 τ 1 -τ 2 2
. Moreover, under this configuration, the zero solution is asymptotically stable. However, the optimal multiplicity of a CIR is two and it is reached if and only if ω 0 = ±tan (α 0 )/τ 1 and τ 2 = -τ 1 α 0 /tan (α 0 )

a 1 = ω 0 cos (ω 0 τ 2 ) sin (ω 0 (τ 1 -τ 2 )) , a 2 = ω 0 cos (ω 0 τ 1 ) sin (ω 0 (τ 2 -τ 1 )) , (5)
where α 0 is a fixed point of α → tan • tan(α) or equivalently, the root of α = tan • tan(α) which is a nonempty set. Now, let us show the above claim. We start by the configuration allowing to a spectral value z 0 with multiplicity three, which is characterized by the solutions of the system of equations

∆(z 0 , τ ) = . . . = ∂ k z,...,z ∆(z 0 , τ ) = 0 for k ∈ {0, 1, 2} and ∂ 3 z,...,z ∆(z 0 , τ ) = 0
where the characteristic function is ∆(z) = z + a 1 e -zτ1 + a 2 e -zτ2 , which can be also written

∆(z) = (z -z 0 )   1 + τ2 e -(z-z 0 ) τ 1 τ1(τ2-τ1) + τ1 e -(z-z 0 ) τ 2 τ2(τ1-τ2) + z 0 (z -z 0 )   .
To prove the asymptotic stability of the zero solution, let assume that there exists z 1 = ζ + j η a root of ∆(z) = 0 such that ζ ≥ 0. Thus, z 1 is a root of the second factor of ∆, which gives: Let us check now that the upper bound of the multiplicity for CIR is two. Such a multiplicity is reached for jω 0 if and only if ∆(jω 0 , τ ) = ∂ z ∆(jω 0 , τ ) = 0. Let us proceed by elimination. The first equation ∆(jω 0 , τ ) = 0 gives, a 1 cos (ω 0 τ 1 ) + a 2 cos (ω 0 τ 2 ) + (ω 0 -a 1 sin (ω 0 τ 1 ) -a 2 sin (ω 0 τ 2 )) j = 0. The vanishing of the real part allows to a 2 = -a 1 cos (ω 0 τ 1 )/cos (ω 0 τ 2 ), which is substituted in the imaginary part. After equating the resulting expression one has a 1 = ω cos (ω τ 2 )/sin (ω τ 1 -ω τ 2 ). The obtained value for a 1 is then itself substituted in a 2 . The final values of a 1 and a 2 are then substituted in the equation ∂ z ∆(jω, τ ) = 0, which provides two equations (real and imaginary terms) in three left free parameters ω, τ 1 and τ 2 . Allowing as expected, to the expression of the crossing frequency ω 0 as well as τ * 2 (τ 1 still being free). Remark. For a given quasipolynomial function, the multiplicity of real roots may reach Polya-Szegö bound, see [START_REF] Boussaada | Computing the codimension of the singularity at the origin for time-delay systems in the regular case: A vandermonde-based approach[END_REF], [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF]. However, the multiplicity of an imaginary crossing root may exceed the dimension of the delay-free system but never reach Polya-Szegö bound. A generalization of the above observation and further analysis illustrate the main result.

1 = 2 i=1 c i e -(z1-z0) τi -1 τ i (z 1 -z 0 ) , where c i = 2 l=1 l =i τ l (τ l -τ i )

C. Rekasius Substitution for a second order system

It is well known that Rekasius substitution is a way for identifying the crossing roots for a given quasipolynomial, [START_REF] Sipahi | Complete stability robustness of third-order {LTI} multiple time-delay systems[END_REF]. It consists in offsetting the exponential type transcendental terms in a given quasypolynomial function by the following rational function: e -zτ = (1 -T z)/(1 + T z), where T is some positive real number. This transformation is accompanied with the relation:

τ = 2 ωc (arctan(ω c T ) + kπ) , k = 0, 1, • • • where ω c ∈ R * +
and ω c denotes the crossing frequency. Il is worthy to note that such a relation associate to each solution of the transformed characteristic function in the auxiliary variable T to a 2π/ω c -periodic countably infinite set of solutions of the original characteristic quasipolynomial function in the variable τ , see for instance [START_REF] Olgac | An exact method for the stability analysis of timedelayed lti systems[END_REF], [START_REF] Munz | Stability analysis of time-delay systems with incommensurate delays using positive polynomials, Automatic Control[END_REF]. However, the question on roots multiplicity-conservation by such a transformation may be more deeply investigated. For instance, by the following example, we show that even for the delay margin2 τ 0 , such a multiplicity is not preserved by Rekasius substitution. Let consider the quasipolynomial function :

f (z, τ ) = a 0,0 + a 0,1 z + z 2 + (a 1,0 + a 1,1 z) e -τ z (6) with a 0,0 = π 2 +4 π 2 -4 , a 0,1 = a 1,0 = -4π π 2 -4 , a 1,1 = -8 π 2 -4 .
First, one knows that the imaginary roots of f are the ones of the following rational function in the variables z and T :

f (z, T ) = P(z, T ) Q(z, T )
where Q(z, T ) = 1 + T 1 z and

P(z, T ) = a 0,0 z + a 0,1 z 2 + z 3 -a 1,0 z -a 1,1 z 2 T + a 0,0 + a 0,1 z + z 2 + a 1,0 + a 1,1 z.
Elementary algebraic manipulations yields to z 0 = j is a root of f when T * = 1. Thus, z 0 = j is a root of f where the associated delays τ 0 = π/2 is easily deduced. Finally, we check that when τ = τ 0 , one has z 0 = j is a double root for f . However, the evaluation of ∂ z f at (z, T ) = (z 0 , T * ) gives: ∂ z f (z 0 , T * ) = 2 (πj -2)/(π + 2). Thus, z 0 = j is a simple root for f when T = T * , leading us to the following: Remark. For a given quasipolynomial function, Rekasius substitution yields to identify all roots on the imaginary axis but does not necessarily preserve their multiplicity. Furthermore, the other limitations of such a substitution arise when at least one of the polynomials (involved in the considered quasipolynomial) depend on z and τ . Which is obviously the case when exploring the multiplicity of CIR.

III. FUNCTIONAL CONFLUENT VANDERMONDE

MATRICES

By functional confluent Vandermonde matrix Φ we designate the matrix defined by:

     Φ =[Φ 1 Φ 2 . . . Φ M ] ∈ M l,δ (R), Φ i =[f (σ i ) f (1) (σ i ) . . . f (di-1) (σ i )], f (σ i ) =g(σ i ).[1 . . . σ l-1 i ] T , for 1 ≤ i ≤ M, (7)
for sufficiently regular function g ∈ C k (R). In the sequel, we are particularly interested by square matrices then it is assumed that M i=1 d i = δ = l. When g(x) = 1 then we are dealing with the so-called confluent Vandermonde matrix. When additionally d i = 1 for i = 1 . . . N then we recover the classical Vandermonde matrix and in this case M = δ since Φ is assumed to be a square matrix. An explicit recursive formula for the LU-factorization for a class of functional confluent Vandermonde matrices characterized by g(x) = x k is available in [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF] or in its expended version [START_REF] Boussaada | On codimension of zero singularities for time-delay systems: A link with vandermonde and birkhoff incidence matrices[END_REF]. The following notations are adopted. Let ξ stands for the vector composed from σ i counting their repetition d i through columns of a given matrix Φ, i.e.

ξ = (σ 1 , . . . , σ 1 d1 , . . . , σ M , . . . , σ M d M ).
It is important to note that the vector ξ, which will be called in the sequel the incidence vector, is a compact way to indicate the so-called incidence matrix associated with Φ, see for instance [START_REF] Gonzalez-Vega | Applying quantifier elimination to the Birkhoff interpolation problem[END_REF]. It characterizes a given square functional confluent matrix Φ defined in [START_REF] Levin | Distribution of zeros of entire functions[END_REF] by indicating its variables as well as their repetitions which indicate the associated orders of derivations through its columns.

In the light of the above notations and under the setting d 0 = 0, without any loss of generality: k) , where 0 ≤ r ≤ ν -1 and α ≤ d r+1 , here ν(k) denotes the index of component of σ associated with ξ k , that is σ ν(k) = ξ k and by κ(k) the order of ξ k in the sequence of ξ composed only by σ ν(k) . Obviously, ν(k) = r + 1 and κ(k) = α. Proposition 1. The functional confluent Vandermonde square matrix Φ given in (7) admits the following LUfactorization with unitary diagonal elements L i,i = 1:

ξ k = ξ d0+...+dr+α = ξ ν(k)-1 l=0 d l +κ(
               L i,1 =σ i-1 1 for 1 ≤ i ≤ δ, U 1,j =Φ 1,j for 1 ≤ j ≤ δ, L i,j =L i-1,j-1 + L i-1,j ξ j for 2 ≤ j ≤ i, U i,j =(κ(j) -1) U i-1,j-1 + U i-1,j σ ν(j) -ξ i-1 for 2 ≤ i ≤ j. ( 8
)
Proof: Here, we follow the same steps of the proof given in [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF], namely, a 2-D recurrence, see also the extended version [START_REF] Boussaada | On codimension of zero singularities for time-delay systems: A link with vandermonde and birkhoff incidence matrices[END_REF].

The determinant of the functional confluent Vandermonde matrix Φ follows directly from ( 8):

Corollary 1. The determinant of the functional confluent Vandermonde matrix Φ is given by det(Φ) = j=1 (U j,j ) where U j,j for 1 ≤ j ≤ δ are defined by:

                     U 1,1 = g(σ 1 ), U j,j = g(σ k+1 ) k l=1 (σ k+1 -σ l ) d l when j = 1 + d k for 1 ≤ k ≤ Ñ -1, U j,j = (j -1 -d k ) U j-1,j-1 when d k + 1 < j ≤ d k+1 for 1 ≤ k ≤ Ñ -1,

IV. CHARACTERIZING THE MULTIPLICITY OF CIR:

COMPUTATIONAL ISSUES Let denote by c k = cos(σ k ω) and s k = sin(σ k ω) and for a given real positif number x, x designate the integer part of x or equivalently the floor function at x.

Lemma 1. An imaginary complex number z = j ω is a root of ∂ k z k ∆(z, τ ) = 0 for k ≥ 0 if and only if the coefficients of P i for 0 ≤ i ≤ Ñ satisfy α k + αk = 0 and β k + βk = 0 (9)
where

                   α 0 = n i=0, i even (-1) i 2 ω i a 0,i and β 0 = n i=1, i odd (-1) i 2 ω i a 0,i , α k = ∂ ω β k-1 and β k = ∂ ω α k-1 ∀k ≥ 1, αk = Ñ i=1 di l=0 a i,l ∂ l c i σ i k ∂σ l i and βk = Ñ i=1 di l=0 a i,l ∂ l s i σ i k ∂σ l i .
Let denote by γ k = (α k + αk ) 2 + (β k + βk ) 2 and the variety As a direct consequence of the above Lemma 1, one has: Proposition 2. An imaginary crossing root z = j ω for Eq (3) is of multiplicity m if and only if the variety V m is non empty and γ m+1 = 0, where

V m = (ω, a 0,0 , . . . , a Ñ ,d Ñ , σ 1 , . . . , σ Ñ ) ∈ R n+δ+1 × R Ñ + , γ k = 0 for k = 0, . . . , m -1}
Remark. The above proposition provide a system of trigonometric-polynomials to solve. This allows to characterize multiple CIR. A standard approche in solving such systems of nonlinear trigonometric polynomials suggest to consider the sin(ωσ i ) and cos(ωσ i ) functions as new variables s i and c i respectively and then intentionally add some equations expressing s 2 i + c 2 i = 1. Thus, Groebner basis approches may be involved to characterize the solutions of the problem, [START_REF] Cox | Ideals, varieties, and algorithms. An introduction to computational algebraic geometry and commutative algebra[END_REF]. Unfortunately, in our case, even though the added algebraic equations the links between the variables s i and c i are provided, but, the link with the variable ω (which is involved in the argument of cos and sin) may be lost. This often introduce inconsistant solutions.

In what follow, we provide a linear algebra framework allowing to overcome the emphasized computational issue.

V. LINKS WITH FUNCTIONAL VANDERMONDE MATRICES

In the light of Lemma 1 as well as the auxiliary results established in section III, we are now able to prove the following proposition. For the sake of simplicity, for a given positive integer m let α m , β m stand for the vectors: as well as the respective functions g A (x) = cos(ωx) and g B (x) = sin(ωx). Proposition 3. An imaginary crossing root z = j ω for Eq (3) is of multiplicity m if and only if its frequency ω satisfy the following linear system:

α m (a 0 , ω) = -(α 0 , . . . , α m-1 ) , β m (a 0 , ω) = -(β 0 , . . . , β m-1 ) . (10) Let =
A m (ω, σ).p = α m (a 0 , ω), B m (ω, σ).p = β m (a 0 , ω). (11)
Proof: The proposition follows directly from Lemma 1 by considering the two ideals consisting (respectively) in the real part α k + αk (respectively the imaginary part β k + βk ) of the successive derivatives of the quasipolynomial function. Which gives rise to the system [START_REF] Hassard | Counting roots of the characteristic equation for linear delay-differential systems[END_REF].

Corollary 2. When the square matrices A and B are non degenerate then the multiplicity of any imaginary crossing root z = j ω 0 for Eq (3) is bounded by . Corollary 3. The functional confluent Vandermonde matrix A (respectively B ) is invertible if and only if ∀ 1 ≤ i = j ≤ Ñ we have ωσ i = π 2 + kπ and σ i = σ j (respectively ωσ i = kπ and σ i = σ j ).

Remark. A necessary and sufficient condition for A and B to be invertible is

σ i = σ j ∀ 1 ≤ i = j ≤ Ñ and ωσ i = kπ 2 , k ∈ N.
Moreover, a given CIR has the multiplicity if the both of the linear systems given by [START_REF] Hassard | Counting roots of the characteristic equation for linear delay-differential systems[END_REF] admit the same real solution. Thus, the above multiplicity bound for a given CIR is sharper than P S Polya-Szegö bound. Indeed,

P S = n + Ñ i=1 (d i + 1).

VI. ILLUSTRATIVE EXAMPLES

A. Scalar equation with two delays

Let us consider the scalar delay equation: ẋ(t) + a 1,0 x(tτ 1 ) + a 0,1 x(t -τ 2 ) + a 0,0 x(t) = 0. Here p = (a 1,0 , a 0,1 ) and Ñ = 2. Following Corollary 1 the optimal multiplicity of a given CIR is bounded by = 2 < P S = 3. Moreover, a necessary condition for a CIR to reach multiplicity two is τ 1 = τ 2 and ωτ i = kπ 2 . Indeed, in such a configuration, we have σ 1,0 = -τ 1 and σ 0,1 = -τ 2 and the functional Vandermonde matrices:

A 2 = cos(σ 1,0 ω) cos(σ 0,1 ω) σ 1,0 cos(σ 1,0 ω) σ 0,1 cos(σ 0,1 ω) , B 2 = sin(σ 1,0 ω) sin(σ 0,1 ω) σ 1,0 sin(σ 1,0 ω) σ 0,1 sin(σ 0,1 ω) , α 2 = -a 0,0 -1 and β 2 = -ω 0 .
Using Proposition 1, one proceeds by LU-factorisation for each of the matrices given in System [START_REF] Hassard | Counting roots of the characteristic equation for linear delay-differential systems[END_REF] and obtains a necessary and sufficient condition for a given CIR to reach its optimal multiplicity, that is to be double:

(a 0,0 σ 0,1 -1)sin (ω σ 1,0 ) = ω σ 0,1 cos (ω σ 1,0 ) (a 0,0 σ 1,0 -1)sin (ω σ 0,1 ) = ω σ 1,0 cos (ω σ 0,1 )
When a 0,0 = 0 then one recovers example II-B.

B. Planar system with two delays

One associates to the generic planar system with two delays the characteristic quasipolynomial function ∆(z, τ ) = P 0 (z) + P 1,0 (z)e -zτ1 + P 0,1 (z)e -zτ2 + P 2,0 (z)e -2zτ1 + P 0,2 (z)e -2zτ2 + P 1,1 (z)e -z(τ1+τ2) . Furthermore, following Corollary 3 the associated functional Vandermonde matrices A 7 (ω, σ) and B 7 (ω, σ) are non singular when σ i = σ j ∀ 1 ≤ i = j ≤ 7 and ωσ i = kπ 2 , k ∈ N. Under the assumption of Proposition 2, a given CIR has at most the multiplicity 7. However, the Polya-Szegö bound is P S = 9.

VII. CONCLUSION

In this paper we addressed the question on the characterization of multiple CIR for time-delay systems as well as their admissible and optimal multiplicities. The considered configuration is generic, allowing to the link with some Hermite interpolation problem through functional confluent Vandermonde matrices with trigonometric inputs. A general formulas for LU-factorization of such matrices is established. To the best of our knowledge, the use of this class of matrices as well as their factorisation could not be found in standard linear algebra textbooks. Finally, we emphasize that our approach may be immediately extended to characterize CIR for neutral systems. ACKNOWLEDGMENT APPENDIX Proof of Lemma 1 We emphasize that the coefficient α k (respectively β k ) is the real part (respectively the imaginary part ) of the polynomial P (k) 0 (jω). Furthermore, α k and β k are multivariate polynomials in ω and in the components of a 0 (the coefficients of P 0 in z). It is also important to note that, α s+1 = β s+1 = 0 for all integer s ≥ n. Consider the quasipolynomial given in (4). A complex number jω is a root of ∆(z, τ ) = 0 if and only if P 0 (jω) + Ñ k=1 P k (jω) e jσ k ω = 0 which gives: a 0,0 + . . . + a 0,n (jω) n + a 1,0 + . . . + a 1,d1 (jω) d1 e jσ1ω + . . . + a Ñ ,0 + . . . + a Ñ ,d Ñ (jω) d Ñ e jσ Ñ ω = 0.

Using Euler formula, one obtains: α 0 + jβ 0 + a 1,0 + . . . + a 1,d1 (jω) d1 (c 1 + js 1 ) + . . . + a Ñ ,0 + . . . + a Ñ ,d Ñ (jω) d Ñ (c Ñ + js Ñ ) = 0. Now, extracting the real and imaginary parts of the above expression allows us to prove [START_REF] Olgac | An exact method for the stability analysis of timedelayed lti systems[END_REF] for k = 0:

                               α 0 + c 1 a 1,0 -ωs 1 a 1,1 -ω 2

Fig. 1 .

 1 Fig. 1. Bifurcation of the double imaginary crossing root defined by (5) up to the standard normalisation τ * 2 = 1, the crossing frequency ω 0 = 11.53849769 and the delay τ * 1 = 0.1436275981 is subject to a positive perturbation ( = 0 . . . 6). The evolution of the spectral value from the nominal value = 0 (yellow) to the final value = 6 (red).

Fig. 2 .

 2 Fig. 2. Bifurcation of the double imaginary crossing root of (6) where the delay τ is subject to a positive perturbation ( = 0 . . . 6). The evolution of the spectral value from the nominal value = 0 (yellow) to the final value = 6 (red).

Ñ

  i=1 (d i +1), we define the functional confluent matrices A m and B m belonging to M m, (R) and characterized by the incidence vector ξ = (σ 1 , . . . , σ 1 d1 , . . . , σ Ñ , . . . , σ Ñ d Ñ )

c 1 a 1 , 2 + 1 a 1 1 a 1

 121111 . . . + ω d1 ∂ d1 c 1 ∂σ d1 ,d1 + . . . + c Ñ a Ñ ,0 + . . . + ω d Ñ ∂ d Ñ c Ñ ∂σ d Ñ Ñ = 0, β 0 + s 1 a 1,0 + ωc 1 a 1,1 -ω 2 s 1 a 1,2 + . . . + ω d1 ∂ d1 s 1 ∂σ d1 ,d1 + . . . + s Ñ a Ñ ,0 + . . . + ω d Ñ ∂ d Ñ s Ñ ∂σ d Ñ Ñ = 0.Consider the first partial derivatives of the quasipolynomial:∂ z ∆(z, τ ) = P 0 (z) + Ñ i=1 (P i (z) + σ i P i (z)) e σiz , ∂ 2 z 2 ∆(z, τ ) = P 0 (z) + Ñ i=1 (P i (z) + 2σ i P i (z) + σ 2 i P i (z))e σiz .By induction we arrive to :∂ k z k ∆(z, τ ) = P that z = jω is a root of ∂ k z k ∆(z, τ ) = 0 for k ≥ 1 and let us show that (9) is satisfied. (c i + js i ) . Thus, α k = -Re( (c i + js i )). Obviously, when k ≥ n then P (k) 0 (jω) = α k + jβ k = 0.The right hand side of the above equality is nothing else than αk . Indeed, l z l (m) z=jω (c i + js i ) z=jω (c i + js i ) = αk since one shows by induction w.r.t to l (for a fixed k) that ∂ l c i σ z=jω (c i + js i ) .

The degree of a given quasipolynomial is defined by the sum of the number of its polynomials plus their respective degrees minus one.

Rekasius substitution asserts that when an imaginary root jω 0 is associated to a given delay τ 0 then it is also associated to each value of the delay among the family {τ k } k∈Z = {τ 0 + k 2π |ω 0 | } k∈Z . But, as it was emphasized in[START_REF] Jarlebring | Invariance properties in the root sensitivity of time-delay systems with double imaginary roots[END_REF], if j ω 0 is double crossing root associated with the delay τ = τ 0 then, j ω 0 is a simple crossing root for any τ k .


Chi Jin
email: chi.jin@l2s.centralesupelec.fr

Silviu-Ilian Niculescu
email: silviu.niculescu@lss.supelec.fr

Islam Boussada
email: islam.boussaada@l2s.centralesupelec.fr

Keqin Gu

Stability Analysis of Control Systems subject to Delay-Difference Feedback

Keywords: Time-Delay Systems, Stability Analysis, Delay-dependent Coefficients, τ-decomposition Approach

In engineering practice delay-difference is often used to approximate the derivatives of output signals for feedback control, which leads to a closed-loop system with delay both in the states and in the system coefficients. Our objective is to find all the delay values below some upper bound that guarantee the exponential stability of the closed-loop system subject to the delay-difference approximation. A method for stability analysis of systems with delay-dependent coefficients developed in our previous work is extended and applied for the systems considered in this paper. The proposed stability analysis procedure is demonstrated through the design of a mobile-robot path-following controller.

INTRODUCTION

Control design based on output measurement is very common in engineering practice, often due to the difficulty in measuring all the state variables. In some occasions a static feedback of the output is not sufficient for stabilizing the system or to ensure a satisfactory performance. There are several strategies to deal with this problem. For instance, when dynamic output feedback is considered, one may design an observer to reconstruct the whole state based on the information of the output. Another strategy commonly taken in practice is to use the time derivative of the output signal y(t) for stabilizing the system, which leads to simpler controllers in comparison with observer-based control design. Since the derivative of the output usually can not be measured directly, it is usually realized through delaydifference:

ẏ(t) ≈ y(t) -y(t -τ) τ ( 1
)
where τ is some positive delay value. In consequence of delaydifference approximation (1) used in the feedback, the closed loop system is a delay system with delay-dependent parameters.

The idea of using delay for stabilization is not new. For instance, a multiple delay framework is developed in [START_REF] Niculescu | Stabilizing a chain of integrators using multiple delays[END_REF] for stabilizing a chain of integrators, while in [START_REF] Yamanaka | Use of multiple time-delays as controllers in IMC schemes[END_REF] multiple delays are used for analysing some internal model control scheme. Bounded control for global stabilization has also been addressed in [START_REF] Mazenc | Global asymptotic stabilization for chains of integrators with a delay in the input[END_REF], where only a single delay is used. Our research differs from the previous ones in that we fix the other parameter of the controller while looking for the range Sponsor and financial support acknowledgment goes here. Paper titles should be written in uppercase and lowercase letters, not all uppercase.

of the delay parameter which guarantees that the closed-loop system is exponentially stable with some pre-specified decay rate.

Systems with delay-dependent coefficients appears can be found in biological systems [START_REF] Fabien | Global Asymptotic Stability and Hopf Bifurcation for a Blood Cell Production Model[END_REF], physical systems [START_REF] Wilmot-Smith | A time delay model for solar and stellar dynamos[END_REF] as well general control system as discussed in this paper. A large amount of research effort has been dedicated to stability analysis of delay systems, see [START_REF] Gu | Stability of timedelay systems[END_REF]; Niculescu (2001); [START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems: An Eigenvalue-Based Approach[END_REF] for comprehensive discussion of the related results. However, research pertaining to systems with delaydependent coefficients is not common in the literature. Berreta and Kuang provided an efficient method for analysing stability of such type of systems with a single delay in [START_REF] Beretta | Geometric stability switch criteria in delay differential systems with delay dependent parameters[END_REF]. Gu et al. (2016) relaxed some of their restrictive conditions and extended their approach for more general delay systems. Given I , a closed delay interval of interest, the method presented in Gu et al. (2016) can be used to find all the sub-sets in I that guarantees asymptotical stability of the system. It can be considered as a generalization of the classical τ-decomposition approach (see for instance [START_REF] Michiels | Stability, Control, and Computation for Time-Delay Systems: An Eigenvalue-Based Approach[END_REF]; [START_REF] Lee | On the τ-decomposition method of stability analysis for retarded dynamical systems[END_REF]), which was proposed originally for delay-systems with fixed coefficients.

We will first specify the form of control law considered in this paper and the characteristic equation of the linearized closed-loop system resulting from the control design. Then we will show that by shifting the variable in the characterize equation, the condition for exponential stability with decay rate α is equivalent to a condition for just asymptotical stability. After that we will make some further extension for the method developed in Gu et al. (2016) so that it can be used for the stability analysis of the control system considered in this paper.

MOTIVATING EXAMPLE

Consider a robot path following problem. A unicycle travelling at a constant speed V is required to follow a straight line. The robot is non-holonomic, so the direction of its translational velocity is always along its heading. The control input u is the derivative of its yaw rate, which reflects the yaw moment applied to the robot. It is easy to see the linearized dynamics of the system is described by

ė = V θ θ = ω ω = u
where e stands for the lateral tracking error, θ is the heading angle of the robot and ω is the yaw rate. Variable u denotes the control input. [Put a figure here to illustrate the model] To stabilize the system, one may first render the subsystem consisting only the variable ω stable. Therefore we define the control input as u = -k 0 ω + v where v is the new control input to be designed. Then it is easy to see that there exits real number k 1 , k 2 such that by choosing v = -k 1 ek 2 θ one can stabilize the system. However in practice θ is not convenient to measure. Noticing θ = V -1 ė and ė ≈ e(t)-e(t-τ) τ , we choose instead the following control law that uses delayed signal:

u = -k 0 ω -k 1 y -k 2 y(t) -y(t -τ) τ (2)
It can be shown that if System (2) can be stabilized by the following control law for some fixed number k

0 , k 1 , k 2 u = -k 0 ω -k 1 y -k 2 ẏ
(3) then it can also be stabilized using (2) for sufficiently small delay. Unfortunately if the delay value is too small, the noise contained in the measurement of y will be greatly amplified and injected into the closed-loop system and thus severely undermine the performance. On the other hand too large value of τ may cause slow convergence, strong oscillation, or even instability. Therefore for practical consideration it is useful to find a set of delay value for (2) such that the closed loop system is exponentially stable with some guaranteed decay rate and then one can choose an appropriate delay value in this set.

PROBLEM STATEMENT

Consider a linear system of the form ẋ = Ax + Bu (4) where x ∈ R n is the system state and u ∈ R is the control input. Suppose that the system admits a unique equilibrium and there exits a set of outputs y = Cx that can be measured for control feedback and the following two conditions holds: Condition 1. The relative degree of y w.r.t u is larger than 1, in other words CB = 0. Condition 2. There exits feedback law of the form u(t) = K 0 y(t) + K 1 ẏ(t)

(5) such that System (4) with ( 5) is asymptotically stable. In practice, it is very often that signal ẏ can not be measured. In this case, delay-difference of y(t) can be used to approximate ẏ(t) and ( 5) thus becomes

u(t) = K 0 y(t) + K 1 y(t) -y(t -τ) τ (6)
where τ > 0 is a constant number. For the closed loop system consisting of ( 4) and ( 6), if the system trajectory x t converges to zero exponentially fast with decay rate α, then we say the system is α-stable. Otherwise the system is α-unstable. For fixed K 0 , K 1 and a given delay interval I = (0, τ u ) as well as some positive decay rate α, we are interested in finding all the subintervals contained in I such that the closed loop system is α-stable for all τ in these subintervals. We require that no eigenvalue of the open loop system has a real part equal to -α.

Discussion

Although in the framework described above the feedback that only utilizes y(t) and ẏ(t), however in practice a variety of stabilization problems can be converted into the form described here by first introducing some extra feedback terms in the control law. The basic idea is as follows. Suppose z = C z x are all the outputs that can be measured. One may first define the control law as u = K z z + v (7) and the system dynamics becomes ẋ = (A + BK z C z )x + Bv (8) Then System (8) with v as the new control input may satisfy our conditions and thus fits our framework. This procedure has already been illustrated by our motivating example. As a matter of fact, the classical PID controller can be analysed in this Framework using the same idea. To generate the integral term, one first introduces a state σ that satisfies σ = y Further define u = K I σ + v, then by setting v = K P y + K D ẏ we realize a PID control.

As noted in [START_REF] Niculescu | Stabilizing a chain of integrators using multiple delays[END_REF], when the open loop system possesses more than a pair of imaginary roots, then it is necessary to introduce multiple delays in order to stabilize. We will leave this issue to our future work and restrict ourselves to control system with a single delay in this paper.

STABILITY ANALYSIS

Characteristic Equation and Stability

Let G(λ ) = 0 be characteristic equation of the open loop system (4) with u ≡ 0, then G(λ ) is a polynomial of λ of order n. The characteristic equation of the delay-free system consisting of (4) and ( 5) takes the form

G(λ ) + G u (λ )(k 1 + k 2 λ ) = 0 (9)
where k 1 , k 2 are two constants depending on K 1 as well as K 2 , and G u (λ ) is also a polynomial of λ and its order is smaller than n -1 due to Condition I. When the delayed feedback ( 6) is applied, it is easy to see that the corresponding characteristic equation of the closed loop system becomes

G(λ ) + G u (λ )(k 1 + k 2 1 -e -λ τ τ ) = 0 (10)
It is known that the closed loop system with control law ( 6) is α-stable if and only if all the roots of (10) in λ admit a real part smaller than -α. Define

P α (λ , τ) = G(λ -α) + (k 1 + k 2 τ -1 )G u (λ -α) Q α (λ , τ) = -τ -1 k 2 e ατ G u (λ -α)
It is easy to see that α-stability of the closed loop system with control law ( 6) is now equivalent to the asymptotical stability of a delay system with the following characteristic equation

P α (λ , τ) + Q α (λ , τ)e -λ τ = 0
(11) Our objective can be stated as finding all the delay subintervals contained in I that guarantee that all the the roots of ( 11) in λ are located on L.H.P.

A General τ-decomposition Approach

Characteristic equation ( 11) corresponds to a delay system with delay-dependent coefficients. Stability analysis of systems with this special feature is not common in the literature. Inspired by the earlier work of Beretta and Kuang [BK], the authors proposed a systematic method for analysing the stability of systems with delay-dependent coefficients, which can be viewed as a generalized τ-decomposition approach. In this subsection we will present some related results.

To begin with, define 13) and ( 14) }. Let T be exactly the set of all critical delay values in I such that (11) admits imaginary roots. We assume that both Φ F and T are finite sets and Φ F ∩ T = {φ }. The elements of these two sets are arranged in the ascending order: Φ F = {τ (i) , i = 1, ..., N F } and τ (1) < ... < τ (N F ) T = {τ i , i = 1, ..., N T } and τ 1 < ... < τ N T Notice that both sets can be empty.

F(ω, τ) = P α ( jω, τ)P α (-jω, τ) -Q α ( jω, τ)Q α (-jω, τ) (12) Consider the equation consisting of F(ω, τ) = 0 (13) ∂ ω F(ω, τ) = 0 (14) Define Φ F = {τ|τ ∈ I , ∃ω, s.t (ω, τ) satisfies (
I is thus decomposed into subintervals I (i) , i = 1, ..., T and I (1) = (0, τ (1) ], I (i) = (τ (i-1) , τ (i) ] for i = 1.

It follows from Gu et al. (2016) that within each I (i) , there exits a fixed number m(i) (m(i) can be zero), such that equation ( 13) with (14) admits exactly m(i) positive solutions in ω. Each of these solutions can be viewed as an analytical function of τ within I (i) , denoted as ω

(i) k (τ), 1 ≤ τ ≤ m(i).
These functions are referred to as the critical frequency functions. Now we can introduce the phase angle function θ

(i) k defined in interval I (i) for ω (i) k as θ (i) k (τ) = ∠P α ( jω (i) k (τ), τ) -∠Q α ( jω (i) k (τ), τ) + ω (i) k (τ)τ + π
(15) where ∠ can be any differentiable function in each I (i) that measures the phase angle of a non-zero complex number. The function θ

(i) k (τ), τ)| + |Q α ( jω (i) k (τ), τ)| = 0
(16) To see the statement above is correct, one only needs to notice that if the L.H.S of the last inequality can become zero for some τ * , then we must have G u ( jω i k (τ * )α) = 0, which further implies that λ * = jω i k (τ * )α is an imaginary root of G(λ ) = 0. In other words, the open loop system admits an eigenvalue with real part -α, a contradiction to our requirement on α.

A necessary and sufficient for λ = jω * to be an imaginary root for τ = τ * ∈ I is that the following equation holds

θ (i) k (τ) = 2lπ (17)
for some integer l, k, and i is the number that satisfies τ * ∈ I (i) . By finding the solutions of (17) for all possible integer i, k and l, one is able to obtain the set T .

It can be shown that for some delay value τ * ∈ I i and for some corresponding imaginary root λ = jω * , λ is a differentiable function of τ in a neighbourhood of τ = τ * . We denote it as λ (τ). To indicate whether this root will be come stable or unstable, we define

Inc(τ * , ω * ) = sign ℜ dλ (τ) dτ τ=τ * (18) Inc(τ * ) = m(i) ∑ k=1 Inc(τ * , ω i k (τ * )) (19)
In this work we assume that for any critical delay τ * the right hand side of the last equation is not zero. Then it is easy to see that if Inc(τ * , ω * ) = 1, then this imaginary root moves to R.H.P as τ increases. It moves to L.H.P if the indicator is -1.

Consequently Inc(τ * ) is just the increased number of unstable roots. It is clear that if this assumption is violated, or in other words at some critical delay the derivative of the root w.r.t delay has zero real part, then more involved definition as well as more advanced root crossing criterion than that used in this paper will be needed. More comprehensive discussion of this issue and the related results have recently been developed in Gu et al. (2016).

Once the indicating functions are know for each critical delay, the total number of unstable roots for any τ ∈ I -T , denoted as N u (τ), can be computed as

N u (τ) = N u (τ l ) + i=L(τ) ∑ i=1 Inc(τ i ) (20)
where L(i) is the largest integer such that τ L(τ) < τ. The following theorem provides a convenient way to compute the indicating function in (18). Theorem 1. Suppose for (11) τ * is a critical delay value and λ = jω * is a corresponding imaginary root. Then there exits integer i, k such that τ * ∈ I (i) and ω * = ω (i) k (τ * ). Assume also 13 and 14 do not hold together for (ω * , τ * ). Then λ as a local function of τ denoted as λ (τ), which satisfies

sgn ℜ d dτ λ (τ * ) = sgn ∂ ∂ ω F(ω * , τ * ) ×sgn d dτ θ (i) k (τ * ) (21)
The last theorem indicates an interesting correlation between the imaginary root cross direction and how the corresponding phase function crosses 2lπ for some integer l. One can easily determine the imaginary root cross direction based on the graph of θ

(i)
k (•) for some integer i and k. We now summarize the analysis procedure in several steps:

Step 1. Solve (13) together with ( 14) subject to ω ≥ 0, τ ∈ I to obtain Φ F as well as τ (i) , i = 0, ...N F . I is thus decomposed into each

I (i) = [r (i-1) , r (i) ].
Step 2. In each I (i) , solve ( 17) to find all the critical delay value τ i , i = 1, ..., N T and thus the set T .

Step 3. Compute Inc(τ i ) for each τ i > 0 by plugging ( 21) into (18) and further uses (19).

Step 4. Now for any interval

I * i = [τ i-1 , τ i ],
we can arbitrarily pick a delay value r in its interior and compute N u (r ) via ( 20), then it follows that N u (I * io ) = N u (r ), where I * io = (τ i-1 , τ i ).

Lower-bounding the Delay Interval

The stability analysis we have presented requires the system coefficients to be continuous on the delay interval I . However for the system considered in this paper the coefficients is not bounded in I as τ appears in the denominator. Fortunately it follows from [cite] that for any real c, all the roots of ( 11) with real part larger than c converge to the roots of ( 5) with real parts also larger than c. The direct implication is the existence of some positive number τ ε such that N u (τ) = N u (0) for all τ ∈ (0, τ ε ] and therefore we only need to carry out the analysis on the interval

I = [τ ε , τ u ].
The issue remaining is how to find such a lower bound delay.

For fixed α, define 11) can be written as

P 1 (s, τ) = G(s -α) + G u (s -α)(k 1 + k 2 (s -α)) P 2 (s, τ) = 1 -e j(s-α)τ τ -s + α then (
P 1 (s, τ) + Q(s)P 2 (s, τ) = 0 (22) where we denote Q(s) = k 2 G u (s -α)
We limit τ ε below some positive constant τ up , and denote

c 0 = e ατ -1, c 1 = τ up -e τ up + 1 τ up ( 23
)
Regarding term P 2 , we have

P 2 ( jω, τ) = 1 -e jωτ -jωτ τ + τ -e τ + 1 τ + j sin(ωτ)(e ατ -1) + 1 -cos(ωτ) τ = P 20 ( jω, τ, c 0 ) + c 1
where the expression for P 20 can be determined according to the last equation and the definition of c 1 .

To show there exits τ ε > 0 such that (11) admits no imaginary solution, it is clear that we only need to show that for τ ∈ [0, τ ε ], the following equation does not hold for real ω

G(ω) = |P 2 ( jω, τ)| (24) where G(ω) = |P 1 ( jω, τ)| |Q( jω)| 2 + c 2 (25)
and constant c 2 can be any positive number to ensure that G(ω) is well-defined. For computing a larger τ ε , it is desirable to choose a very small c 2 . In the case where Q( jω) = 0 for ω ≥ 0, c 2 can be set to 0.

Condition 1 indicates that Ord(P 1 ) -Ord(Q) ≥ 2. We can also pick small τ up to make c 0 as small as possible, which also renders |∂ ω P 20 (0, τ, c 0 )| as small as we wish. Then it is easy to see that there exit positive numbers 26) if there exists any, otherwise let ω 1 = 0. Let ω 2 be the maximal solution of G 2 (ω) = (k 2 ω) 2 (27) if there exits any real solutions in ω, otherwise let ω 2 = 0. We note that both ( 26) and ( 27) are polynomial equations and thus easy to solve using numerical methods. Define the piecewise linear function

k 1 , k 2 ≥ k 1 such that G(ω) -c 0 > k 1 ω for ω ≥ 0 and |P 20 ( jω, τ, c 0 )| ≤ k 2 ω for ω ≥ 0. Let ω 1 (τ) be the minimal positive solution of |P 20 ( jω, τ, c 0 )| 2 = (k 1 ω) 2 (
L τ (ω) = k 1 ω 0 ≤ ω ≤ ω 1 (τ) k 2 ω ω 1 (τ) < ω (28)
it is easy to see that τ 1 ω 1 (τ 1 ) = τ 0 ω 1 (τ 0 ) for any τ 1 , τ 0 . Proposition 2. There exits positive numbers τ ub and τ ε ≤ τ ub such that G(ω) > |P 2 ( jω, τ)| (29) For any τ ∈ (0, τ ε ], (11) has no real solution in ω. Furthermore, if (29) holds for ω ∈ [0, ω 2 ], then it also holds for all nonnegative ω.

For any given delay value τ * , the last proposition allows us to check whether τ * can be taken as τ ε by only checking a bounded interval of ω. If (29) holds for τ = τ * ω ∈ [0, ω 2 ], then τ ε can be set as τ * .

Proof. Throughout the proof ω only assumes non-negative value and τ is always positive. We first point out two facts that can check straightforwardly

|P 20 ( jω, τ, c 0 )| < L (τ) (ω), τ > 0 (30) and for 0 < τ 0 < τ L (τ 0 ) (ω) < L (τ) (ω) (31)
We first consider the case ω 1 (1) = 0. This implies that for any 24) can not hold and therefore (11) does not admit imaginary roots for any delay value. Now consider the case ω 1 (1) > 0, we pick τ ε = ω 1 (1) ω 2 . Then L τ ε takes the form

τ > 0 G(ω) -c 1 > L 1 (ω) = L τ (ω) ≥ |P 2 (ω, τ)| and thus G(ω) > |P 2 (ω, τ)| which indicates that (
L τ (ω) = k 1 ω 0 ≤ ω ≤ ω 2 k 2 ω ω < ω 2 then it is straightforward to verify that G(ω) -c 1 > L τ ε (ω).
Since we have as well

L τ ε (ω) ≥ |P 2 ( jω, τ)|, it follows that G(ω) -c 1 > |P 20 ( jω, τ, c 0 )| 2
and thus (24) can not hold. As a result (11) does not have real solutions in ω for τ ∈ (0, τ ε ].

APPLICATION

Now we apply the stability analysis method developed so for for the unicycle model (2) and assume the possible delay value ranges in I = (0, 0.5]. We start with designing the control law of the form (5). Various tools for controlling LTI systems can be employed. For instance, to minimize a quadratic cost function in the time domain, one can use LQR to optimally determine the control parameters. Since the specific way to determine control parameters is irrelevant to our stability analysis procedure, we simple set the eigenvalues of the closed loop system as -2 and -1.5 ± 4 j. Accordingly we obtain the following control parameters: k 0 = 5, k 1 = 97/4, k 2 = 73/2 and the characteristic equation of the closed loop system is

λ 3 + k 0 λ 2 + k 1 λ + k 2 = 0
(32) to shift the eigenvalues of the closed loop system to λ 1 = -1.5 + 4 j, λ 2 = -1.5 -4 j, λ = -2 From the real part of the eigenvalues we deduce that the decay rate of the control system without delay is 1.5. Suppose we require that when the delay-difference approximation is used for feed-back control, the decay rate is no less than 1, hence we pick α = 1. Now replace λ in (32) with sα, we arrive at the explicit expression of (11) for this particular system s 3 + a 0 s 2 + a 1 s + a 2 (τ) + a 3 (τ)e -sτ = 0 (33) where a 0 = k 0 -3, a 1 = 3 -2k 0 , a 2 = k 1 /τ -1 and a 3 = -e τ k 1 /τ from which the expression of P α ( jω, τ) and Q α ( jω, τ) is derived as

P α = (-ω 3 + a 1 ω) j + a 2 -a 0 ω 2 Q α = a3
We pick τ up = 0.05 and define

G(ω) = (-ω 3 + k 1 ω) 2 + (k 2 -k 0 ω 2 ) 2 k 2 (34)
It can be verified that one can pick L 1 = 1.8. Solving the polynomial equation G(ω) 2 = (L 1 ω + c 1 (0.05)) 2 (35) we find the maximal real solution in ω is approximately 7.5752. Therefore we pick ω c = 7.5752 and we only need to ensure the graph of G(ω) has no intersection with that of |P 2 ( jω, τ ε )| for ω ∈ [0, ω c ]. The graphs of G and P 2 for τ = 0.05 is plotted in Fig. 1. Hence we will pick τ ε = 0.05 and will only be concerned about the delay interval I ε = [0.05, 0.5].

F(ω, τ) = ω 6 + (a 2 0 -2 * a 1 )ω 4 + (a 2 1 -2a 0 a 2 )ω 2 +a 2 2 -a 2 3 (36) ∂ ω F(ω, τ) = 6ω 5 + 4(a 2 0 -2a 1 )ω 3 +2(a 2 1 -2a 0 a 2 )ω (37)
Solving the previous two equations together numerically using MATLAB for τ with the constraint ω ∈ R and τ ∈ I , we find only one solution τ (1) ≈ 0.1143. Hence I is decomposed as 1) ] and I (2) = [τ (1) , τ u ]. It can be verified that m(1) = 2 and m(2) = 0. The value of ω

I = I (1) ∪ I (2) and I (1) = [τ ε , τ (
(1) 1 and ω (2) for various τ ∈ I (1) numerically based on (36) and the graph of these two functions are depicted in Fig. 2.

Subsequently we obtain the function θ

(1) 1 (τ) and θ

(1) 2 (τ), the graph of which is presented in Fig. 3.

We see that the only intersection between the phase curves corresponding to θ (1) i , i = 1, 2 and the horizontal lines whose value equals 2lπ for integer l takes place only at τ 1 ≈ 0.0915 when l = 0. Therefore two imaginary roots appear when τ = τ 1 . Since ∂ ω F(ω (1) 2 (τ 1 ), τ 1 ) > 0 and the derivative of θ

(1) 2 (τ 1 ) is positive, we conclude that (33) the two imaginary roots crosses toward the imaginary axis toward R.H.P as τ goes through τ 1 increasingly. No roots of (33) crosses the imaginary roots in Interval I (2) as m(2) = 0. Now it follows from (20) N u (τ) = 0 for τ ∈ [τ ε , τ 1 ) and N u (τ) = 2 for τ ∈ (τ 1 , τ u ]. Consequently, we can claim that for α = 1 the system is α-stable for τ ∈ (τ 1 , τ u ] and not α-stable for τ ∈ (τ 1 , τ u ]. Simulation is carried out using Simulink. We investigate the evolution of the signal e αt e(t) over time by setting τ = 0.09 and 0.11 respectively. It is clear that if the system is α-stable, then e αt e(t) must be convergent, otherwise it diverges.

CONCLUSION

We have addressed the stability analysis for control scheme that uses delay-deference for approximating the derivative of output signals. We regard the delay as a design parameter and for any given bounded delay interval of interest, we are concerned with finding all the sub-sets of delay-values contained in the interval such that the system is exponentially stable with guaranteed convergence rate. An approach developed in our previous work for systems with delay-dependent coefficients is modified and applied to solve this problem. The stability analysis procedure is illustrated with a mobile robot path tracking problem. Simulation shows that our conclusion drawn from frequency-domain analysis is consistent with the time-domain response of the robot trajectory.

  well defined because Condition 1 and Condition 2 guarantee that |P α ( jω

Fig. 1 .

 1 Fig. 1. The Graphs of function G and |P 2 | have no intersection for τ ε = 0.05 It follows that

  Fig. 2. Frequency functions.

  Fig corresponds to τ = and we observe the convergence of e αt e(t) converges and Fig corresponds to τ = and the divergence of e αt e(t) is observes. These results are all consistent with our previous analysis.

Fig. 4 .Fig. 5 .

 45 Fig. 3. Phase functions.