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

Guilherme Mazanti
email: guilherme.mazanti@l2s.centralesupelec.fr

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

Padé Approximation and Hypergeometric Functions: A Missing Link with the Spectrum of Delay-Differential Equations

Keywords: Padé approximation, delay-differential equations, exponential stability, Kummer hypergeometric functions

It is well known that rational approximation theory involves degenerate hypergeometric functions and, in particular, the Padé approximation of the exponential function is closely related to Kummer hypergeometric functions. Recently, in the context of the study of the exponential stability of the trivial solution of delay-differential equations, a new link between the degenerate hypergeometric function and the zeros distribution of the characteristic function associated with linear delay-differential equations was emphasized. Such a link allowed the characterization of a property of time-delay systems known as multiplicity-induced-dominancy (MID), which opened a new direction in designing low-complexity controllers for time-delay systems by using a partial pole placement idea. Thanks to their relations to hypergeometric functions, we explore in this paper links between the spectrum of delay-differential equations and Padé approximations of the exponential function. This note exploits and further comments recent results from [I. Boussaada, G.

INTRODUCTION

On beyond of the arithmetic theory and the number theory, the theory of continued fractions is strongly related to the theory of analytic functions; it is involved in topics such as the study of definite integrals, power series and the summation of divergent series, see for instance [START_REF] Perron | Die Lehre von den Kettenbrüchen[END_REF]; [START_REF] Wall | Analytic theory of continued fractions[END_REF]. Continuous fractions theory goes back to the pioneering work of T. S. Stieltjes mainly elaborated at the end of the 19th century. It appears that several links exist between continuous fractions theory with concepts and methods such as the well known Padé rational approximation and the ϵ-algorithm. As a matter of fact, the ϵ-algorithm is a transformation of the partial sums of a series into the corresponding Padé quotients or the corresponding continued fractions. In fact, the ϵ-algorithm is an old but still relevant way to improve the convergence rate of slowly convergent sequences, see, for instance, [START_REF] Graves-Morris | The epsilon algorithm and related topics[END_REF] and references therein. Moreover, a strong connection between such an algorithm and the Padé approximation has been underlined by P. Wynn in the early 60s, [START_REF] Wynn | On the convergence and stability of the epsilon algorithm[END_REF]. Through his results on continuous fractions theory, Wynn's work covered several problems in analytic functions theory, such as the distribution of zeros of a class of special functions called the Kummer (de-generate) hypergeometric functions or, equivalently, the Wittaker hypergeometric functions, see, e.g., [START_REF] Wynn | On the zeros of certain confluent hypergeometric functions[END_REF].

More contemporary contributions have pointed out to a missing link existing between the zeros of such a class of special functions with the spectrum distribution of an appropriate Delay-Differential Equation (DDE) including one pointwise delay, see for instance Mazanti et al. (2021a); Boussaada et al. (2022a). Interestingly, this link allowed deriving an exhaustive characterization of a property for DDE 1 called Generic Multiplicity Induced Dominancy (GMID), see [START_REF] Boussaada | The generic multiplicity-induced-dominancy property from retarded to neutral delay-differential equations: when delay-systems characteristics meet the zeros of Kummer functions[END_REF]. Indeed, the GMID is a special case of a more general property called Multiplicity Induced Dominancy (MID). Such a property consists of determining the conditions under which a given multiple complex zero of a quasipolynomial is dominant. It should be noted that the GMID property asserts that a DDE spectral value admitting the maximal multiplicity is necessarily the corresponding spectral abscissa. However, multiple roots with intermediate admissible multiplicities may or may not be dominant. Thanks to this property, a consistent control strategy is proposed in [START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: analysis and application in control design[END_REF]; [START_REF] Balogh | Conditions for stabilizability of time-delay systems with real-rooted plant[END_REF], which consists in assigning a root with an admissible multiplicity once appropriate conditions guaranteeing its dominancy have been established. Furthermore, the MID property may be used to tune standard controllers. For instance, in [START_REF] Ma | PID control design for firstorder delay systems via MID pole placement: Performance vs. robustness[END_REF] it is applied to the systematic tuning of the stabilizing PID controller of a first-order plant. Here, we aim at assigning dominant multiple real roots with admissible codimensions.

Beyond the reminder of the concepts and links discussed, the contribution of the present note is threefold. First, to highlight the existing links between roots of generic maximal multiplicity for a quasipolynomial, Kummer confluent geometric functions, and Padé approximations of the exponential function. Second, the use of Padé approximation allows simplifying some newly obtained analytic functions results, thus opening some new perspectives in a partial pole placement methodology for linear timeinvariant DDE. Finally, it allows to slightly correct some of the results derived in [START_REF] Saff | On the zeros and poles of Padé approximants to e z[END_REF] where the links between Kummer functions and Padé approximations are underlined.

The remaining of the paper is organized as follows: Section 2 recalls the definitions and properties of some concepts studied in the paper. Next, Section 3 is dedicated to the main result. Section 4 presents further comments on the distribution of the Kummer and the Whittaker zeros distribution. Some concluding remarks end the paper.

C + = {s ∈ C | ℜ(s) > 0}. 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 .

PREREQUISITES

This section provides a brief presentation of the definitions and results that shall be of use in the sequel concerning Padé rational approximation and degenerate hypergeometric functions.

On the Padé rational approximation

The Padé rational approximation can be seen as a generalization of the Taylor expansion of a function in which one wishes to approximate some function f around a given point (which is typically normalized to 0) by a rational function, i.e., by a function R expressed as R = Q P , where Q and P are polynomials. We start this subsection by a precise definition of what is meant by Padé approximation in this paper. Definition 2.1. Let f : U → C be an analytic function defined on a neighborhood U ⊂ C of 0. Given (m, n) ∈ N 2 , we say that a pair of polynomials (Q, P ) with deg Q ≤ m and deg

P ≤ n is a Padé approximation of f of order (m, n) if P (z)f (z) -Q(z) = O(z n+m+1 ) as z → 0.
Classical results on Padé approximations can be found, for instance, in [START_REF] Baker | Padé approximants, volume 59 of Encyclopedia of Mathematics and its Applications[END_REF] and (Brezinski, 2002, Chapter 4), and we also refer to [START_REF] Brezinski | History of continued fractions and Padé approximants[END_REF] for a historical presentation of this topic. For every (m, n) ∈ N 2 , the Padé approximation (Q, P ) of an analytic function f defined on a neighborhood of 0 exists and is unique up to a multiplicative constant. If n = 0, the corresponding Padé approximation (Q, 1) is such that Q coincides with the Taylor-Maclaurin expansion of f of order m.

Padé approximations are widely used in rational approximation theory (see, for instance, [START_REF] Szegő | Über eine Eigenschaft der Exponentialreihe[END_REF]; [START_REF] Dieudonné | Sur les zéros des polynomes-sections de e x[END_REF]; [START_REF] Saff | On the zeros and poles of Padé approximants to e z[END_REF]; [START_REF] Baratchart | Rational interpolation of the exponential function[END_REF] and the references therein). In particular, several applications involve Padé approximations of the exponential function z → e z and some of the above works highlighted links between such approximations and degenerate hypergeometric functions. As reported in (Perron, 1977, p246), the Padé approximation of order (m, n) of the exponential function is the pair

(Q n,m , P n,m ) defined for z ∈ C by            P n,m (z) = n k=0 (n + m -k)!n!z k k!m!(n -k)! , Q n,m (z) = m k=0 (n + m -k)!(-z) k k!(m -k)! , (1)
and the rational approximation

R n,m (z) = Qn,m(z) Pn,m(z) satis- fies R n,m (z) -----→ n,m→∞ e z ,
uniformly on compact subsets of C. Further, as can be found in [START_REF] Perron | Die Lehre von den Kettenbrüchen[END_REF], the Padé remainder e z -R n,m (z) satisfies

P n,m (z) (e z -R n,m (z)) = z n+m+1 m! 1 0 e tz (1 -t) m t n dt,
(2) which, as detailed in Section 2.2 below, is closely related to Kummer hypergeometric functions (see for instance [START_REF] Saff | On the zeros and poles of Padé approximants to e z[END_REF] and the references therein). We also highlight the fact that, from a control theory point of view, Padé approximations of the exponential function have been exploited in approximating delay systems via finite-dimensional systems (see, for instance, [START_REF] Glader | Approximation of delay systems-a case study[END_REF] for a deeper discussion on the approximation of the parametric family of first-order stable delay systems). Remark 2.2. The expressions of P n,m and Q n,m provided above use a different normalization with respect to the expressions in [START_REF] Perron | Die Lehre von den Kettenbrüchen[END_REF]: we have chosen here the normalization consisting in requiring the denominator P n,m to be monic, i.e., the coefficient of the term z n is 1.

Degenerate hypergeometric functions

We present in this section the definition and main properties of Kummer and Whittaker hypergeometric functions. The presentation in this section follows that of Boussaada et al. (2022a), and we start by providing the definition and main properties of Kummer hypergeometric functions (for further details, see, for instance, [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) and references therein). For an historical perspective, the reader is referred to Kampé de Fériet (1937). Definition 2.3. Let a, b ∈ C and assume that b is not a nonpositive integer. The Kummer confluent hypergeometric function Φ(a, b,

•) : C → C is the entire function defined for z ∈ C by the series Φ(a, b, z) = ∞ k=0 (a) k (b) k z k k! . (3)
The series in (3) converges for every z ∈ C, and the function Φ(a, b, •) satisfies the Kummer differential equation

z ∂ 2 Φ ∂z 2 (a, b, z) + (b -z) ∂Φ ∂z (a, b, z) -aΦ(a, b, z) = 0. (4)
The equation ( 4) admits two linearly independent solutions, which sometimes are both called Kummer confluent hypergeometric functions. In the present paper, we are concerned only with the solution given by (3).

Kummer functions admit the integral representation

Φ(a, b, z) = Γ(b) Γ(a)Γ(b -a) 1 0 e zt t a-1 (1 -t) b-a-1 dt, (5)
for every a, b, z ∈ C such that ℜ(b) > ℜ(a) > 0, where Γ denotes the Gamma function. For all complex numbers a, b, z such that b is not a nonpositive integer, Kummer functions also satisfy the relation

Φ(a, b, z) = e z Φ(b -a, b, -z), (6)
M k,l (z) = e -z 2 z 1 2 +l Φ( 1 2 + l -k, 1 + 2l, z). (8)
If 1 2 +l is not an integer, the function M k,l is a multi-valued complex function with branch point at z = 0. Whenever 2l is not a negative integer, the nontrival roots of M k,l coincide with those of Φ( 1 2 +l-k, 1+2l, •) and M k,l satisfies the Whittaker differential equation

φ ′′ (z) = 1 4 - k z + l 2 -1 4 z 2 φ(z). (9)
Similarly to the Kummer differential equation ( 4), other solutions of the Whittaker differential equation ( 9) are also known as Whittaker functions in other works, but they will not be used in this paper. Notice also that, since M k,l is a nontrivial solution of the second-order linear differential equation ( 9), any nontrivial root of M k,l is necessarily simple. The roots of Whittaker functions satisfy the following immediate symmetry property. Proposition 2.5. Let k, l ∈ C and assume that 2l is not a negative integer.

If z ∈ C \ {0} is a nontrivial root of M k,l , then -z is a root of M -k,l .
The following result from Boussaada et al. (2022a), whose proof relies on the Green-Hille transform of (9) (see, e.g., [START_REF] Hille | Oscillation theorems in the complex domain[END_REF]), provides information on the location of nontrivial zeros of Whittaker functions. Proposition 2.6. Let k, l ∈ R be such that 2 l -1 ≥ 0.

(a) If k = 0, then all nontrivial roots z of M k,l are purely imaginary.

(b) If k > 0, then all nontrivial roots z of M k,l satisfy ℜ(z) > 0. (c) If k < 0, then all nontrivial roots z of M k,l satisfy ℜ(z) < 0. (d) If k ̸ = 0, then all nontrivial roots z of M k,l satisfy 4k 2 ℑ(z) 2 -4(l 2 -k 2 ) -1 ℜ(z) 2 > 0. ( 10
)
Moreover, in all cases, all non-real roots

z of M k,l satisfy |z| > √ 4 l 2 -1.
As an immediate consequence of (8) and Proposition 2.6, we have the following result on the location of zeros of Kummer functions with real parameters, which was stated in Boussaada et al. (2022a).

Corollary 2.7. Let a, b ∈ R be such that b ≥ 2. (a) If b = 2a, then all nontrivial roots z of Φ(a, b, •) are purely imaginary. (b) If b > 2a, then all nontrivial roots z of Φ(a, b, •) satisfy ℜ(z) > 0. (c) If b < 2a, then all nontrivial roots z of Φ(a, b, •) satisfy ℜ(z) < 0. (d) If b ̸ = 2a, then all nontrivial roots z of Φ(a, b, •) satisfy (b -2a) 2 ℑ(z) 2 -(4a(b -a) -2b) ℜ(z) 2 > 0.
Moreover, in all cases, all non-real roots z of Φ(a, b, •) satisfy |z| > b(b -2).

Delay-differential equations and their characteristic equations

Consider the linear time-invariant delay-differential equation

y (n) (t) + n-1 k=0 a k y (k) (t) + m k=0 α k y (k) (t -τ ) = 0, (11)
where y(•) is the real-valued unknown function, τ > 0 is the delay, and a 0 , . . . , a n-1 , α 0 , . . . , α m are real coefficients. Such an equation is said to be of retarded type if m < n (i.e., if the highest order of derivation appears only in the non-delayed term y (n) (t)), or of neutral type if m = n. We refer to [START_REF] Hale | Introduction to functional differential equations[END_REF]; [START_REF] Michiels | Stability, control, and computation for time-delay systems: An eigenvaluebased approach[END_REF] for more information on delay-differential equations and for proofs of the major facts stated in the sequel.

Characterizing the stability of the zero solution of (11) in terms of the coefficients a 0 , . . . , a n-1 , α 0 , . . . , α m is in general a challenging question. Exponential stability can be characterized through the roots of the characteristic function of (11), which is the entire function ∆ :

C → C defined for s ∈ C by ∆(s) = s n + n-1 k=0 a k s k + e -sτ m k=0 α k s k . ( 12
)
More precisely, the zero solution of ( 11) is exponentially stable if and only if there exists γ < 0 such that, for every root s of ∆, we have ℜ(s) ≤ γ. Rightmost roots of ∆ thus determine the asymptotic behavior of the system and its rate of exponential stability.

It should be mentioned that the function ∆ from ( 12) is a particular case of a quasipolynomial. Quasipolynomials have been extensively studied in the literature (see, e.g., [START_REF] Bellman | Differentialdifference equations[END_REF], [START_REF] Hale | Introduction to functional differential equations[END_REF]) and, in particular, a result presented in (Pólya and Szegő, 1998, Part Three, Problem 206.2) implies that any root s of ∆ has multiplicity at most D P S = n + m + 1.

The number D P S is known as the degree of ∆, and such a bound on the multiplicity of a root of ∆ is known as Pólya-Szegő bound.

MAIN RESULTS

The main result of this paper is the following, which highlights the links between roots of maximal multiplicities of ∆, Kummer confluent hypergeometric functions, and Padé approximations of the exponential function. More precisely, our main result characterizes the case where ∆ admits a root of maximal multiplicity D P S through a suitable Kummer hypergeometric function by using the results on the Padé approximations of the exponential function recalled in Section 2.1. Theorem 3.1. Consider the quasipolynomial ∆ given by ( 12) and let s 0 ∈ R. Let P and Q be the polynomials defined by

P (s) = s n + n-1 k=0 a k s k , Q(s) = m k=0 α k s k
and define the quasipolynomial ∆ and the polynomials P and Q by setting, for z ∈ C,

∆(z) = τ n ∆ s 0 + z τ , (13)
P (z) = τ n P s 0 + z τ , (14)
Q(z) = e -s0τ τ n Q s 0 + z τ ( 15
∆(z) = n!z n+m+1 (n + m + 1)! Φ(m + 1, n + m + 2, -z). ( 16
) (d)
The coefficients a 0 , . . . , a n-1 , α 0 , . . . , α m of ∆ are given by

                   a k = (-1) n-k n! n j=k j k m+n-j m s j-k 0 j! τ n-j for k ∈ 0, n -1 , α k = (-1) n-1 e s0τ m j=k (-1) j-k (m + n -j)! s j-k 0 k! (j -k)! (m -j)! τ n-j for k ∈ 0, m . ( 17
)
In addition, if the above equivalent assertions are satisfied and m ≤ n, then s 0 is the rightmost root of ∆, i.e., we have ℜ(s) ≤ ℜ(s 0 ) for every s ∈ C such that ∆(s) = 0.

Proof. Clearly, thanks to the definition of ∆, assertion (a) is equivalent to requiring 0 to be a root of multiplicity D P S of ∆. Since ∆ is an analytic function and D P S is the maximal possible multiplicity of any of its roots, the latter fact is equivalent to requiring that ∆(z) = O(z n+m+1 ) as z → 0. Using that ∆(z) = P (z) + e -z Q(z), the latter is equivalent to

P (z)e z + Q(z) = O(z n+m+1
) as z → 0, which, thanks to the definition of Padé approximation, is equivalent to (b).

Using the uniqueness of the Padé approximation of order (m, n) with a monic denominator, we deduce that (b) is equivalent to requiring that -Q = Q n,m and P = P n,m , where P n,m and Q n,m are defined in (1), i.e., it is equivalent to requiring that, for every z ∈ C,

∆(z) = P n,m (z) -e -z Q n,m (z).
In particular, equivalence with (d) follows after straightforward but long computations by using (13) (see, e.g., (Mazanti et al., 2021a, Lemma 4.2) for detailed computations in the case m = n-1, which can be easily generalized to any m).

Finally, if (b) is satisfied, then, combining (2), ( 5), and (6), we get that

e z ∆(z) = P n,m (z)(e z -R n,m (z)) = n!z n+m+1 (n + m + 1)! Φ(n + 1, n + m + 2, z) = e z n!z n+m+1 (n + m + 1)! Φ(m + 1, n + m + 2, -z),
yielding (c). Conversely, if (c) holds, then, using ( 7) and ( 16), we deduce that

P (z)e z + Q(z) = e z ∆(z) = O(z n+m+1
), yielding (b) thanks to the definition of Padé approximation. Remark 3.2. The equivalence between items (a), (c), and (d) of Theorem 3.1 was already stated and proved in [START_REF] Boussaada | The generic multiplicity-induced-dominancy property from retarded to neutral delay-differential equations: when delay-systems characteristics meet the zeros of Kummer functions[END_REF] by exploiting the integral representation (5). By exploiting the Padé approximation of the exponential function and in particular (2), Theorem 3.1 provides the additional equivalence with (b) and allows for a much simpler proof of the equivalences between (a)-(d).

Note that, by considering the first equation in ( 17) with k = n -1, one obtains the simple and interesting relation between s 0 , τ , and a n-1 given by

s 0 = - a n-1 n - m + 1 τ . (18)
As remarked in [START_REF] Boussaada | The generic multiplicity-induced-dominancy property from retarded to neutral delay-differential equations: when delay-systems characteristics meet the zeros of Kummer functions[END_REF], if any of the equivalent assertions of Theorem 3.1 is satisfied, then s 0 is the unique real root of ∆, and, more precisely, it is the unique root of ∆ on the horizontal strip {s ∈ C | |ℑ(s)| < 2π τ } of the complex plane. In addition, the assumption in Theorem 3.1 of requiring s 0 to be real is justified by the fact that nonreal roots of ∆ cannot have a multiplicity equal to the Pólya-Szegő bound D P S , since any root s 0 of ∆ attaining the maximal multiplicity D P S necessarily satisfies ( 18), and thus it will be real since a n-1 is real.

As an immediate consequence of Theorem 3.1, if ( 17) is satisfied for some s 0 ∈ R, then the trivial solution of ( 11) is exponentially stable if and only if a n-1 > -n (m+1) τ .

Note that Theorem 3.1 provides necessary and sufficient conditions for a real number s 0 to be a root of maximal multiplicity of the quasipolynomial ∆ from ( 12). The main result of this section states that, under those conditions, s 0 is necessarily a dominant root of ∆.

ADDITIONAL COMMENTS ON THE LOCATION OF ROOTS OF WHITTAKER FUNCTIONS

Proposition 2.6 is actually a corrected version of a result by G. E. Tsvetkov, (Tsvetkov, 1941, Theorem 7). In [START_REF] Tsvetkov | On the roots of Whittaker functions[END_REF], the author provides statements of results concerning the location of zeros of Whittaker functions without the corresponding proofs, pointing only to the technique by E. Hille from [START_REF] Hille | Oscillation theorems in the complex domain[END_REF]. It turns out, as highlighted in (Boussaada et al., 2022a, Counterexample 3.1), that there are counterexamples to the statement of (Tsvetkov, 1941, Theorem 7). Proposition 2.6, proved in Boussaada et al. (2022a), uses the technique suggested by G. E. Tsvetkov in his paper, based on the Green-Hille transform, to provide a corrected version of (Tsvetkov, 1941, Theorem 7).

Some works in the literature have relied on (Tsvetkov, 1941, Theorem 7) to provide additional properties on the location of roots of Whittaker functions. This is the case, for instance, of (Saff and Varga, 1978, Proposition 3.2), which studies Whittaker functions motivated by the previously described links with Padé approximations of the exponential function. More precisely, the statement of [START_REF] Saff | On the zeros and poles of Padé approximants to e z[END_REF], Proposition 3.2) is the following. Proposition 4.1. [START_REF] Saff | On the zeros and poles of Padé approximants to e z[END_REF]). Let k ∈ R, l > 0, and z be a nontrivial zero of M k,l .

(a) If k > 0, then ℜ(z) > 2k and ℑ(z

) 2 > ℜ(z) 2 4l 2 -4k 2 -1 4k 2 . (b) If k < 0, then ℜ(z) < 2k and ℑ(z) 2 > ℜ(z) 2 4l 2 -4k 2 -1 4k 2 . (c) If k = 0, then ℜ(z) = 0 and ℑ(z) 2 > 4l 2 -1.
It turns out that the counterexample to (Tsvetkov, 1941, Theorem 7) from (Boussaada et al., 2022a, Counterexample 3.1) also provides a counterexample to items (a) and (b) of Proposition 4.1. For the sake of completeness, let us detail the counterexample.

Let l > 0 and take k = l + 3 2 . If z is a nontrivial root of M k,l , then, by (8), z is a nontrivial root of Φ(-1, 1 + 2l, •). From (3), we have Φ(-1, 1+2l, z) = 1-z 1+2l , and its unique root is z = 1+2l. In particular, ℜ(z) = 1+2l < 2k = 3+2l, and thus Proposition 4.1(a) is not verified. In addition, as a consequence of Proposition 2.5, z = -1 -2l is the unique nontrivial root of M -k,l and, since ℜ

(z) = -1 - 2l > -2k = -3 -2l, Proposition 4.1(b) is not verified.
In view of Proposition 4.1 and Proposition 2.6, a natural question is whether one may provide, in the case k > 0, a lower bound on the real part of nontrivial roots of M k,l as a function of k. A numerical exploration of such a question was provided in Boussaada et al. (2022a) by considering the root z = 1 + 2l of M k,l for k = l + 3 2 . Note that such a root exists whenever l > -1 2 and it is real and simple. Hence, for every l > -1 2 , there exists an interval I l containing 1 + 2l and a curve k → z l (k) ∈ R defined on I l such that, for every k ∈ I l , z l (k) is a real root of M k,l with z l (l + 3

2 ) = 1 + 2l. The article Boussaada et al. (2022a) provides the numerical computation of these curves for l ∈ -1 4 , 0, 1 4 , 1 2 , 3 4 , 1 , which we reproduce in Figure 1. The black dots (k, z) correspond to k = l + 3 2 and the root z = 1 + 2l. .

The particular case n = 1, Ξ 1 = ζ ∈ R tan ζ 2 = ζ 2 , had already been identified in [START_REF] Mazanti | Effects of roots of maximal multiplicity on the stability of some classes of delay differentialalgebraic systems: The lossless propagation case[END_REF].

CONCLUDING REMARKS

This paper discusses the existing links between roots of generic maximal multiplicity for a quasipolynomial, Kum-mer confluent geometric functions and Padé approximations of the exponential function. More precisely, we have shown that in the case of a real characteristic root having the maximum multiplicity of a simple quasipolynomial, the corresponding polynomials define a pair of Padé approximations. Moreover, the multiplicity of the root is equal to the Pólya-Szegő bound and it is dominant in the sense that it explicitly defines the spectral abscissa of the spectrum of the corresponding delay-differential equations. As a byproduct of the analysis, we have slightly corrected some of the results from the literature on the Padé approximation.

Notations:

  Throughout the paper, the following notations are used: N * denotes the set of positive integers and N = N * ∪ {0}. The set of all integers is denoted by Z and, for a, b ∈ R, we denote a, b = [a, b] ∩ Z, with the convention that [a, b] = ∅ if a > b. For a complex number s, ℜ(s) and ℑ(s) denote its real and imaginary parts, respectively. The open left and right complex halfplanes are the sets C -and C + , respectively, defined by C -= {s ∈ C | ℜ(s) < 0} and

  and we have in addition Φ(a, b, z) = 1 + O(z) as z → 0. (7) Kummer confluent hypergeometric functions have close links with Whittaker functions, defined as follows (see, e.g., Olver et al. (2010)). Definition 2.4. Let k, l ∈ C and assume that 2l is not a negative integer. The Whittaker function M k,l is the function defined for z ∈ C by

  ) Then the following assertions are equivalent. (a) The real number s 0 is a root of multiplicity D P S = n+ m + 1 of ∆. (b) The pair of polynomials (-Q, P ) forms a Padé approximation of order (m, n) of the exponential function z → e z . (c) The quasipolynomial ∆ is given by

