Silviu-Iulian Niculescu

Islam Boussaada

Counting Characteristic Roots of Linear Delay Differential Equations. Part II: from Argument Principle to Rightmost Root Assignment Methods

This chapter briefly presents some "user-friendly" methods and techniques (frequency-domain approaches) for the analysis and control of linear dynamical systems in presence of delays. The presentation is as simple as possible, focusing more on the main intuitive (algebraic, geometric) ideas to develop theoretical results, and their potential use in practical applications. To fix better the ideas, scalar and second-order examples are largely discussed. Next, a particular attention will be paid to the existing links between the maximal allowable multiplicity of the characteristic roots and the spectral abscissa of the dynamical system. The underlying property -multiplicity induced dominancy -will be particularly useful in constructing low-complexity controllers by partial pole placement. Such an idea is particularly exploited in vibration control.

Introduction

Propagation, transport and heredity represent some of the typical phenomena of physical, biochemical processes and properties of living organisms that can be modeled by using delays in their mathematical representation. Understanding the ways the delays and the other parameters may affect the models' dynamics is a problem of recurring interest during the last decades. This chapter completes the contribution by [START_REF] Niculescu | Counting characteristic roots of linear delay-differential equations. part I: frequency-sweeping stability tests and applications[END_REF] (same volume) devoted to the stability analysis of linear dynamical time-delay systems described by Delay Differential Equations (DDEs) in frequency-domain. More precisely, we will present different (frequencydomain) techniques to count the characteristic roots located in the right complex half-plane with particular emphasis on the Pólya-Szegö results1 and the so-called multiplicity induced dominancy and its potential applications in control engineering.

In the 30s, the development of feedback amplifiers is at the origin of most of frequency-response methods in control engineering and signal processing. In this frame, [START_REF] Nyquist | Regeneration theory[END_REF] proposed a graphical method for the analysis of Single-Input/Single-Output (SISO) Linear Time-Invariant (LTI) systems to determine the (asymptotic) stability of the closed-loop system based on the transfer function of the open-loop system. Due to its simplicity, such a method became extremely popular, and it is at the origin of most of the graphical (stability) approaches and tests in frequency-domain in the open literature, see, for instance, [START_REF] Macfarlane | Frequency-response methods in control systems[END_REF] and the references therein. It should be mentioned that the Cauchy's argument principle is a simple and elegant way to prove Nyquist criterion. However, inspired by some ideas about the propagation of sinusoidal signals through the systems, Nyquist came up with a different proof. It should be mentioned that the seminal work of [START_REF] Tsypkin | Stability of systems with delayed feedback[END_REF] includes an extension of the Nyquist criterion in the case of dynamical systems with input delay and the proof is based on the Cauchy's argument principle.

In the 50s-70s, except for some natural extensions of Nyquist diagrams to delay dynamical systems, the Mikhailov criterion2 has been the source of further developments (see, for instance, the discussion in [START_REF] Barker | Mikhailov stability criterion for time-delayed systems[END_REF] and the references therein). For some connections between Mikhailov and Nyquist tests applied to delay systems we mention the almost forgotten paper by [START_REF] Chen | A general frequency stability criterion for multi-input-output, lumped and distributed-parameter systems[END_REF]. Other graphical criteria include the well-known rootlocus methods3 and the Satche's diagrams (see, e.g. [START_REF] Satche | Discussions on "stability of linear oscillating systems with constant time lag[END_REF], where the latter can be interpreted as a variant of the Nyquist criterion, sometimes called the dual root-locus methods. For a pedagogical presentation as well as some extensions of the root locus methods and Nyquist criterion to deal with dynamical systems including one delay in the input/output channel, we refer to [START_REF] Krall | Stability techniques for continuous linear systems[END_REF]. Finally, for an historical perspective on the frequency-response methods in control area, we refer to [START_REF] Macfarlane | Frequency-response methods in control systems[END_REF].

At the end of the 70s, using the Cauchy's argument principle, [START_REF] Stépán | On the stability of linear differential equations with delay[END_REF] proposed a simple criterion to count the unstable roots of a linear DDE. The simplicity of the method has been illustrated in [START_REF] Stépán | Retarded Dynamical Systems: Stability and Characteristic Functions[END_REF] by constructing stability charts for several examples from mechanical engineer-ing and represented by DDEs. The underlying idea has been extended by [START_REF] Hassard | Counting roots of the characteristic equation for linear delay-differential systems[END_REF] to also cover some cases where multiple characteristic roots exist on the imaginary axis. This method is briefly discussed in Section 2 and applied to a second-order DDE including a single delay.

Exploration of some earlier ideas present, e.g., in [START_REF] Pinney | Ordinary difference-differential equations[END_REF] and the seminal work [START_REF] Hayes | Roots of the transcendental equation associated with a certain difference-differential equation[END_REF] demonstrated that the spectral values of systems of DDEs reaching their maximal possible multiplicity tend to be dominant, in what has come to be known as the Multiplicity Induced Dominancy (MID) property. The MID property turns out to open some interesting perspectives in designing low-complexity feedback control laws for delay systems. Since these seminal works, many research efforts have been devoted to the characterization of the classes of systems for which such a property is valid and to its exploration in applications to stabilize delay systems. This chapter addresses some of the most recent results in this direction.

To the best of the authors' knowledge, it should be mentioned that most general result involving the MID property so far is that of [START_REF] Boussaada | The generic multiplicityinduced-dominancy property from retarded to neutral delay-differential equations: When delay-systems characteristics meet the zeros of Kummer functions[END_REF], which shows that the MID property holds for retarded and neutral DDEs of an arbitrary order with a single delay and highlights the potential of applicability of such a property in the design of control feedback laws by using a partial pole placement methodology. Such results are obtained by exploring links4 between spectra of systems of DDEs and roots of a family of confluent hypergeometric functions, known as Kummer functions.

The contributions of the chapter are fourfold. First, our intention is to complete the presentation of [START_REF] Niculescu | Stability, delays and multiple characteristic roots in dynamical systems: A guided tour[END_REF] by introducing other frequency-domain methods useful to count the number of unstable characteristic roots. Second, we wish to shed some light on various methods, techniques and ideas developed by the authors of this contribution to deal with the case of multiple characteristic roots and, in particular, the so-called MID property. Finally, we are interested to emphasize the way that such a property can be applied to some control problems and, more precisely, in vibration control.

The chapter is organized as follows: some preliminary results and prerequisites are introduced in Section 2. Section 3 is devoted to the characterization of the multiplicity of the characteristic roots in terms of structured matrices. Next, Section 4 introduces the MID property and gives some insights on the resulting rightmost spectral value assignment on comprehensive examples, opening perspectives in control design. A parametric characterization of the MID property in the case of a second-order linear DDE is provided in Section 5. Next, the effect of the maximal multiplicity of a spectral value on the distribution of the remaining spectrum is discussed in Section 6, where some interesting links with Kummer hypergeometric functions are emphasized. A software dedicated to the ensuing partial pole placement, called P3δ is described in 7. Section 8 illustrates the use of the proposed control strategy in the active damping of vibrations occurring in a flexible mechanical structure. Some notes and comments end the chapter.

Notations. In this chapter, we use the following standard notations5 : Z (N) denote the set of (non negative) integers, R (C) denotes the set of real (complex) numbers. Next, R + (R -) denotes the set of strictly positive (negative) real numbers and R * = R \ {0}. For a given complex number λ, (λ) ( (λ)) denote its real (imaginary) part. The open left (right) complex half-plane is the set

C -(C + ) defined by C -= {λ ∈ C | (λ) < 0} (C + = {λ ∈ C | (s) > 0})
. For a given complex number λ the modulus is defined by |λ| = (λ) 2 + (λ) 2 . For a given discrete set χ, card(χ) denotes the number of elements of χ. The symbol sgn designates the sign of a real valued function f at a given argument. Given two vectors u

= (u 1 , . . . , u n ) T ∈ R n and v = (v 1 , . . . , v m ) T ∈ R m we adopt the nota- tion: (u/v) = (u 1 , . . . , u n , v 1 , . . . , v m ) T . For a set K of complex numbers, denote by M m,n (K) (M n (K)) the set of m × n (n × n) matrices with co- efficients in K. Given matrices M 1 ∈ M m,n1 (R) and M 2 ∈ M m,n2 (R) the notation [M 1 M 2 ]
designates the matrix obtained by concatenation. Next, for a vector u ∈ R n , u 2 denotes its 2-norm and for a matrix A ∈ M n (R), its 2-norm is defined by:

A 2 = sup u =0 Au 2 / u 2 . Given k, n ∈ N with k ≤ n, the binomial coefficient n k is defined as n k = n! k!(n-k)!
and this notation is extended to k, n ∈ Z by setting n k = 0 when n < 0, k < 0, or k > n. For α ∈ C and k ∈ N, (α) k is the Pochhammer symbol for the ascending factorial, defined inductively as (α) 0 = 1 and (α) k+1 = (α + k)(α) k . Finally, deg(P ) denotes the degree of P ∈ R[x].

Preliminaries and Prerequisites

Let Ω ⊂ R np be an open set and consider two matrix functions A