Fig. 1 .

 1 Fig. 1. Real root z(k) of M k,l satisfying z(l + 3 2 ) = 1 + 2l for six different values of l.An inspection of Figure1leads to the conjecture that the maximal interval I l on which z l is defined isI l = (l + 1 2 , +∞) and that z l (k) → 0 as k → +∞ and z l (k) → +∞ as k → l + 1 2 .In particular, if this conjecture is true, then one cannot expect to correct Proposition 4.1(a) by replacing the term 2k by any function of k which remains lower bounded as k → +∞.Finally, we notice the remarkable case depicted in[START_REF] Boussaada | The generic multiplicity-induced-dominancy property from retarded to neutral delay-differential equations: when delay-systems characteristics meet the zeros of Kummer functions[END_REF], where one is able to characterize the zeros of the Whittaker function M 0,n+ 1 2 for n ∈ N as {iζ | ζ ∈ Ξ n }, where Ξ n is the set of ζ ∈ R satisfying


Keywords: Delay, stability, controller design, Python toolbox, GUI, online software

This paper presents the software entitled "Partial Pole Placement via Delay Action," or "P3δ" for short. P3δ is a Python software with a friendly user interface for the design of parametric stabilizing feedback laws with time-delays for dynamical systems. After recalling the theoretical foundation of the so-called "Partial Pole Placement" methodology we propose as well the main features of the current version of P3δ. We illustrate its use in feedback stabilization of several control systems operating under time delays.

New Features of P3δ Software.

Insights and Demos ⋆ Islam Boussaada * , * * Guilherme Mazanti * Silviu-Iulian Niculescu * Ayrton Hammoumou * * Titouan Millet * * Jayvir Raj * * Julien Huynh * * 1. INTRODUCTION Time delays often occur in controlling dynamical systems, mainly due to the time required for acquiring, propagating, or processing information. It is commonly accepted that a delay in a control loop induces instability, oscillations and bad performance of the overall scheme. For instance, since the 30s, [START_REF] Callender | Timelag in a control system[END_REF], [START_REF] Hartree | Time-lag in a control system[END_REF] showed the difficulty of handling delays in control loops for second-and third-order linear time-invariant (LTI) dynamical systems and one of the natural ideas was to compensate it by pre-correction (see, e.g., [START_REF] Porter | Introduction to servomechanisms[END_REF]).

However, as noted and briefly discussed in [START_REF] Sipahi | Stability and stabilization of systems with time delay: limitations and opportunities[END_REF], in some cases, the delay can have a stabilizing effect. Furthermore, it has been emphasized in [START_REF] Suh | Proportional minus delay controller[END_REF] and [START_REF] Atay | Balancing the inverted pendulum using position feedback[END_REF] that one may replace the classical proportional-derivative controller by a proportionaldelayed controller, using the so-called "average derivative action" due to the delay. Finally, regarding the beneficial effect of the delay, closed-loop stability may be guaranteed for some dynamical systems subject to input delays precisely by the existence of such delays, as pointed out in [START_REF] Niculescu | Delay effects on output feedback control of dynamical systems[END_REF] and the references therein in controlling oscillators by delayed output feedback.