0 , A 1 : Ω → M n (R), - → p ∈ Ω → A j ( - → p ) (j = 0, 1)
, that are assumed to be sufficiently smooth. With these notations, consider a dynamical system described by a DDE of retarded type with a single delay of the form:

ξ = A 0 ( - → p )ξ(t) + A 1 ( - → p )ξ(t -τ ), (1)
where ξ = (ξ 1 , . . . , ξ n ) T ∈ R n and ( -→ p , τ ) ∈ Ω × R + denote the statevector and the system's parameters, respectively. Assume further that, for a fixed delay value τ ∈ R + , for almost all -→ p ∈ Ω, under appropriate initial conditions6 , the system above is properly defined7 . It is well known that the asymptotic behavior of the solutions of ( 1) is determined from the location of the spectrum8 χ of the characteristic matrix function9 and defined by M :

C × Ω × R + → C n×n with M(λ, - → p , τ ) = λ I -A 0 ( - → p ) -A 1 ( - → p ) e -τ λ . (2)
Furthermore, these eigenvalues are the roots of the associated characteristic function which is the quasipolynomial

Q : C × Ω × R * + → C Q(λ, - → p , τ ) = det M(λ, - → p , τ ). ( 3
)
Define now the sets χ + , χ 0 , χ -⊂ χ as follows:

χ + = {λ ∈ C, Q(λ, p, τ ) = 0, (λ) > 0}, χ -= {λ ∈ C, Q(λ, p, τ ) = 0, (λ) < 0} and χ 0 = {λ ∈ C, Q(λ, p, τ ) = 0, (λ) = 0}
which give a partition of the spectrum χ of the DDE (1) with respect to the imaginary axis. More precisely, the spectrum χ can be divided into χ = χ + ∪ χ 0 ∪ χ -. This chapter focuses on the computation of card(χ + ), that is the number of unstable characteristic roots. If the only parameter of interest is the delay τ , than Q(λ, p, τ ) simply rewrites as Q(λ, τ ) and (1) becomes:

ξ = A 0 ξ(t) + A 1 ξ(t -τ ), (4)
with A 0 , A 1 constant real valued matrices in M n (R). The proof of the proposition below is given in [START_REF] Michiels | Stability and stabilization of time-delay systems[END_REF].

Proposition 2.1. If λ is a characteristic root of system (4), then it satisfies

|λ| ≤ A 0 + A 1 e -τ λ 2 .
(5)

The inequality (5) combined with the triangular inequality provides a generic envelope curve around the characteristic roots corresponding to the DDE (1), see for instance [START_REF] Mori | Stability of ẋ(t) = Ax(t) + Bx(t -τ )[END_REF] for further insights on spectral envelopes for DDEs with a single delay.

As in the previous chapter by [START_REF] Niculescu | Counting characteristic roots of linear delay-differential equations. part I: frequency-sweeping stability tests and applications[END_REF], most of the examples considered in the sequel are represented by a particular class of quasipolynomial functions given by: Q(λ, p, τ ) = P 0 (λ, p) + P 1 (λ, p) e -λτ , where deg(P 0 ) > deg(P 1 ), and a particular attention will be paid to scalar and second-order DDEs. Notice that some generic result on the location of associated spectral values for arbitrary polynomials P 0 and P 1 satisfying the condition deg(P 0 ) ≥ deg(P 1 ) can be found in [START_REF] Boese | Stability with respect to the delay: On a paper of K. L[END_REF]; see also [START_REF] Michiels | Stability and stabilization of time-delay systems[END_REF]. Furthermore, in most of the cases, the parameter vector p is explicitly defined by the coefficients of the polynomials P 0 and P 1 , assumed to be independent one another.

In the following paragraphs, we present a characterization of the admissible/maximal spectral values' multiplicities10 and describe their effect on the stability of the trivial solution as well as the corresponding exponential decay rate, see for instance [START_REF] Mori | On an estimate of the decay rate for stable linear delay systems[END_REF]. As a byproduct of the approach, the main steps of an algorithmic procedure are proposed. We think that such an idea can be extended to more general quasipolynomials.

On Integration Contours for Quasipolynomials Corresponding to DDEs of Retarded Type

Cauchy's argument principle is a basic complex analysis property11 widely used in the stability analysis of LTI dynamical systems, see for instance [START_REF] Marden | Geometry of Polynomials[END_REF]; [START_REF] Ahlfors | Complex Analysis[END_REF]. Roughly speaking, it relates the difference between the number of zeros and poles of a meromorphic function f in a simply connected domain D ⊂ C to a contour integral of the function's logarithmic derivative on the boundary ∂D; note that, this is also equal to the winding number of the curve ∂D with respect to f . Theorem 2.2 [START_REF] Ahlfors | Complex Analysis[END_REF]. Let D be a simply connected region with boundary Γ (piecewise smooth and oriented anti-colockwise). Let f be a meromorphic function in an open set containing the closure D with poles p 1 , . . . , p l and zeros z 1 , . . . , z m counted according to their multiplicity, none of which belonging to the closed curve Γ. Then

1 2iπ Γ f (s) f (s) ds = Z -P, g 2 g 3 g 1 R R -R x ω 0 s 0 R R g 1 g 2 g 3 g 4 -R
x ω 0 Figure 1. LTI dynamical system; generic contours for applying the argument principle to investigate qualitative properties. (Left) Standard Bromwich contour usually adopted for asymptotic stability investigation;

(Right) a contour often used to investigate ω-stability or the dominancy of a given root. Figures are borrowed from [START_REF] Boussaada | Multiplicity-induced-dominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF].

where Z and P designate respectively the number of zeros and the number of poles of f enclosed by Γ.

In particular, note that if f is a holomorphic function12 and Γ is a closed piecewise C 1 curve then

1 2iπ Γ f (s) f (s) ds = Z.
Moreover, standard contours like the modified Bromwich contour are often used to explore the asymptotic stability of dynamical systems' solutions, as it allows to count the zeros with (s) > 0 (see, e.g., Figure 1). Other contours can be chosen based on the inherent properties of the considered characteristic function (Figure 2). Practically, owing to Proposition 2.1, a generic supremum bound for the real and imaginary parts of the roots of function ( 30) is established. Then, one may define an integration contour Γ = ∪ l k=1 g k which is not passing through zeros of Q and is taken as a counterclockwise closed curve. Hence, the contour integral over Γ is defined as the sum of the integrals over the directed smooth curves that make up Γ, as illustrated in Figure 2. Boussaada and Niculescu (2018) for applying the argument principle to investigate the dominancy of multiple roots; dashed red curve: generic spectrum envelope given by Proposition 2.1 for the quasipolynomial Q(λ, p, τ ) = λ + a 0 + a 1 e -λτ under conditions λ 0 = -a0τ +1 τ and a 1 = e λ 0 τ τ with a 0 = τ = 1. The figure is borrowed from [START_REF] Boussaada | Multiplicity-induced-dominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF].

Since we are concerned with zeros of quasipolynomials13 , the following result14 is a direct consequence of the argument principle applied to such functions.

Theorem 2.3 [START_REF] Bellman | Differential-difference equations[END_REF]. If f is an analytic function inside and on a closed contour Γ, and is not zero on Γ, then the number of zeros (counting their respective multiplicities) of f within Γ is equal to 1/2π times the variation of the argument of f (s) as s moves once around Γ in the counterclockwise sense.

Example 2.4. Consider the scalar DDE:

ξ(t) + a ξ(t) + b ξ(t -τ ) = 0, (6)
where p = (a, b) T and τ are such that τ > 0 and 0 < b < a. To study the asymptotic stability of the trivial solution of (6), one investigates the distribution of zeros of the characteristic function given by:

Q(λ, p, τ ) = λ + a + b e -λ τ . ( 7
)
To apply the argument principle, consider the standard Bromwich contour illustrated in Figure 1. First, it is easy to observe that Q does not admit imaginary roots. Indeed, by assuming Q(i ω, p, τ ) = 0, one gets:

a = -b cos(ω τ ) ω = b sin(ω τ ), (8)
where the first condition in ( 8) is inconsistent since a > b. Next, as emphasized in Proposition 2.1, in C + , the spectrum is finitely bounded by an appropriate envelope, so that Q does not admit zeros on the semi-circle C R when R → ∞. Furthermore, since, in C + , the zeros of Q(λ, p, τ ) coincides with the ones of Q(λ, p, τ )/(λ + a) = 1 -Q(λ, p, τ ) where Q(λ, p, τ ) = -b e -λ τ λ+a , then one may investigate the values of λ for which Q(λ, p, τ ) = 1. Notice that the limiting position as R → ∞ of the contour Γ under the mapping Q(λ, p, τ ) is called the amplitude-phase characteristic, see for instance [START_REF] El | Introduction to the theory and application of differential equations with deviating arguments[END_REF].