The above reasons explain, in part, an abundant literature on such topics, such as [START_REF] Gu | Stability of time-delay systems[END_REF]; [START_REF] Michiels | Stability, control, and computation for time-delay systems: An eigenvaluebased approach[END_REF]; [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF] and the references therein. It should be mentioned that time-delay dynamical systems are infinite-dimensional systems and there exist several ways to represent their dynamics. In the sequel, their dynamics are represented by delay-differential equations ⋆ The authors wish to acknowledge the work of the full P3δ development team, which, in addition to the authors, also include Yoann Audet, Thomas Charbonnet, Honoré Curlier, Pierre-Henry Poret, Achrafy Said Mohamed and Franck Sim. The development of P3δ was also made possible thanks to the work of the Cyb'Air Association.

(DDEs). For an introduction on DDEs, we refer to [START_REF] Hale | Introduction to functional differential equations[END_REF].

The software1 P3δ was introduced in Boussaada et al. (2020a[START_REF] Boussaada | New features of P3delta software: Partial pole placement via delay action[END_REF]. The main intention of the authors is to help users interested in stability analysis and stabilization of dynamical LTI systems in the presence of a single delay in closed loop. P3δ makes use of the so-called multiplicity-induced-dominancy (MID) and partial pole placement methods. More precisely, MID is an unexpected spectral property stating that, in some cases, the characteristic root with maximal multiplicity defines the spectral abscissa of the corresponding characteristic function, i.e., the rightmost root of the spectrum. In control, MID opened a new and interesting perspective relying on the idea of a partial pole placement by reinforcing the use of the delay as a control parameter. For a deeper discussion on these methods, we refer to [START_REF] Amrane | On qualitative properties of low-degree quasipolynomials: further remarks on the spectral abscissa and rightmost-roots assignment[END_REF]; [START_REF] Bedouhene | Real spectral values coexistence and their effect on the stability of time-delay systems: Vandermonde matrices and exponential decay[END_REF]; Boussaada and Niculescu (2018); [START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF][START_REF] Boussaada | Further remarks on the effect of multiple spectral values on the dynamics of time-delay systems. Application to the control of a mechanical system[END_REF]; [START_REF] Mazanti | Multiplicity-induced-dominancy for delay-differential equations of retarded type[END_REF]Mazanti et al. ( , 2020a,b),b); Boussaada and Niculescu (2016a,b); Boussaada et al. (2016Boussaada et al. ( , 2018a)); [START_REF] Ma | PID control design for firstorder delay systems via MID pole placement: Performance vs. robustness[END_REF].