To apply the argument principle or equivalently the argument variation on Q(λ, p, τ ), one needs to calculate the number of circuits of the amplitudephase characteristic around the point λ = 1 (rather than λ = 0). Note also that under the mapping Q(λ, p, τ ), the semi-circle C R when R → ∞ reduces to a point. So that, one needs only to inspect the mapping Q(λ, p, τ ) through the imaginary axis traversed in the negative direction. When R → ∞ the limiting characteristic of Q(λ, p, τ ) is nothing but the rational transformation -b/(λ+a), which transforms the imaginary axis into the circle of center λ = -b/2a and radius b/a. Since a, b > 0, then the mapping Q(λ, p, τ ) does not make any circuit around λ = 1, so that Q has no zeros in C + .

Another way to obtain the same conclusion of Example 2.4 is to employ the Rouché's Lemma.

Theorem 2.5 [START_REF] Marden | Geometry of Polynomials[END_REF]. If f and g are analytic function inside a simple closed Jordan curve Γ and if they are continuous on Γ with

|f (λ)| < |g(λ)|,
then the function f + g has the same number of zeros of g inside Γ.

Let us revisit Example 2.4 by applying Rouché's Lemma. To do so, let us define the functions g(λ) = λ + a and f (λ) = b e -λτ . Since 0 < b < a, the condition of Rouché's Lemma are satisfied on the contour Γ defined as the standard Bromwich contour presented in Figure 1, that is

|b e -λτ | < |λ + a|.
Hence, Q defined by ( 7) does not admit any roots in C + .

When a given quasipolynomial function corresponding to a DDE of retarded type has no roots on the imaginary axis, then the seminal Stépán's work [START_REF] Stépán | On the stability of linear differential equations with delay[END_REF] gives a new formalism and an easy procedure to characterize the exact number of unstable roots, that is the roots located in χ + , see also Stépán (1989, Theorems 2.15-2.16). The proof of such a result is mainly based on the argument principle. Motivated by a potential application in bifurcation theory, the main theorem from [START_REF] Hassard | Counting roots of the characteristic equation for linear delay-differential systems[END_REF] 15 relaxes the assumption card(χ 0 ) = ∅. Thus, it emphasizes the link between card(χ + ) and card(χ 0 ), both taking into account the multiplicity.

Theorem 2.6 (Hassard, 1997, p. 223). Consider the quasipolynomial function Q defined by

Q(λ, τ ) = P 0 (λ) + N k=1 P k (λ) e -τ k λ . with τ = (τ 1 , . . . , τ N ) T . Let ρ 1 , . . . , ρ r be the positive roots of R(y, τ ) = (i n Q(i y, τ ))
, counted by their multiplicities and ordered so that 0 < ρ 1 ≤ . . . ≤ ρ r . For each j = 1, . . . , r there exist τ j such that Q(i ρ j , τ j ) = 0, assume that the multiplicity of iρ j as a zero of Q(λ, τ ) is the same as the multiplicity of ρ j as a root of R(y, τ ). Then card(χ + ) is given by the formula:

card(χ + ) = n -card(χ 0 ) 2 + (-1) r 2 sgn S (µ) (0, τ ) + r j=1 sgn S(ρ j , τ j ),
where µ designates the multiplicity of the zero spectral value of Q(λ, τ ) = 0 and S(y, τ ) = (i -n Q(iy, τ )). Furthermore, card(χ + ) is odd (respectively, even) if

Q (µ) (0, τ ) < 0 (Q (µ) (0, τ ) > 0)
. If R(y, τ ) = 0 has no positive zeros (in y), set r = 0 and omit the summation term in the expression of card(χ + ). If s = 0 is not a root of the characteristic equation, set µ = 0 and interpret S (0) (0, τ ) as S(0, τ ) and Q (0) (0, τ ) as Q(0, τ ).

In order to illustrate the above result, we count the number of unstable roots of the following quasipolynomial with delay-dependent coefficients:

Q(λ, τ ) = λ 2 - 4 λ τ + 6 τ 2 - 6 τ 2 + 2 λ τ e -λτ . (9)
Notice that (9) admits a characteristic root at the origin with multiplicity four.

Characterizing imaginary roots of the quasipolynomial. In order to apply the argument principle on the standard Bromwich contour, allowing to count the roots of the quasipolynomial (9) in C + , a deflation which eliminates the roots on the imaginary axis is required. To do so, we first investigate nonzero imaginary roots of (9). Assuming that there exists ω > 0 such that λ = iω is a root of ( 9), we define

R(ω) = (i -2 Q(iω, τ )) and S(ω) = (i -2 Q(iω, τ )), so that      R(ω) = 6 cos (ω τ ) τ 2 + 2 ω sin (ω τ ) τ + ω 2 - 6 τ 2 , S(ω) = 2 ω cos (ω τ ) τ -6 sin (ω τ ) τ 2 + 4 ω τ . ( 10
)
As a consequence, for any λ = iω, root of ( 9), we have R(ω) = S(ω) = 0. Then, some algebraic manipulations allow to eliminate the trigonometric functions as follows:

       cos (ω τ ) = - 1 2 7 ω 2 τ 2 -18 ω 2 τ 2 + 9 , sin (ω τ ) = - 1 2 ω τ ω 2 τ 2 -18 ω 2 τ 2 + 9 .
Using the standard trigonometric identity cos 2 (ω τ ) + sin 2 (ω τ ) = 1, one obtains exclusively 0 and the non-vanishing solutions ω = ± 3 i τ , which are discarded since we are dealing with positive frequencies. By using the notations from the Theorem above where card(χ 0 ) stands for the number of imaginary roots of ( 9) and κ designate the multiplicity of 0. In our case, card(χ 0 ) = κ = 4. Hence, the deflated function which is integrated on the Bromwich contour B is given by

Q(λ, τ ) = Q(λ, τ ) λ 4 ,
It has the same nontrivial zeros encircled by B as the quasipolynomial Q.

Characterizing ω k the positive roots of R and exploring sgn(S(ω k )). Consider now the positive roots of R. If ω > 0 is a root of the first transcendental function of (10), then

6 cos (ω τ ) τ 2 + 2 ω sin (ω τ ) τ + ω 2 - 6 τ 2 = 0.
Multiplying the last equation by τ 2 reduces the problem to the search of zeros of an univariate function:

F (ρ) = 6 cos (ρ) + 2 ρ sin (ρ) + ρ 2 -6 (ρ = ω τ ).
which is equivalent to state that, for a given positive delay τ , there exists a unique frequency ω ∈ R + such that ρ = ω τ and R(ω , τ ) = 0.

By similarity, define G(ρ) = τ 2 S(ω). Then,

G(ρ) = 2 ρ(2 + cos (ρ)) -6 sin (ρ)
which is positive in the interval (π, 2π) since sin (ρ) < 0. In particular, G(ρ , τ ) > 0 and thus S(ω ) > 0.

Counting the number of roots dominating zero. We are able now to compute the number of unstable roots of ( 9) based on the computation of the argument variation. As a matter of fact, combining all the above collected information on the behavior of both the real R and the imaginary S parts of the quasipolynomial function ( 9), we apply directly the theorem above. Namely, Z designates the number of roots dominating 0:

Z = n -card(χ 0 ) 2 + 1 2 (-1) r sgn S (κ) (0) + r j=1 (-1) j-1 sgn S(ω j ), ( 11
)
where n is the order of the system, κ is the multiplicity of 0 as a root of (9), K is the total number of roots of ( 9) on the imaginary axis and r is the number of positive real roots of R.

In our case, n = 2, κ = card(χ 0 ) = 4, r = 1 and S (4) (ω) = 2 sin (ω τ ) τ 2 +2 ω cos (ω τ ) τ 3 which vanishes at ω = 0.
As a result, we easily obtain that Z = 0 owing to formula (11). Hence, the dominancy of 0 as a root of ( 9) is established; that is no roots with positive real part exist for the quasipolynomial Q.

Pólya-Szegö Theorem: Counting Quasipolynomial Roots

The following result was first introduced and claimed in the problems collection published in 1925 by G. Pólya and G. Szegö. In the fourth edition of their book16 , the authors emphasize that the proof was obtained by N. Obreschkoff in 1928 using the argument principle, see [START_REF] Obreschkoff | Nullstellen linearer kombinationen von exponentialfunktionen[END_REF]. Such a result gives a bound for the number of quasipolynomial's roots in any horizontal strip. As a consequence, a bound for the number of quasipolynomial's real roots can be easily deduced.

Proposition 2.7 ( Pólya and Szegő, 1972, pp. 144)

. Let τ 1 , . . . , τ N denote real numbers such that τ 1 < τ 2 < • • • < τ N , and d 1 , . . . , d N positive integers satisfying d 1 ≥ 1, d 2 ≥ 1, . . . , d N ≥ 1 and define D such that D = d 1 + d 2 + • • • + d N -N.
Let f i,j (λ) the function f i,j (λ) = λ j-1 e τi λ , for 1 ≤ j ≤ d i and 1 ≤ i ≤ N . For α, β two real numbers, let P S be the number of zeros of the complex valued function

f (λ) = 1≤i≤N,1≤j≤di c i,j f i,j (λ) contained in the horizontal strip {λ ∈ C, α ≤ (λ) ≤ β}. Assuming that 1≤k≤d1 |c 1,k | > 0, . . . , 1≤k≤d N |c N,k | > 0, then (τ N -τ 1 ) (β -α) 2 π -D + 1 ≤ P S ≤ (τ N -τ 1 ) (β -α) 2 π + D + N -1.
Setting α = β = 0, the above Proposition allows to P S ≤ D + N -1 where D stands for the sum of the degrees of the polynomials involved in the quasipolynomial function f and N designates the associated number of polynomials. As explained in the sequel, in the case of complete polynomials, this gives a sharp bound.