The software P3δ covers DDEs of retarded or neutral type2 with a single time delay, under the form

y (n) (t) + a n-1 y (n-1) (t) + • • • + a 0 y(t) + b m y (m) (t -τ ) + • • • + b 0 y(t -τ ) = 0, (1)
under appropriate initial conditions, where τ > 0 is the positive delay, y is the real-valued unknown function, n and m are nonnegative integers with n ≥ m, and a 0 , . . . , a n-1 , b 0 , . . . , b m are real coefficients.

To address the stability analysis of LTI DDEs, the software P3δ relies on spectral methods (see, e.g., [START_REF] Hale | Introduction to functional differential equations[END_REF]; [START_REF] Michiels | Stability, control, and computation for time-delay systems: An eigenvaluebased approach[END_REF]), which consist on the study of the complex roots of a characteristic function of the system. The characteristic function ∆ : C → C of (1) is given by

∆(s) = s n + n-1 k=0 a k s k + e -sτ m k=0 b k s k , (2)
and ( 1) is exponentially stable if and only if the spectral abscissa γ = sup{Re s | ∆(s) = 0} satisfies γ < 0.

P3δ considers that b 0 , . . . , b m are free parameters. In its "Generic" mode, a 0 , . . . , a n-1 are also assumed to be free and τ is fixed while, in its "Control-oriented" mode, a 0 , . . . , a n-1 are assumed fixed and τ can be assumed either free or fixed. The "Control-oriented" mode is usually suitable for control applications, as illustrated in Section 3.