Characterizing Multiplicity Using Structured Matrices

Matrices arising from a wide range of problems in mathematics and engineering typically display a characteristic structure. Exploiting such a structure is the mean to the design of efficient algorithms, see for instance [START_REF] Bini | A fast algorithm for approximate polynomial GCD based on structured matrix computations[END_REF]. The discussion below is a crossroad between the investigation of a class of such structured matrices originally involved in multivariate interpolation problems 17 and the estimation of the upper bound for the codimension of spectral values of linear DDEs (which are the zeros of the corresponding characteristic quasipolynomial ). The aim of this section is threefold: firstly, it emphasizes the link between the above two quoted issues. Secondly, it shows that the codimension of the zero as well as purely imaginary spectral values of a given DDE are characterized by some algebraic properties of an appropriate functional Birkhoff matrix. Finally, it shows the effectiveness of the proposed constructive approach by exploring the generic settings as well as investigating some specific but significant sparsity patterns. In both cases, symbolic algorithms for LU-factorization are presented for some novel classes of Birkhoff matrices. It is worth mentioning that such an attempt can be exploited for further classes of Birkhoff matrices and should be of interest in some linear algebra problems involving structured matrices as well as in applications including polynomial interpolation. The proofs of the results presented in this section can be found in Boussaada and Niculescu (2016a,b).

Singularity Codimension May Exceed the Model Order:

Bogdanov-Takens Singularity

Consider first a simple scalar DDE with one delay representing a biological model discussed in [START_REF] Cooke | Stability analysis for a vector disease model[END_REF] describing the dynamics of disease propagation. Namely, the infected host population x(t) is governed by:

ẋ(t) + a 0 x(t) + a 1 x(t -τ ) -a 1 x(t -τ ) x(t) = 0,
where a 1 > 0 designates the contact rate between infected and uninfected populations and it is assumed that the infection of the host recovery proceeds exponentially at a rate -a 0 > 0; see also [START_REF] Ruan | Delay differential equations in single species dynamics[END_REF] for more insights on the modeling and stability results. The linearized system is given by

ẋ(t) + a 0 x(t) + a 1 x(t -τ ) = 0, ( 12
)
where p = (a 0 , a 1 ) T with (a 0 , a 1 , τ

) ∈ R 2 × R + , then the associated charac- teristic (transcendental) function Q becomes Q(λ, p, τ ) = λ + a 0 + a 1 e -λτ . (13)
Zero is a spectral value for (12) if, and only if, Q vanishes at zero which is equivalent to a 0 + a 1 = 0. The computations of the first derivatives of (13) with respect to λ give:

Q (λ, p, τ ) = 1 -τ a 1 e -λτ , Q (λ, p, τ ) = τ 2 a 1 e -λτ .
If additionally a 1 = 0, then the codimension of the zero spectral value is at most two18 since the algebraic multiplicity 2 is ensured for τ * = 1/a 1 > 0, a 0 = -a 1 and Q (0, p , τ ) = 0. It should be mentioned that Q includes 3 parameters (a 0 , a 1 , τ ) and the codimension cannot exceed 2. It is easy to observe that the Pólya-Szegö bound P S = 2 if if one chooses19 α = β = 0.

Codimension Counting: a Vandermonde/Birkhoff-Based Framework

Consider now that the DDE (4) includes N discrete (constant) delays where the vector ζ = (ζ 1 , . . . , ζ n ) ∈ R n denotes the state-vector, then the DDE reads

ζ = N k=0 A k ζ(t -τ k ), ( 14
)
where τ k , k = 1, . . . , N are strictly increasing positive constant delays such that τ 0 = 0 and τ = (τ 1 , . . . , τ N ) T , and A k ∈ M n (R) for k = 0, . . . , N . In this case, the characteristic function of system ( 14) reads

Q(λ, p, τ ) = det λ I -A 0 - N k=1 A k e -τ k λ , (15)
where p denotes the entries of the matrices A k , for k = 0, . . . , N . For a given delay vector τ and vector parameter p, the corresponding spectrum χ can be split into

χ = χ + ∪ χ 0 ∪ χ -where χ + = {λ ∈ C, Q(λ, p, τ ) = 0, (λ) > 0}, χ -= {λ ∈ C, Q(λ, p, τ ) = 0, (λ) < 0} and χ 0 = {λ ∈ C, Q(λ, p, τ ) = 0, (λ) = 0}.
Let us start by setting a new parameterization for the quasipolynomial function (15) of the DDE ( 14) and defining some useful notations adopted through this section. Some straightforward computations give the following formal expression of the quasipolynomial function ( 15)

Q(λ, p, τ ) = P 0 (λ, p) + M k ∈S N,n P M k (λ, p) e σ M k λ , (16)
where σ M k = -M k τ and S N,n is the set of all the possible row vectors

M k = (M k 1 , . . . , M k N ) belonging to {1, . . . , n} N such that 1 ≤ M k 1 + . . . + M k N ≤ n.
Furthermore, by running the index from 1 to the cardinality ÑN,n := card(S N,n ), ( 16) can be written in the following compact form:

Q(λ, p, τ ) = P 0 (λ, p) + ÑN,n k=1 P k (λ, p) e σ k λ . ( 17
)
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. In this case, one has M 2 = (0, 1, 0) and Ñ3,2 = 9. A generic property of DDEs ( 14) allows considering P 0 as a monic polynomial of degree n in λ and the polynomials P M k satisfying deg(

P M k ) = n - N s=1 M k s ≤ (n -1) ∀M k ∈ S N,n .
In the sequel, P 0 (•, p) will be called the delay-free polynomial and the quasipolynomial function λ → ÑN,n k=1 P M k (λ, p) e σ k λ will be called the transcendental part of the quasipolynomial.

Next, define a j,k as the coefficient of the monomial λ k for the polynomial P M j , 1 ≤ j ≤ ÑN,n , and denote P M 0 = P 0 . Thus, a 0,n = 1 and a j,k = 0

∀k ≥ d j = n- N s=1 M j s .
Here, d j -1 is nothing but the degree of P M j . Furthermore, we denote by a 0 = (a 0,0 , a 0,1 , . . . , a 0,n-1 ) T the vector of the coefficients of the delay-free polynomial and by a j = (a j,0 , a j,1 , . . . , a j,dj -1 ) T the vector of the coefficients of the polynomial associated to the auxiliary delay σ j for 1 ≤ j ≤ ÑN,n . Next, set the delay auxiliary vector

σ = (σ 1 , σ 2 , . . . , σ ÑN,n ) and a = (a 1 / a 2 / • • • / a ÑN,n ) T .
This leads to a new parameterization of the quasipolynomial Q(λ, p, τ ) = Q(λ, a, τ ). Finally, let us denote by Q (k) the k-th derivative of Q with respect to the variable λ. We say that zero is an eigenvalue of algebraic multiplicity m ≥ 1 for ( 14) at a = a * and τ

= τ * if Q(0, a * , τ * ) = Q (k) (0, a * , τ * ) = 0 for all k = 1, . . . , m -1 and Q (m) (0, a * , τ * ) = 0. We assume also in what follows that σ k = σ k for any k = k where k, k ∈ S N,n
. Indeed, if for some value of the delay vector τ there exists some k = k such that σ k = σ k , then the number of auxiliary delays and the number of polynomials is reduced by considering a new family of polynomials P satisfying PM k = P M k + P M k .

In the sequel, D q = ÑN,n k=0 deg(P M k ) + ÑN,n -1 will designate the degree of the transcendental part of the quasipolynomial 20 . Now, to characterize the structure of a given quasipolynomial function one needs to introduce a vector V, called incidence vector, vector that reproduces the data on the vanishing components of the vector a defined above. Thus, V is a sparsity patterns indicator for the transcendental part of the quasipolynomial. To do so, we finally introduce the symbol to indicate the vanishing of a given coefficient of the transcendental part of the quasipolynomial. To illustrate the above notions, consider the following quasipolynomial function Q(λ, a, τ ) = P 0 (λ, a) + (a 1,0,0,0 + a 1,0,0,1 λ) e -λτ1 +a 0,1,0,2 λ 2 e -λτ2 + a 0,0,1,1 λe -λτ3 .

In this case, σ i,j,k = -(i j k) • τ where i, j, k = 0, 1. More precisely, σ 1,0,0 = -τ 1 , etc. According to the above considerations, deg(P 0 ) = n ≥ 3, and the transcendental part of ( 18) is characterized by the incidence vector

V = (x 1 , x 1 , , , x 2 , , x 3 ).
Namely, the first two components of V indicate that P M 1 is a complete polynomial with deg(P M 1 ) = 1, the three components , , x 2 , indicate that a 0,1,0,0 = a 0,1,0,1 = 0 and P M 2 is lacunary with deg(P M 2 ) = 2 and , x 3 indicate that a 0,0,1,0 = 0 and P M 3 is lacunary with deg(P M 3 ) = 1.

In the spirit of the definition of functional confluent Vandermonde matrices introduced by [START_REF] Ha | A note on the determinant of a functional confluent Vandermonde matrix and controllability[END_REF], the functional Birkhoff matrices are defined as follows:

Definition 3.1. The square functional Birkhoff matrix Υ is associated to a sufficiently regular real function : R → R and an incidence matrix E (or equivalently an incidence vector V) and is defined by concatenation

Υ = [Υ 1 Υ 2 . . . Υ M ] ∈ M δ (R), where Υ i = [κ (ki 1 ) (x i ) κ (ki 2 ) (x i ) . . . κ (ki d i ) (x i )],
such that k i l ≥ 0 for all (i, l) ∈ {1, . . . , M } × {1, . . . , d i } and

M i=1 d i = δ with κ(x i ) = (x i )[1 . . . x δ-1 i ] T , for 1 ≤ i ≤ M.
Remark 3.2. In the sequel, if n = deg(P 0 ) let η = n + 1. Then, we consider (x i ) = x η i . By analogy to the Birkhoff interpolation problem, the non degeneracy of the matrices Υ will be a fundamental assumption for investigating the codimension of the zero spectral values for DDEs.

Remark 3.3. When η = 0, the matrix Υ is nothing else but the standard Birkhoff matrix and thus (x i ) = 1. If, in addition, V does not contain 's then we recover the confluent Vandermonde matrix [START_REF] Ha | A note on the determinant of a functional confluent Vandermonde matrix and controllability[END_REF]). The particular case d i = 1 for i = 1, . . . , M corresponds to the standard Vandermonde matrix and, in this case, M = δ since Υ is assumed to be a square matrix.

Thanks to this formalism based on functional Birkhoff matrices initially presented in Boussaada and Niculescu (2016a), one is able to characterize the codimension of imaginary spectral values of DDEs of retarded type.

Codimension of Zero Singularities of DDEs

The next result gives a bound for the zero spectral value, bound that takes into account the system structure. However, it should be noted that the proof is constructive and it exploits the existing links between the multiplicity of the zero singularity and Birkhoff matrices, see, for instance, Boussaada and Niculescu (2016a). Additionally, it gives the values of the system parameters guaranteeing an admissible multiplicity for the zero spectral value. Finally, the result holds even when the delay associated polynomials are sparse.

Proposition 3.4. The following assertions hold:

(i) The multiplicity of the zero root for the generic quasipolynomial function (17) cannot be larger than P S = D + ÑN,n , where D is the sum of degrees of the polynomials involved in the quasipolynomial and ÑN,n + 1 is the number of the associated polynomials. Moreover, such a bound is reached if, and only if, the parameters of (17) satisfy simultaneously:

a 0,k = - i∈S N,n a i,k + k-1 l=0 a i,l σ i k-l (k -l)! , 0 ≤ k ≤ P S -1. ( 19
)
(ii) Consider a quasipolynomial function (17) containing at least one incomplete polynomial for which we associate an incidence vector V Ẽ .

When the associated functional Birkhoff matrix Υ Ẽ is nonsingular, then the multiplicity of the zero root for the quasipolynomial function (17) cannot be larger than n plus the number of nonzero coefficients of the polynomial family (P M k ) M k ∈S N,n .

Remark 3.5. In the generic case, the Pólya-Szegö bound P S is completely recovered by the first assertion of Proposition 3.4. But, its advantage consists in providing the parameter values ensuring any admissible multiplicity of the zero singularity. The proof of Proposition 3.4 provides a constructive linear algebra alternative for identifying such a bound.

Remark 3.6. Obviously, the number of non-zero coefficients of a given quasi-polynomial function is bounded by its degree plus its number of polynomials. Thus, the bound proposed in Proposition 3.4 ii) is sharper than P S , even in the generic case, that is all the parameters of the quasipolynomial are left free, these two bounds are equal. Indeed, in the generic case, that is when the number of the left free parameters is maximal, the Pólya-Szegö bound P S = D + ÑN,n = n + D q + ÑN,n which is nothing else than n plus the number of parameters of the polynomial family (P M k ) M k ∈S N,n .

Remark 3.7. When the matrix Υ Ẽ is singular, one keeps the generic Pólya-Szegö bound P S .

Remark 3.8. The above proposition can be interpreted as follows. Under the hypothesis:

Q(iω) = 0 ⇒ ω = 0 (H)
(that is, all the imaginary roots are located at the origin), the dimension of the projected state on the center manifold associated with zero singularity for equation ( 17) is less or equal to its number of nonzero coefficients minus one. Indeed, under (H), the codimension of the zero spectral value is equal to the dimension of the state on the center manifold since, in general, the dimension of the state on the center manifold is nothing but the sum of the dimensions of the generalized eigenspaces associated with the spectral values having a zero real part.

Remark 3.9. For a given quasipolynomial function, the multiplicity of real roots may reach the Polya-Szegö bound, see Boussaada et al. (2014); Boussaada and Niculescu (2014). However, the multiplicity of an imaginary crossing root may exceed the dimension of the delay-free system but never reach the Polya-Szegö bound, see, e.g., Boussaada and Niculescu (2016b).

Multiple Induced Dominancy and Partial Pole Assignment: Comprehensive Examples

In this section, we focus on three case studies. The first one corresponds to the simplest DDE by which one can exhibit and explain in a comprehensive way the stability induced when forcing some particular spectral value to be multiple. The second, which apparently illustrates the limitation of such a property, gives an in-depth insight into the conditions of application of such a property. The last example concerns the stabilization of a second-order oscillator controlled by a delayed-output feedback. An appropriate choice of the gain and of the delay (seen as a parameter of the system) allows assigning a dominant characteristic root in C with multiplicity 3.

Exponential Decay Rate of a Scalar DDE With a Single Delay

Let us revisit the scalar DDE with one delay given in (6):

ξ(t) + a ξ(t) + b ξ(t -τ ) = 0, ( 20
)
where b > 0, . It was shown in Boussaada et al. (2016) that for a given positive delay, equation ( 20) admits a double spectral value at λ = λ 0 if, and only if,

λ 0 = - aτ + 1 τ and b = e λ0 τ τ . ( 21
)
In addition, it was emphasized that λ 0 is the rightmost root, and that if λ 0 < 0 then the zero solution of system ( 20) is asymptotically stable. Now, recall that λ = λ 0 is a spectral value of ( 20) if, and only if, λ 0 is a root of the characteristic function:

Q(λ, p, τ ) = λ + a + b e -λτ . ( 22
)
The main ingredient of the dominancy proof of λ 0 is an integral equation which cannot be satisfied for any spectral value λ with (λ) > λ 0 . Namely, it was shown that if b satisfies ( 21), then the characteristic function reads:

Q(λ, p, τ ) = (λ -λ 0 ) 1 - 1 0 e -τ (λ-λ0) t dt . ( 23
)
As a matter of fact, if λ 1 = ζ + iη = λ 0 is a root of ( 23) then λ 1 is a root of its second factor. Hence, we obtain

1 = 1 0 e -τ (ζ-λ0) t dt.
But, e -τ (ζ-λ0) t < 1 for ζ -λ 0 > 0 and 0 < t < 1, thereby exhibiting the dominancy of λ 0 .

Remark 4.1. The rightmost root λ 0 corresponding to equation ( 22), where system ( 21) is satisfied, varies in the interval λ 0 ∈ (-∞, -a). Figure 3 illustrates the behavior of the rightmost root with respect to the variation of τ .

Multiple Spectral Values for DDEs Systems Are Not Necessarily Dominant

The problem of stabilization of a chain of integrators is considered in [START_REF] Niculescu | Stabilizing a chain of integrators using multiple delays[END_REF] where a single integrator can be stabilized by a single delay state-feedback. Indeed, a positive gain guarantees the closed-loop stability of the system free of delay, and, by continuity, there exists a (sufficiently small) delay in the output preserving the stability of the closed-loop system. However, the situation is completely different for a chain of integrators of order n when n > 1. For instance, consider the time-delay system characterized by the following quasipolynomial function:

Q(λ, p, τ ) = λ 2 + α e -τ λ . ( 24
)
It can be checked that the maximal admissible multiplicity is 2 and it can be attained if, and only if,

α = -4 e -2 τ 2 , λ = - 2 τ . ( 25
)
However, the main result from [START_REF] Niculescu | Stabilizing a chain of integrators using multiple delays[END_REF] asserts that either n distinct delays or a proportional+delay compensator with n -1 distinct delays are sufficient to stabilize a chain including n integrators. In [START_REF] Kharitonov | Static output feedback stabilization: necessary conditions for multiple delay controllers[END_REF], a like assertion is shown to be also necessary to stabilize the chain of n integrators. Hence, in our case, either 2 distinct delays or a proportional+delay are necessary and sufficient to stabilize the double integrator. In conclusion, there exists at least a spectral value for (24) with a positive real part. As a result, λ 0 = -2 τ , while being a multiple root, it is not dominant. Indeed, consider ( 24)-( 25) with τ = 1, that is