The strategy used by P3δ to stabilize a time-delay system is to tune its free parameters to assign finitely many roots while also guaranteeing that the rightmost root on the complex plane is among the chosen ones. Two main properties define such a strategy: (i) first, assigning a real root of maximal multiplicity and proving that this root is necessarily the rightmost root of the characteristic quasipolynomial, a property which has been named multiplicityinduced-dominancy, or MID for short, and (ii) second, assigning a certain amount of real roots, typically equally spaced for simplicity, and proving that the rightmost root among them is also the rightmost root of the characteristic quasipolynomial, a property which has been named coexisting real roots-induced-dominancy, or CRRID for short, in [START_REF] Amrane | On qualitative properties of low-degree quasipolynomials: further remarks on the spectral abscissa and rightmost-roots assignment[END_REF]; [START_REF] Bedouhene | Real spectral values coexistence and their effect on the stability of time-delay systems: Vandermonde matrices and exponential decay[END_REF].

The MID property for (1) was shown, for instance, in [START_REF] Boussaada | Further remarks on the effect of multiple spectral values on the dynamics of time-delay systems. Application to the control of a mechanical system[END_REF] in the case (n, m) = (2, 0), in [START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF] in the case (n, m) = (2, 1) (see also Boussaada and Niculescu (2018)), in [START_REF] Mazanti | Multiplicity-induced-dominancy for delay-differential equations of retarded type[END_REF] In all the above cases, the maximal multiplicity of a real root or, equivalently, the maximal number of coexisting simple real roots is n+m+1. The idea to exploit the nature of (real or complex) open-loop roots in control design was proposed for second-order systems in [START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF] and extended for arbitrary order systems with real-rooted plants in [START_REF] Balogh | Towards an MID-based delayed design for arbitrary-order dynamical systems with a mechanical application[END_REF][START_REF] Balogh | Conditions for stabilizability of time-delay systems with real-rooted plant[END_REF].

As mentioned earlier, P3δ allows for the parametric design of stabilizing feedback laws with time-delays by exploiting the MID and CRRID properties briefly presented above. The present paper describes the new functionalities of P3δ and provides some illustrative examples for its use.

NEW FEATURES OF P3δ

New features of the online version

After a first iteration of the software, P3δ Online was completed with new features (Figure 1) that enrich the software and improve user experience. The one-click online version of P3δ is still hosted on Binder, which provides a personalized computing environment directly from a GitHub repository. The P3δ team continued to develop the online software in Python, using the dynamic Jupyter Notebook format and with an user interface built using interactive widgets from Python's ipywidgets module.

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

3
P3δ Online is based on the program of the executable version and enriched with exclusive features to this version.

The online software includes features from the "Generic MID", "Control-oriented MID", and "Generic CRRID" modes of P3δ described in [START_REF] Boussaada | New features of P3delta software: Partial pole placement via delay action[END_REF]Boussaada et al. ( , 2020a).

In the "Generic MID" mode, the online version of the software returns the spectrum distribution as well as a normalized quasipolynomial which admits a root of multiplicity n + m + 1 at the origin4 . In the "Controloriented MID" mode, the online version of software returns the admissibility region, a normalized quasipolynomial which admits a root of multiplicity m+2 at the origin, and an illustration of the bifurcation of the root of multiplicity m + 2 with respect to variations of the delay τ .