Q(λ, p, τ )| τ =1 = λ 2 -4e -(λ+2) . ( 26
)
Figure 4. Sparsity-induced loss of dominancy for the multiple spectral value. Each intersection between the solid blue/dashed red curves corresponds to a spectral value of ( 26). For τ = 1, the dominancy of λ 1 ≈ 0.557 with respect to the double root at λ 0 = -2 is illustrated. The roots' distribution is illustrated using the QPmR toolbox from [START_REF] Vyhlídal | Mapping based algorithm for large-scale computation of quasi-polynomial zeros[END_REF]. The figure is borrowed from [START_REF] Boussaada | Multiplicity-induced-dominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF].

As illustrated in Figure 4, the dominancy property is lost since λ 1 ≈ 0.557 is a root of ( 26). This is justified by the sparsity (the vanishing of some coefficients of the corresponding P 0 polynomial) of ( 26).

Stabilizing an Oscillator via a Delayed Output-Feedback

Consider the stability of the trivial solution corresponding to the control system ξ(t) + c 1 ξ(t) + c 0 ξ(t) = u(t),

φ(t) = γ ξ(t),
where φ(t) is the system output and u(t) = β γ ξ(t -τ ) is the control law, which amounts to studying the roots of the quasipolynomial function

Q(λ, p, τ ) = λ 2 + a 1 λ + a 0 + β e -τ λ .
Using the standard linear change of variables λ → c1 λ 2 , one obtains the normalized characteristic function

Q(λ, p, τ ) = λ 2 + 2 λ + a 0 + α e -λτ (27)
where α = 4 c 2 1 β, τ = c1 2 τ and a 0 = 4 c0 c1 2 . If α = 0, the spectral abscissa is minimized at a 0 = 1 which corresponds to the rightmost root located at λ 0 = -1, see for instance [START_REF] Kirillov | Robust stability at the swallowtail singularity[END_REF]. By exploiting the delay effect, the following proposition proved in Boussaada et al. (2018) asserts that the solution's decay rate can be further improved by decreasing the corresponding rightmost root. Assume that a 0 > 1, then we have: Theorem 4.2 (Boussaada et al., 2018). The following properties hold:

i) The multiplicity of any given root of the quasipolynomial ( 27) is bounded by 3. ii) The quasipolynomial (27) admits a real spectral value at s 0 = -1 -1 τ with algebraic multiplicity 3 if, and only if,

τ = 1 a 0 -1 with α = - 2 e -(1+τ ) τ 2 . ( 28
)
iii) If equations (28) are satisfied then s = s 0 is the rightmost root of (27).

Remark 4.3. If equalities (28) are satisfied then the trivial solution of the second-order equation ξ(t) + 2 ξ(t) + a 0 ξ(t) = α ξ(t -τ ) is asymptotically stable with ξ(t) ≈ e s0t .

Parametric MID for Second-Order Systems

Second-order linear systems capture the dynamic behavior of many natural phenomena and have found numerous applications in a variety of fields, such as vibration and structural analysis. Stabilization of solutions to such a reduced order model represents a standard test bench to approve new paradigms and methodologies in control design.

The problem we consider in the sequel is the characterization of the admissible multiplicities of spectral values and their effect on stability of the DDE

ξ(t) + a 1 ξ(t) + a 0 ξ(t) + α 0 ξ(t -τ ) + α 1 ξ(t -τ ) = 0
It should be mentioned that such a system depicts several unexpected properties and in particular, for each choice of the coefficients, the delay parameter induces a stabilizing/destabilizing behavior. In other words, increasing the delay value may be beneficial for the system's dynamics 21 . ¿From a control theory viewpoint, the aim is to construct an appropriate delayed-state-feedback controller given by u(t) = -α 0 ξ(t -τ ) -α 1 ξ(t -τ ) able to stabilize the following dynamical system:

ξ(t) + a 1 ξ(t) + a 0 ξ(t) = u(t). (29)
In frequency-domain, the latter generic control problem yields the following characteristic quasipolynomial function:

Q(λ, p, τ ) = P 0 (λ) + P 1 (λ) e -τ λ (30)
where P 0 (λ) = λ 2 + a 1 λ + a 0 and P 1 (λ) = α 1 λ + α 0 . ¿From an algebraic geometry viewpoint, it is consistent to set up such a control problem via the discriminant of the characteristic polynomial related to the open-loop (uncontrolled) equation. Namely, the complex/real nature of the roots of the polynomial P 0 (λ) has a strong effect on the characterization of the controller's gains and delay enabling a desired fast stabilization using the MID property. In the sequel, it shall be emphasized that the discriminant of P 0 (∆ = a 2 1 -4a 0 ) defines an efficient and necessary criterion to exhibit the potential applicability and the limitations of the MID methodology, and our analysis splits following the sign of the discriminant ∆. The following results were presented in [START_REF] Boussaada | Multiplicity-induced-dominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF].

Open-Loop Systems With One Oscillating Mode

One oscillating mode as solution of the uncontrolled equation (two conjugate complex roots) corresponds, from a purely algebraic viewpoint, to a characteristic polynomial with a strictly negative discriminant ∆. In such a case, the following theorem gives a bound for quasipolynomial root's multiplicity and provides the explicit MID-based controller's gains and delay.

Theorem 5.1. Consider the quasipolynomial (30). Then, the following assertions hold: i) the multiplicity of any given root of the quasipolynomial (30) is bounded by 4, and this can only be achieved on the real axis and under negativity of ∆. ii) the quasipolynomial (30) admits a real spectral value at λ = λ 0 with algebraic multiplicity 4 if, and only if,

λ 0 = - a 1 + √ -2 ∆ 2 ,
and the system's parameters satisfy

τ = 2 - 2 ∆ , α 0 = 5 ∆ -a 1 √ -2 ∆ 4 e s0τ , α 1 = - √ -2 ∆ 2 e s0 τ . (31)
iii) if (31) is satisfied, then λ = λ 0 is the spectral abscissa corresponding to (30). iv) if (31) is satisfied then the trivial solution of the closed-loop system (29) is asymptotically stable if, and only if, either a 1 ≥ 0 and a 0 > a 2 1 4 or a 1 < 0 and a 0 > 3 a 2 1 8 .

Open-Loop Systems With Non Oscillating Modes

Consider now the case where the uncontrolled equation admits two real spectral values. As expected, under such a configuration, the MID baseddesign cannot be applied by exploiting the maximal multiplicity which is equal to 4, see Theorem 5.1. The following theorem gives a bound for the quasipolynomial roots' multiplicity in the case when ∆ ≥ 0. In addition, it explicitly provides the MID-based controller's gains and delay by exploiting a lower multiplicity.

Theorem 5.2. Consider the quasipolynomial (30). Then the following assertions hold: i) if the discriminant ∆ ≥ 0, then the multiplicity of any given root of the quasipolynomial (30) is bounded by 3. ii) for an arbitrary positive delay τ , the quasipolynomial (30) admits a real spectral value at λ = λ ± with algebraic multiplicity 3 if, and only if,

λ ± = -τ a 1 -4 ± √ 8 + τ 2 ∆ 2 τ ,
and the system parameters satisfy:

       α 0 = a 1 λ ± + a 1 2 2 - ∆ 2 + 6 a 1 + 10 λ ± τ + 6 τ 2 e λ±τ , α 1 = 2 λ 0 + a 1 + 2 τ e λ±τ .
iii) if ( + ) (respectively ( -)) is satisfied, then λ = λ + is the spectral abscissa corresponding to (30) (respectively λ -cannot be the spectral abscissa corresponding to (30)). Furthermore, for an arbitrary delay τ the multiple spectral value at λ -is always dominated by a single real root λ 0 . iv) if ( + ) is satisfied then the trivial solution is asymptotically stable if, and only if, τ satisfies the following conditions

     τ ∈ (0, τ -) when a 0 < 0, or τ ∈ (0, τ -) ∪ (τ + , ∞) when a 0 > 0 and a 1 < 0,
where τ ± = -a1± √ ∆+2 a0 a0 .

Remark 5.3. The second and third assertions of Theorem 5.2 hold also for negative discriminant -8 τ 2 < ∆ < 0. Moreover, when ∆ = -8 τ 2 the triple root at λ ± becomes the quadruple root prescribed in Theorem 5.1.

6 The Generic MID Property Some works have shown that, for some classes of dynamical systems represented by DDEs of retarded type, a real root of maximal multiplicity is necessarily the rightmost root. Such a a property is called Generic MID, or GMID for short. This link between maximal multiplicity and dominance has been suggested in [START_REF] Pinney | Ordinary difference-differential equations[END_REF] after the study of some simple, loworder cases, but without any attempt to address the general case. Recently, a more general result on the GMID property from Mazanti et al. (2021a) for generic retarded DDEs of order n with delayed (polynomial) "term" of order n -1, which relies on links between quasipolynomials with a real root of maximal multiplicity and the Kummer's confluent hypergeometric function in terms of the location of the characteristic roots. The GMID property was also extended to neutral DDEs of orders 1 and 2 in [START_REF] Ma | Multiplicity-Induced-Dominancy extended to neutral delay equations: Towards a systematic PID tuning based on Rightmost root assignment[END_REF]; [START_REF] Benarab | The MID property for a second-order neutral time-delay differential equation[END_REF]; Mazanti et al. (2021b), as well as to the case of complex conjugate roots of maximal multiplicity in [START_REF] Mazanti | Spectral dominance of complex roots for single-delay linear equations[END_REF].

Degenerate Hypergeometric Functions

M (a, b, λ) = ∞ k=0 (a) k (b) k λ k k! . ( 32
)
Furthermore, the series in (32) converges for every λ ∈ C. As presented in [START_REF] Buchholz | The confluent hypergeometric function with special emphasis on its applications[END_REF]; [START_REF] Erdélyi | Higher transcendental functions[END_REF]Olver et al. (2010), the function M (a, b, •) satisfies Kummer's differential equation

λ ∂ 2 M ∂λ 2 (a, b, λ) + (b -λ) ∂M ∂λ (a, b, λ) -aM (a, b, λ) = 0. ( 33
)
Other solutions of (33) are usually also called Kummer's confluent hypergeometric functions, but they are not used in this chapter. Note that the Kummer's confluent hypergeometric function M admits an integral representation, which can be found, for instance, in [START_REF] Buchholz | The confluent hypergeometric function with special emphasis on its applications[END_REF] The main result on confluent hypergeometric functions used in Mazanti et al. (2021a) to prove the GMID property is the following one on the location of the roots of some particular functions. It was proved in [START_REF] Wynn | On the zeros of certain confluent hypergeometric functions[END_REF] by using a continued fraction expansion of the ratio of two such functions.

Proposition 6.1. Let a ∈ R be such that a > -1 2 . 1. If λ ∈ C is such that M (a, 2a + 1, λ) = 0, then (λ) > 0. 2. If λ ∈ C is such that M (a + 1, 2a + 1, λ) = 0, then (λ) < 0.

Spectral Values of Maximal Multiplicity Are Dominant

Consider the quasipolynomial function

Q : C → C defined for λ ∈ C by Q(λ, τ ) = λ n + n-1 k=0 a k λ k + e -λτ n-1 k=0 α k λ k . ( 34
)
The following result proved in Mazanti et al. (2021a) is the following characterization of real roots of maximal multiplicity of Q, their dominance and the corresponding consequences for the stability of the trivial solution of

y (n) (t) + n-1 k=0 a k y (k) (t) + n-1 k=0 α k y (k) (t -τ ) = 0, ( 35
)
Theorem 6.2 (Mazanti et al., 2021a). Consider (34) and let λ 0 ∈ R.

1. The number λ 0 is a root of multiplicity 2n of Q if and only if, for every k = 0, . . . , N -1,

             a k = n k (-λ 0 ) n-k + (-1) n-k n! n-1 j=k j k 2n -j -1 n -1 λ j-k 0 j!τ n-j , α k = (-1) n-1 e λ0τ n-1 j=k (-1) j-k (2n -j -1)! k!(j -k)!(n -j -1)! λ j-k 0 τ n-j .
(36) 2. If (36) is satisfied, then λ 0 is a strictly dominant root of Q.

If

ζ(t) + 2 ξ ω 0 ζ(t) + ω 2 0 ζ(t) = c(t)
where ω 0 > 0 and 0 < ξ < 1 stand respectively for the oscillator natural frequency and the damping factor. Consider now the controller c as a proportional-derivative delayed-controller:

c(t) = -b 0 ζ(t -τ ) -b 1 ζ(t -τ ).
Thus, the closed-loop characteristic function is given by:

Q(λ, p, τ ) = λ 2 + 2 ξ ω 0 λ + ω 0 2 + (b 0 + b 1 λ) e -τ λ .
Assume the natural frequency ω 0 = 1 and the damping factor ξ = 1/2, which corresponds to an open-loop plant with a complex-conjugate pair λ ± OL = -1/2 ± i √ 3/2. Then, the closed-loop plant corresponds to the following characteristic quasipolynomial

Q(λ, p, τ ) = λ 2 + λ + 1 + (b 0 + b 1 λ) e -τ λ . ( 37
)
Forcing the existence of a triple spectral value suggests that

λ ± CL = -1/2 -2 τ -1 ± 1/2 τ -1 -3 τ 2 + 8 (38)
are the only admissible roots. As a matter of fact, those triple spectral values are defined if, and only if, the controller's gains are such that:

b 0 = 6 + 2 + λ ± CL τ 2 + 10 λ ± CL + 6 τ e λ ± CL τ τ 2 , b 1 = e λ ± CL τ 2 λ ± CL τ + τ + 2 τ . (39)
It follows that if (39) is satisfied and the triple root at λ + CL is selected then the spectral abscissa corresponds to λ + CL as illustrated in Figure 5. Notice that the assignment of the triple root λ = λ + CL is possible only in 3), see for instance [START_REF] Boussaada | Multiplicity-induced-dominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF]. This fact is illustrated in Figure 6.

Active Vibration Control in a Mechanical Flexible Structure

The problem of active vibration damping of thin mechanical structures is a topic that has received a great attention by the control community since several years (see, for instance, Tliba, 2011 and the references therein), espe-Figure 6. MID stabilizability region defined by (-∞, λ M ], in which the assignment of a triple negative dominant root of (37) is possible, where in such a case λ M corresponds to the lowest (in absolute value) assignable decay rate. cially when actuators and sensors are based on piezoelectric materials. For deformable mechanical structures, piezoelectric materials are used as strain sensors or strain actuators. Many works have concerned the vibration control problem of the "Euler-Bernoulli beam" equipped with one rectangular piezoelectric actuator and sometime, another one, identical and collocated, but used as sensor. See for example [START_REF] Chen | Optimal sensor design and control of piezoelectric laminate beams[END_REF]; [START_REF] Banks | Proper orthogonal decomposition-based control of transverse beam vibrations: experimental implementation[END_REF] where one edge of the beam is clamped whereas the other remains free. Other works dealt with the problem of vibration control for laminated rectangular plates [START_REF] Kögl | Analysis of smart laminates using piezoelectric mitc plate and shell elements[END_REF] or complex plate like structures [START_REF] Tliba | Active vibration damping of a smart flexible structure using piezoelectric transducers: H ∞ design and experimental results[END_REF].

Modeling of the Vibrating Beam

Consider now the flexible structure depicted in Fig. 7. It is an aluminiumbased beam, embedded in a mobile support. The mobile support is subjected to an acceleration, denoted by w in the sequel, and it is moving along the z axis. This flexible beam is equipped with two piezoelectric patches made with lead zirconate titanate (also called PZT). One of them is used as an actuator and the other works as a sensor. These patches are supposed to be rigidly bounded on the beam, one on each side, located at the clamped edge. The whole device is called thereafter as a piezo-actuated beam. It can be deformed by the application of a voltage, denoted by u, across the actuator. The sensor delivers an electrical voltage which corresponds to a measure, denoted by y, of the local deformation under the piezoelectric patch. Very often, this equipped mechanical structure is partly described by the in-plane Euler-Bernoulli Partial Differential Equation (PDE) that suffers from the lack of precision in describing the electro-mechanical interactions between the passive structure and the piezoelectric components. Indeed, these latter are often withdrawn in the computation of the eigenfrequencies [START_REF] Tliba | Contrôle actif des vibrations dans des structures mécaniques minces instrumentés de transducteurs piézoélectriques[END_REF] of the whole structure. Nevertheless, such a structure obeys to fundamental equations of continuum mechanics in 3D space [START_REF] Gérardin | Mechanical vibrations: theory and application to structural dynamics[END_REF], involving computations of gradient of displacement vector and divergence operator applied to strength tensor. When completed with Neumann and Dirichlet boundary conditions, the fundamental equations give several PDEs that are coupled, thus hardly or impossible to solve analytically. Then, for controller design purposes, one naturally turns toward numerical methods in order to get the input-to-output dynamical models [START_REF] Komzsik | What every engineer should know about computational techniques of finite element analysis[END_REF].

Using Finite Element Modeling (FEM), one obtains the LTI system described by the following ordinary differential equations:

M qq q(t) + D qq q(t) + K qq q(t) = M qw w(t) -K qu u(t) (40)
y(t) = K qy q(t) (41) z(t) = F zw w(t) -F zu u(t) -F zq q(t) -F zv q(t) (42)
where w(t) ∈ R is the absolute acceleration (m/s 2 ) of the movable support along axis z, z(t) ∈ R is the relative acceleration (m/s 2 ) of the free end, derived from the equations of motion, u(t) ∈ R is the piezoelectric voltage (V ) across the actuator (control signal), y(t) ∈ R is the piezoelectric sensor voltage (V ) across the sensor (measured output signal y(t)). Moreover, the terms M qq , D qq , M qw , K qq , K qu , K qy , F zw , F zu , F zq and F zv are all matrices derived from the assembly step of the FEM such that the obtained approximation is with several thousands degrees of freedom.