The first new feature of P3δ Online is the addition of a "Home" tab that contains information about the project and the software settings. Two settings are currently available: the appearance of notifications and the activation of the software limits.

In the "Generic MID" mode, the spectral distribution analysis now provides two new forms of the output equation, Factorized integral equation and Hypergeometric factorization. In the "Control-oriented MID" mode, it is now possible to choose the value of the discretization step as well as the number of iterations for the τ sensitivity plot.

Finally, the P3δ team has developed a new feature to export the results obtained in the software in the format of a report automatically written in a PDF file. The user has the choice of which calculation mode results they want to export and can save the report directly on their computer.

Assignment admissibility region

Given a system, the definition of valid hyper-parameters (s 0 and τ ) can sometimes be difficult. In order to help users getting a better idea of the admissibility ranges of s 0 and τ , P3δ has an admissibility region plotting feature which will be described in this section. Note that this feature is only available when using "Control-oriented MID" as it currently is the only mode with constraints on (s 0 , τ ). More precisely, given a 0 , . . . , a n-1 , the admissibility region is defined as the set of pairs (s 0 , τ ) ∈ R × (0, +∞) for which there exist real coefficients b 0 , . . . , b m such that s 0 is a root of ∆ of multiplicity at least m + 2 when the delay is τ . To compute such a region, we first express the coefficients (b i ) 0≤i≤m in terms of (a i ) 0≤i≤n-1 , s 0 , and τ by using the m + 1 equations ∆ (k) (s 0 ) = 0, k ∈ {0, . . . , m}.

As those equations are linear in the m + 1 variables b 0 , . . . , b m , this linear system admits a unique solution. We then replace these expressions of (b i ) 0≤i≤m in the equation ∆ (m+1) (s 0 ) = 0, obtaining an algebraic relation between s 0 , τ , and the coefficients (a i ) 0≤i≤n-1 . By construction, this is a necessary and sufficient condition for s 0 to be a root of multiplicity at least m + 2 of ∆ and, since the coefficients (a i ) 0≤i≤n-1 are known in this mode, this algebraic equation characterizes the admissibility region.

In P3δ, the computations leading to this admissibility region is done by following the previously mentioned steps in a symbolic way using the sympy package. Only the part of the admissibility region in the rectangle [s 0,min , 0] × [0, τ max ] is displayed, where s 0,min < 0 and τ max > 0 are values selected by the user.

APPLICATIONS OF P3δ

To illustrate P3δ, we present some applications of the software.

Harmonic oscillator

Consider a controlled harmonic oscillator described by y ′′ (t) + y(t) = u(t),

(3) where y(t) ∈ R is the instantaneous state of the oscillator available to measurement, the control input u(t) corresponds to the applied force, and the coefficients of the equation were normalized. We assume that the control input is given by a delayed proportional-derivative controller u(t) = -βy(t -τ ) -αy ′ (t -τ ), (4) where α and β are the coefficients of the controller and τ > 0 is the delay. The characteristic equation of the closed-loop system is thus ∆(s) = s 2 + 1 + (β + αs)e -sτ .

(5) The polynomial corresponding to the non-delayed term in ( 5) is of degree n = 2, while that corresponding to the delayed term is of degree m = 1. The degree of the quasipolynomial ∆ is thus n + m + 1 = 4.

Let us use P3δ in order to place a root of multiplicity m + 2 = 3 at some s 0 ∈ R. One should first input into P3δ the values of the known coefficients of the quasipolynomial ∆ from (5), i.e., the coefficients of the polynomial corresponding to the non-delayed term. After introducing the data corresponding to (5) into P3δ, one obtains the plot of the admissibility region, as in Figure 2. Based on Figure 2, one observes that the method used by P3δ allows for the stabilization of the system if the delay is at most approximately 1.4. After inputting the value of the delay, P3δ will compute the root s 0 and the coefficients α and β, and trace the spectrum distribution of (5) with the derived values. In the present example, the corresponding screen of P3δ is shown in Figure 3, where, for τ = 1, one obtains s 0 = -1, α = 0, and β ≈ -0.7358.

Fig. 3. Computed values and spectrum distribution for the harmonic oscillator.

Inverted pendulum

As a second example of application of P3δ, we consider the stabilization of an inverted pendulum from Figure 4, in which a stick of mass m and length ℓ is placed over a cart and can rotate freely around the attachment point A.

The cart can move along a rail and an external force F is applied on the cart, and we assume that the mass of the cart is negligible compared to the mass of the stick. We denote by φ the angle between the stick and the upward vertical direction. In this case, following [START_REF] Molnar | Calculation of the critical delay for the double inverted pendulum[END_REF], the linearization of the dynamics of φ around the unstable equilibrium φ = 0 is

φ ′′ (t) -mgℓ 2I φ(t) = ℓ 2I F (t), (6)
where I = 1 12 mℓ 2 is the moment of inertia of the stick and g is the gravitational acceleration. For the numerical application, we will consider m = 10 kg and ℓ = 10 m, in which case I ≈ 83.33 kg m 2 . We wish to stabilize the origin of ( 6) by a delayed PD feedback of the angular position, i.e., we wish to apply yielding a closed-loop system with characteristic equation

F (t) = -K p φ(t -τ ) -K d φ ′ (t -τ ),
∆(s) = s 2 -mgℓ 2I + Kpℓ 2I + K d ℓ 2I s e -sτ . (7)
Note that, here, the delay τ can be seen as a design parameter, together with K p and K d .

By introducing the data corresponding to (7) into P3δ "Control-oriented MID" mode, we obtain the admissibility plot from Figure 5. Based on this plot, we decide to assign a root of multiplicity 3 of ( 7) at s 0 = -5. P3δ then gives the result shown in Figure 6. From the computed values of the coefficients shown in Figure 6, one deduces that the delay corresponding to s 0 = -5 is τ ≈ 0.112 and the coefficients of ( 7) should be

Transonic flow in a wind tunnel

Consider now a wind tunnel in which a cold fluid is set to motion at high speed. The control of the velocity of such a fluid around an equilibrium state can be described by the system κm