After producing a FEM, a modal analysis is performed to the undamped motion equation ( 40). It consists in finding the eigenstructure of M qq q(t) + K qq q(t) = 0, and using the eigenvalues and eigenvectors as a new basis allowing the description of ( 40)-( 42) into a new system of ODEs, still linear but decoupled, involving a new state vector, called vector of modes. The advantage of this form is that it allows to build a model in statespace form devoted to analysis, (43) below, having its order sufficiently small to describe the dynamical behavior within a low-frequency bandwidth. Furthermore, a reduced-order system, devoted to the controller synthesis step, is also available and differs from the analysis one by the presence of direct feedthrough terms between the outputs {z, y} and the inputs {w, u}. Further details can be found in Boussaada et al. (2018). The frequency responses for the analysis and the reduced-order models are shown in Fig. 8 and the shapes of the first three bending modes in Fig. 9. . By applying the MID property and the ensuing procedure, one can assign λ 0 as a rightmost root of multiplicity equal to 3 such that λ 0 -192.359.

Notice that in such a case the numerical values for the parameters of the controller in (44) are given by: (blue) and in closed-loop (red) when the disturbance w is a rectangular impulse (black), say like a shock. We also put the time response of the control signal u that exhibits a peak of magnitude roughly equal to -60 V which is reasonable for this application.

Notes and Comments

In studying the exponential stability the trivial solution of linear DDEs of retarded type, this chapter addressed several frequency-domain techniques to count the roots of the corresponding characteristic functions in the complex right half-plane with a particular emphasis on the MID property and its potential application in control engineering. Illustrative examples and one application in vibration control show the effectiveness of the method. The main results of this chapter have been reported in Boussaada and Niculescu (2016a)

Figure 2 .

 2 Figure2. Solid blue curve: simplified contour used inBoussaada and Niculescu (2018) for applying the argument principle to investigate the dominancy of multiple roots; dashed red curve: generic spectrum envelope given by Proposition 2.1 for the quasipolynomial Q(λ, p, τ ) = λ + a 0 + a 1 e -λτ under conditions λ 0 = -a0τ +1

Figure 3 .

 3 Figure 3. (Left) Distribution of the spectrum corresponding to equation (22) and conditions (21) for a = τ = 1 and b = e -2 . The roots' distribution is illustrated using the QPmR toolbox from Vyhlídal and Zitek (2009). (Right) Rightmost root corresponding to equation (22) and conditions (21) as a function of the delay τ for a = 1 and b = e -(τ +1) τ . Figures are borrowed from Boussaada et al. (2020).

  Let a, b ∈ C and assume that b is not a nonpositive integer. The Kummer's confluent hypergeometric function M (a, b, •) : C → C is the entire function defined for λ ∈ C by the series

0 e

 0 ;[START_REF] Erdélyi | Higher transcendental functions[END_REF] Olver et al. (2010). Namely, let a, b ∈ C and assume that (b) > (a) > 0. ThenM (a, b, λ) = Γ(b) Γ(a)Γ(b -a) 1 λt t a-1 (1 -t) b-a-1 dt,where Γ denotes the Gamma function.

Figure 5 .

 5 Figure 5. P3δ interface exhibiting the design of a stabilizing delayed PD controller in the case of the characteristic function (37). (Left) Illustration of the roots location in the case τ = τ * ≈ 0.24 which by using (38) corresponds to λ + CL = -3. (Right) Closed-loop response corresponding to the history function ϕ(t) = 3 for all t ∈ [-τ * , 0].

Figure 7 .

 7 Figure 7. Sketch of the piezo-actuated flexible beam, clamped at one edge. The figure is taken from Boussaada et al. (2017).

Figure 8 .Figure 9 .

 89 Figure 8. Frequency responses of the analysis (-) and reduced order (.-) models. The figure is borrowed from Boussaada et al. (2017).

Figure 10 .

 10 Figure 10. Feedback control structure. The figure is borrowed from Boussaada et al. (2017).

nFigure 11 .

 11 Figure 11. Time responses of the measured output y (Top) and of the controlled output z (Bottom).

  (36) is satisfied, then the trivial solution of (35) is exponentially stable if and only if a n-1 > -n 2 τ .Based on the recent results by the authors on the MID property for DDEs, a Python software for the parametric design of stabilizing feedback laws with time-delays, called "Partial Pole Placement via Delay Action" (P3δ 22 for short), has been developed. P3δ also implements other features, which are detailed in Boussaada et al.. In the sequel, we revisit the problem of damping solutions of the standard oscillator by using a delayed feedback, and the proposed results are illustrated by using P3δ:

	7 Software: Partial Pole Placement via Delay Action

  (characterizing the codimension of the zero singularity) Mazanti et al. (2021a) (characterizing the generic MID characterization for DDEs of retarded type with a single delay) and Boussaada et al. (2017) (applica-tion in vibration control). It should be noted that a guided tour of existing methods to analyze multiple characteristic roots (including the MID property mentioned in this chapter) can be found in Niculescu et al. (2021).

More precisely, the estimation of the upper bound of the number of the characteristic roots of exponential polynomials located on a horizontal strip.

Based on the same Cauchy's argument principle.

Whose origins date back to the work of Evans in the late 1950s,[START_REF] Evans | Control system synthesis by root locus method[END_REF];Mac- Farlane (1979).

Such links have been previously highlighted inMazanti et al. (2021a) for a particular class of systems and then generalized in[START_REF] Boussaada | The generic multiplicityinduced-dominancy property from retarded to neutral delay-differential equations: When delay-systems characteristics meet the zeros of Kummer functions[END_REF].

It should be mentioned that, for reading facility, some specific notations are introduced later on.

Belonging to the Banach space of continuous functions C([-τ, 0], R n ).

In terms of existence of solutions for the corresponding Cauchy problem.

that is the set of eigenvalues of the corresponding characteristic matrix function

depending on the parameters -→ p and τ

It should be noted that, even in the simple case of second-order linear DDEs of retarded type, the corresponding 5 parameters induce a large computational effort to perform such an analysis.

The reader may find the proof of this theorem in any standard textbook on the theory of analytic functions.

In our case, we deal with quasipolynomial functions.

Quasipolynomials are analytic functions.

Called also the argument variation.

Which is strongly inspired from Stépán's result.

More precisely, Pólya and Szegő (1972, Problem 206.2, page 144 and page 347).

Namely, the well-known Birkhoff interpolation problem.

Bogdanov-Takens singularity corresponds to equilibrium point admitting a double spectral value at zero for which corresponds a single eigenvector. Near a Bogdanov-Takens singularity a codimension 2 bifurcation may occur in two dimensional parameters space

.19 Indeed, in this case, we have N = 2, τ 1 = 0, τ 2 = τ and d 1 = 1, d 2 = 0.

In other words, the sum of the degrees of the polynomials involved in the quasipolynomial plus the number of polynomials involved minus one is called the degree of a given quasi-polynomial. Further discussions on such a notion can be found in[START_REF] Wielonsky | A Rolle's theorem for real exponential polynomials in the complex domain[END_REF].

Such a property does not hold for scalar systems with one delay, for which increasing the delay leads to instability, see for instance[START_REF] Walton | Direct method for TDS stability analysis[END_REF];[START_REF] Niculescu | Delay Effects on Output Feedback Control of Dynamical Systems[END_REF].

The software is freely available for download on https://cutt.ly/p3delta, where installation instructions, video demonstrations, and the user guide are also available.

Acknowledgements

This work is supported by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissement d'Avenir" program, through the iCODE project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. The authors warmly thank their collaborators Tamas Insperger (BME, Budapest), Guilherme Mazanti (L2S, University Paris-Saclay) and Sami Tliba (L2S, University Paris-Saclay).

mode, one gets the following numerical data for the reduced order model N wz (λ) = -1.572 λ 2 -0.767 λ + 0.114, N uz (λ) = 0.040 λ 2 + 0.019 λ -0.002, N wy (λ) = -0.047 λ 2 -0.023 λ -24664.720, ψ(λ) = λ 2 + 0.487 λ + 59495.866, N uy (λ) = 0.082 λ 2 + 0.040 λ + 5472.410.

An interesting control objective would be to damp the peak of resonance of the first bending mode, by using an output feedback controller, without affecting the vibrating modes that were neglected in the reduced order model. The result presented in the next section as well as its proof can be found 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], see also Boussaada et al. (2018).

Vibration Damping

Consider now the piezo-actuated system inserted in the typical output feedback control structure of Fig. 10, with a zero-reference signal and an input disturbance w corresponding to a rectangular impulse signal. The control problem consists in damping the vibrations due to the first mode when the mobile support is subjected to a shock-like disturbance. We define the output feedback control law u(λ) = C(λ, p, τ ) y(λ) involving the following reduced-complexity controller given in Laplace domain by

where N (λ, p, τ ) = n 0 + n r0 e -τ λ and D(λ, p, τ ) = d 0 + d r0 e -τ λ and p = (n 0 , n r0 , d 0 , d r0 ). By applying the inverse Laplace transform, it can be easily shown that this control law writes in time-domain as:

which is an output feedback control law based on proportional actions plus delayed proportional actions. For the sake of clarity, by omitting the variable λ in the polynomials, the closed loop relation between the disturbance w and the controlled output z is given by:

It can be proven that ψ divide the polynomial (N uz N wy -N wz N uy ) so that it can be removed from the numerator and the denominator. In this last