′ (t) + m(t) = kθ(t -τ 0 ), θ ′′ (t) + 2ζωθ ′ (t) + ω 2 θ(t) = u(t), (8)
where m denotes the deviation of the Mach number of the fluid with respect to its equilibrium state, κ and k are constants depending on the characteristics of the fluid and the desired equilibrium state, θ is the angle of a guide vane driving the velocity of the fluid, τ 0 > 0 is a delay depending only on the temperature of the fluid, ζ and ω are parameters of the dynamics of the guide vane angle, and u is a control input. The above model comes from [START_REF] Armstrong | An application of multivariable design techniques to the control of the National Transonic Facility[END_REF] and its stabilization was previously discussed in [START_REF] Mazanti | Multiplicity-induced-dominancy for delay-differential equations of retarded type[END_REF] under the assumption that one may choose ζ and ω. We consider here a more realistic situation in which ζ and ω are fixed, and a feedback law of the form u(t) = -βm(t -τ 1 ) -α 0 θ(t -τ 0 -τ 1 ) -α 1 θ ′ (t -τ 0 -τ 1 ), where τ 1 ≥ 0 is a new delay, which can be seen as a design parameter and should be at least equal to the delay for measuring m. For simplicity, we denote τ = τ 0 + τ 1 . Inserting this control law into (8) and performing straightforward algebraic manipulations, one deduces that θ verifies the third-order differential equation

θ ′′′ (t) + 2ζω + 1 κ θ ′′ (t) + ω 2 + 2ζω κ θ ′ (t) + ω 2 κ θ(t) + α 1 θ ′′ (t -τ ) + α 0 + α1 κ θ ′ (t -τ ) + α0+βk κ θ(t -τ ) = 0, and hence the closed-loop characteristic quasipolynomial ∆ of (8) is ∆(s) = s 3 + 2ζω + 1 κ s 2 + ω 2 + 2ζω κ s + ω 2 κ + (γ 2 s 2 + γ 1 s + γ 0 )e -sσ , (9)
where γ 2 = α 1 , γ 1 = α 0 + α1 κ , and γ 0 = α0+βk κ .

As a numerical application, we consider the linearization around the steady state with Mach number 0.84 and air temperature 166 K, in which case, as reported in [START_REF] Armstrong | An application of multivariable design techniques to the control of the National Transonic Facility[END_REF], the system parameters are κ = 1.964 s, k = -0.670 36 rad -1 , τ 0 = 0.33 s, ζ = 0.4368, and ω = 3.292 rad s -1 . We can hence insert the known coefficients of (9) into P3δ, as shown in Figure 7, and obtain the admissibility plot from Figure 8. Based on Figure 8, we choose in this example s 0 = -2.94675 as a root of multiplicity 5 of the quasipolynomial ∆, and P3δ provides the output shown in Figure 9. In particular, we deduce the numerical value τ ≈ 0.4140 s, yielding τ 1 ≈ 0.0840 s, as well as γ 2 ≈ 0.04167, γ 1 ≈ -1.9005, and γ 0 ≈ 1.9993, yielding β ≈ -2.9908, α 0 ≈ 1.9217, and α 1 ≈ 0.04167.

CONCLUDING REMARKS AND FUTURE WORK

By exploiting the properties of the MID and the CRRID, P3δ allows for the design of feedback control laws for real

  in the case of any positive integer n and m = n -1 (see alsoMazanti et al. (2020a)), and more recently in[START_REF] Boussaada | The generic multiplicity-induced-dominancy property from retarded to neutral delay-differential equations: When delay-systems characteristics meet the zeros of Kummer functions[END_REF] for arbitrary n ≥ m. It was also studied for neutral systems of orders 1 and 2 in[START_REF] Ma | PID control design for firstorder delay systems via MID pole placement: Performance vs. robustness[END_REF];[START_REF] Benarab | The MID property for a second-order neutral time-delay differential equation[END_REF] and extended to complex conjugate roots of maximal multiplicity in[START_REF] Mazanti | Spectral dominance of complex roots for single-delay linear equations[END_REF]. The CRRID property was shown, for instance, in[START_REF] Amrane | On qualitative properties of low-degree quasipolynomials: further remarks on the spectral abscissa and rightmost-roots assignment[END_REF] in the cases (n, m) = (2, 0) and (n, m) = (1, 0), and in[START_REF] Bedouhene | Real spectral values coexistence and their effect on the stability of time-delay systems: Vandermonde matrices and exponential decay[END_REF] in the case of any positive integer n and m = 0.

Fig. 1 .

 1 Fig. 1. Illustration of the different new features of P3δ Online.

Fig. 2 .

 2 Fig. 2. Data input and admissibility plot for the harmonic oscillator.

Fig. 4 .

 4 Fig. 4. Inverted pendulum on a cart.

Fig. 5 .

 5 Fig. 5. Admissibility plot for the inverted pendulum.

Fig. 6 .

 6 Fig. 6. Computed values and spectrum distribution for the inverted pendulum.

Kpℓ2I≈

  11.53 and K d ℓ 2I ≈ 4.4898, yielding the parameters K p ≈ 192.16 and K d ≈ 74.83.

Fig. 7 .

 7 Fig. 7. Data input for the transonic flow model.

Fig. 8 .

 8 Fig. 8. Admissibility plot for the transonic flow model.

Fig. 9 .

 9 Fig. 9. Computed values and spectrum distribution for the transonic flow model. applications. The main novelty of the version discussed in the paper is the improvement of the graphic interface for the online version of the software. Inspired by the "Control-oriented MID" mode, the team is working on a "Control-oriented CRRID" mode. In addition, an executable version for macOS is currently under development.

An acronym for Partial Pole Placement via Delay Action

For the classification of DDEs, we refer to[START_REF] Hale | Introduction to functional differential equations[END_REF].

Interested readers may also contact directly any of the authors of the paper.

In this case, the multiplicity coincides with the degree of the quasipolynomial ∆, and represents the maximal allowable multiplicity.

ACKNOWLEDGMENTS

This work is partially supported by a grant 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-IDEX0003-02. The authors also acknowledge the support of Institut Polytechnique des Sciences Avancées (IPSA).