Yacine Bouzidi
email: yacine.bouzidi@inria.fr

Adrien Poteaux
email: adrien.poteaux@univ-lille.fr

Real curve analysis and stability of time-delay systems

Keywords: Time delay systems, curve analysis, stability, certified tests, symbolic-numeric methods

In this work, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose dynamics are defined in the frequency domain by quasi-polynomials of the form f (s, τ ) = m j=0 p j (s) e -j τ s . We propose a new approach for studying the stability with respect to the delay: we determine the asymptotic behavior of the roots of the quasi-polynomial near to the critical pairs of f (s, τ ) by analysing the intersection of an associated algebraic space curve with sufficiently refined 3D real boxes. Compared to the existing methods, our method does not require any Puiseux series computation.

INTRODUCTION

In this paper, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose dynamics are defined by the following state-space representation

ẋ(t) = m k=0 A k x(t -k τ ) (1)
where τ ∈ R + := {τ ∈ R | τ ≥ 0}, m ∈ N and A 0 , . . . , A m are constant matrices with entries in a field K (e.g. Q or R). The characteristic function of (1) is defined by a quasipolynomial of the form

f (s, τ ) = m j=0 p j (s) e -j τ s , (2)
where the p j 's are polynomials in the complex variable s with coefficients in K.

Such a system is said to be asymptotically stable if all the zeros of the quasi-polynomial f (s, τ ) have negative real parts, i.e., f (s, τ ) = 0 for all s in the closed right half-plane C + := {s ∈ C | (s) ≥ 0}. Considering the delay τ as a parameter, the stability property will obviously depend on the value of τ and a classical question is then to study the asymptotic stability with respect to the delay value τ .

A classical approach for studying the asymptotic stability depending on the delay is based on the analysis of the socalled critical pairs of f (s, τ ), that is the pairs (ω ∈ R, τ ∈ R + ) such that f (i ω, τ ) = 0 ; see for instance [START_REF] Gu | Stability of Time-Delay Systems[END_REF], [START_REF] Li | Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays[END_REF], [START_REF] Marshall | Time-Delay. Systems: Stability and Performance Criteria with Applications[END_REF], [START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF], [START_REF] Olgac | An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems[END_REF] and the references therein.

If such critical pairs exist, the stability is derived from the way these critical imaginary roots behave under small variations of the time-delay τ . As a result, this approach yields a decomposition of the τ -axis into a set of disjoint intervals in which the considered system is either stable or unstable.

Algorithmically, the above approach can be decomposed in two independent problems. First, the computation of the critical pairs of f (s, τ ), and second, for each critical pair (s 0 , τ 0 ), the study of the asymptotic behavior of the root s 0 with respect to a small variation of the delay τ 0 .

There exist numerous methods for computing the critical pairs of f (s, τ ) (see for instance [START_REF] Bouzidi | Computer algebra methods for the stability analysis of differential systems with commensurate time-delays[END_REF], [START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF], [START_REF] Olgac | An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems[END_REF]). In [START_REF] Bouzidi | Computer algebra methods for the stability analysis of differential systems with commensurate time-delays[END_REF] for example, based on an algebraic formulation for the critical pairs of f (s, τ ), the problem is reduced to solving a system of algebraic equations in two variables and classical symbolic methods are used to obtain certified numerical approximations for the critical pairs. Such approximations consist in a pair of intervals isolating the two coordinates s and τ of the critical pairs.

For the study of the asymptotic behavior of critical imaginary roots, several methods have been developed for the case of simple imaginary roots (see, e.g., [START_REF] Gu | Stability of Time-Delay Systems[END_REF], [START_REF] Li | Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays[END_REF], [START_REF] Marshall | Time-Delay. Systems: Stability and Performance Criteria with Applications[END_REF], [START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF], [START_REF] Olgac | An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems[END_REF] and the references therein). For multiple imaginary roots (MIR), recent methods based on Puiseux series was developed in [START_REF] Bouzidi | Computer algebra methods for the stability analysis of differential systems with commensurate time-delays[END_REF], [START_REF] Li | Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays[END_REF][START_REF] Li | On computing puiseux series for multiple imaginary characteristic roots of lti systems with commensurate delays[END_REF]. These methods require in particular the symbolic computation of the Puiseux series of f (s, τ ) at all the critical pairs.

Contributions and outline of the paper. In this paper, we assume that the critical pairs of f (s, τ ) are given as certified numerical approximations (using for instance the method in [START_REF] Bouzidi | Computer algebra methods for the stability analysis of differential systems with commensurate time-delays[END_REF]) and we focus on the analysis of the asymptotic behavior of s at these critical pairs. We propose a new approach that avoids any Puiseux series computation at the critical pairs. Instead, the asymptotic behavior is derived through the analysis around the critical pairs of the real branches of a space algebraic curve obtaind by local approximation of the quasi-polynomial f (s, t). This analysis is done while avoiding as much as possible symbolic computations, using interval arithmetic instead. More precisely, given a critical pair (ω 0 , τ 0 ), our approach proceeds in the following way:

(1) We compute an upper bound m for the multiplicity of (ω 0 , τ 0 ) by evaluating the successive derivatives of f up to the point that we can certify a non-zero evaluation. See Section 3.1. (2) We consider f (s + ω 0 , τ + τ 0 ) as en element of C [[s, τ ]] by replacing the exponential by their local Taylor series. Note that its coefficients are given with error bounds. To make computations possible, we truncate these series modulo some powers of s and τ , taking at least m + 1 as a truncation bound. This provides a polynomial p(s, τ ) that is an approximation of the studied quasi-polynomial (note that approximations are coming from both numerical errors and the truncation process). See Section 3.2. (3) Finally, we study the local behaviour of the real curve defined by {p r (s r , s i , τ ) = p i (s r , s i , τ ) = 0} where p r and p i are respectively the real and the imaginary part of p(s r + i s i , τ ). More precisely, we isolate the curve in a 3D box that contains only the germ at (0, 0, 0) of the real curve; then, one can conclude about the stability of the system by looking at the intersections of real branches with this box.

Organisation of the paper. For a sake of clarity, we briefly recall in Section 2 the method presented in [START_REF] Bouzidi | Computer algebra methods for the stability analysis of differential systems with commensurate time-delays[END_REF] for the computation of the critical pairs of f (s, t). Then, in Section 3, we show how a quasipolynomial can be locally approximated at a given point by a polynomial while preserving the local topology of the latter. Using this approximated curve, we present in Section 4 a strategy to compute the asymptotic behavior of this curve at a given point by looking at the intersection of the latter with some specific 3D boxes. This enables us to study the stability of the system.

CRITICAL PAIRS OF A QUASI-POLYNOMIAL

As mentioned in the introduction, to study the asymptotic stability of a system of the form (1), we first need to compute the critical pairs of the corresponding quasipolynomial f (s, τ ). These critical pairs are the intersection points of the curve defined by the quasi-polynomial, f (s, τ ) = 0, with the imaginary axis in s, that is

{(ω ∈ R, τ ∈ R + ) | f (i ω, τ ) = 0}
We now briefly recall the method given in [START_REF] Bouzidi | Computer algebra methods for the stability analysis of differential systems with commensurate time-delays[END_REF] for the certified computation of these critical pairs. The main idea is to use the so-called Rekasius transformation, that allows to recover the infinite number of zeros of f (i ω, τ ) from the real zeros of a zero-dimensional algebraic bivariate system. More precisely, we replace the terms e -τ i ω in the quasipolynomial f (i ω, τ ) by the rational fraction 1-γ i ω 1+γ i ω with γ ∈ R. Cleaning the denominators, one obtain a complex polynomial of the form R(ω, γ) + i I(ω, γ). Computing the critical pairs then amounts to identifying the real zeros of the following bivariate polynomial system: R(ω, γ) = 0,

I(ω, γ) = 0, (3)
and from each solution (ω, γ) of (3) the critical delays can then be obtained by:

τ k = 2 ω (arctan(ω γ) + k π), k ∈ Z. (4)
The method proposed in [START_REF] Bouzidi | Computer algebra methods for the stability analysis of differential systems with commensurate time-delays[END_REF] for computing the real solutions of the above algebraic system consists first in computing a univariate representation of the solutions (5), which is roughly speaking, a one-to-one mapping between the solutions of the system (3) and the root of a univariate polynomial f (t) (see [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF] for details), then, numerical approximations for the real solutions are obtained from this representation.

f (t) = 0 ω = g ω /g 0 γ = g γ /g 0 (5)
From the above representation, one can compute intervals with rational endpoints that isolate the coordinates ω and γ of the solutions. Moreover, by substituting the obtained intervals in (4), one also obtain isolating intervals for the corresponding delays τ . Another advantage of the previous representation is that it allows to control the precision of the numerical approximation of the solutions. Indeed, given an isolating interval for a root of the polynomial f (t), one can refine this interval up to a desired precision while it still contains the root which results to refined intervals for ω, γ and τ . As we will see in the sequel, this feature is very useful in the description and the design of our algorithm.

LOCAL APPROXIMATION

An upper bound on the multiplicities of critical pairs

An important information that we need for the description of our algorithm is the multiplicities of the critical pairs as roots of the quasi-polynomial f (s, τ ). Computing explicitly such multiplicities is however a challenging problem which is yet not solved (for detail, see [START_REF] Boussaada | Tracking the algebraic multiplicity of crossing imaginary roots for generic quasipolynomials: A vandermonde-based approach[END_REF] and references therein). Therefore, instead of exact multiplicities, we only compute upper bounds on the latter, which will be sufficient for our purpose. Let first recall the definition of the multiplicity of a critical pair of a quasi-polynomial f (s, τ ). Definition 1. A root (ω 0 , τ 0 ) of the quasi-polynomial f (i ω, τ ) is said to be of multiplicity k if and only if

f (ω 0 , τ 0 ) = ∂f ∂s (ω 0 , τ 0 ) = • • • = ∂ k-1 f ∂s k-1 (ω 0 , τ 0 ) = 0 and ∂ k f ∂s k (ω 0 , τ 0 ) = 0
According to the previous definition, if the evaluation at the point (ω 0 , τ 0 ) can be achieved in an exact way, a straightforward approach for computing the multiplicity of a critical pair (ω 0 , τ 0 ) is to evaluate it on the successive derivatives of f (s, τ ) until a non-zero value is obtained.

The key point of this paper is to avoid as much as possible to run costly symbolic computations. In particular, we are not going to compute the multiplicities of the critical pairs, but only a bound for each of them. To do so, we will just evaluate the successive derivatives ∂f k ∂ k s on the given 2D box (using interval arithmetic -see [START_REF] Moore | Methods and applications of interval analysis[END_REF] for details), until we can conclude that the m-th derivative is non zero in this box. In practice, this strategy will often provide the correct multiplicity, but over-bound might happen, so that we cannot assume to know the correct multiplicity.

In particular, if we get a bound m = 1 for a critical pair, then this multiplicity is correct and we can conclude that the system is not stable.

Local approximation

Given a critical pair (i ω 0 , τ 0 ), the quasi-polynomial f (s, τ ) can be considered in K[[s -i ω 0 , τ -τ 0 ]]. Thus, it can be approximated by a polynomial via a truncation of the powers of (s -i ω 0 ) and (τ -τ 0 ). In practice, we only have to use the Taylor series of the exponential to convert our quasi-polynomial in series. More precisely, it is well known that e x can be approximated by

n-1 k=0 x k k!
with an error bounded by x n n! e x . We thus can approximate f (s + i ω 0 , τ + τ 0 ) as m j=0 p j (s + i ω 0 )e -j i ω0τ0 e -j i ω0τ e -j τ0 s e -j s τ by replacing each last three exponentials by their truncated Taylor series (with x being respectively -j i ω 0 τ , -j τ 0 s and -j s τ ). This provides us a polynomial p(s, τ ) for which there are two different approximations:

• the truncation we just made (for which we have a good majoration of the error), that can be improved by increasing the truncation order, • the approximations we have for τ 0 and ω 0 (p having coefficients defined by interval), that we can improve on demand from the previous computations.

Note that we will always take a truncation bound n greater than the bound m we have for the multiplicity of the critical pair.

ASYMPTOTIC BEHAVIOR

As shown in Section 2, one can compute a set of disjoint 2D boxes in R 2 that isolate all the critical pairs of f (s, τ ). Moreover, thanks to the symbolic univariate representation of these pairs (5), the corresponding boxes can be refined up to an arbitrary precision. Also, for a given critical pair (ω 0 , τ 0 ) of f (s, τ ), Section 3 provides a local approximation of the quasi-polynomial f (s, τ ) by a two variables polynomial p(s, τ ) that preserves its local topology.

Given such a polynomial p(s, τ ), as well as a critical pair of p (a zero of the form (i ω 0 , τ 0 )), the goal of this section is to show how to analyze the asymptotic behavior of this critical pair. As a consequence, provided that this polynomial results from an appximation of a quasi-polynomial, this will naturally provide the asymptotic behavior of the critical pair as a zero of the quasi-polynomial.

Asymptotic behavior through curve intersection in R 3

Let us consider a polynomial p(s, τ ) where s ∈ C and τ ∈ R. In the following, we set s = x + i y, where x, y ∈ R.

For simplicity, we also denote to the variable τ by z.

Substituting s = x+i y as well as τ = z in the polynomial p yields p(x + i y, z) = g(x, y, z) + i h(x, y, z) where g, h are polynomials in three variables with real coefficients. With this notation, we define a critical pair of p as a point (0, y 0 , z 0 ) ∈ R 3 such that p(i y 0 , z 0 ) = 0, or equivalently, as the intersection point of the axis x = 0 with curve C defined by C := g(x, y, z) = 0 h(x, y, z) = 0 and whose coordinate in y (resp. z) is y 0 (resp. z 0 ). In the latter case, and to avoid confusion with critical points of the curve, we refer to it as a zero-crossing point. Moreover, for simplicity, we assume in the sequel that the curve C does not admit a whole component lying in a plan orthogonal to the x-axis. This translates into the fact that for any fixed value x 0 of x, the algebraic system {g(x 0 , y, z) = h(x 0 , y, z) = 0} admits a finite number of complex zeros.

Given a zero-crossing point (0, y 0 , z 0 ), the objective in the following is to analyse the asymptotic behavior of the coordinate x = 0 under a small variation of the coordinate z 0 . We suppose that this zero-crossing point (0, y 0 , z 0 ) is given by a 3D box in R 3 that isolates it from other zero-crossing points of the curve. Such a box can be obtained for example by extending the 2D box obtained from the representation (3) on the x-axis.

First, we need to introduce the notion of x-critical point of C. A point (x, y, z) of C is called x-critical if it is a critical point of the following projection map

Π x : C → R (x, y, z) → x
Algebraically, a point of the curve C is an x-critical point if and only if it is a zero of the polynomial ∂g ∂y ∂h ∂z -∂g ∂z ∂h ∂y . In the following, for a sake of briefness, we denote to this polynomial simply by d(x, y, z).

As a consequence, we can identify two different kinds of zero-crossing points of C. Definition 2. Let (0, y 0 , z 0 ) be a zero-crossing point of C. This point is called x-critical if and only if d(0, y 0 , z 0 ) = 0. Otherwise, it is called regular.

Given a zero-crossing point (0, y 0 , z 0 ) of C, if the latter is a regular zero-crossing points, then the analysis of the asymptotic behavior near to it is straightforward. Indeed, at the neighborhood of these kind of points, the curve C consists in one regular branch that passes through the point (0, y 0 , z 0 ). As a consequence, the variation of x can be deduced from the sign of d(0, y 0 , z 0 ). As (0, y 0 , z 0 ) is not an x-critical point, one can compute a box that contains the latter (in fact, it is simply a 2D box in y and z), and refine it until the interval resulting from the evaluation d(B) does not contain zero.

For x-critical zero-crossing points, analyzing the asymptotic behavior of C is more involved. In broad words, our approach consists in computing a 3-D box B that isolate this point from other x-critical points of the curve C. Then, this box B is refined until all the intersections of C with the boundary of B are located on the sides that are orthogonal to the x-axis. This guarantees that all the branches that intersect the boundary of the refined box pass through the considered zero-crossing point. The behavior is then easily deduced from these intersections. Definition 3. A box B containing a point p := (0, y 0 , z 0 ) of C is called isolating for p if and only if all the real branches of C in B pass through the point p. Theorem 1. Let B be a box that contains one x-critical zero-crossing point (0, y 0 , z 0 ) of C and such that the two following conditions hold • B does not contain any other x-critical point of C • The only intersection points of C with B are located on the sides of B that are orthogonal to the x-axis.

Then, B is an isolating box for (0, y 0 , z 0 ).

Proof. This is a direct consequence of the fact that B contains a single zero-crossing point of C and no other xcritical points. Indeed, since all the branches of C intersect B on the sides that are orthogonal to the x-axis, if one of these branches does not pass through (0, y 0 , z 0 ), this implies that either it intersects the axis x = 0 at y 1 = y 0 , or it passes through an x-critical point, which contradicts our assumption. Remark 2. Figures 1 and 2 below show the projection of two non-isolating 3D boxes on the (x -z)-plan. The first situation corresponds to a box B that contains an x-critical zero-crossing point (lying on the branch depicted in blue) as well as another zero-crossing point (lying the branch depicted in red). The second situation corresponds to a box B that contains two different x-critical points of C among which one is a zero-crossing point. According to the previous theorem, given a box that contains an x-critical zero-crossing point (0, y 0 , z 0 ) and no other x-critical points of C, it is sufficient to refine the latter until all the intersections with C lie on the sides that are orthogonal to the x-axis, in order to get an isolating box for (0, y 0 , z 0 ). This is always possible since we have assumed that there is no branch of C lying in a plan orthogonal to the x-axis.

Computational aspects

Let (0, y 0 , z 0 ) be an x-critical zero-crossing point of C and suppose given a 3D box B that isolate the point (0, y 0 , z 0 ) from other zero-crossing points of C.

A first computational step consists in refining the box B until (0, y 0 , z 0 ) is the only x-critical point of C that lies in the refined box. To do so, a straighforward approach consists in computing disjoint boxes that isolate the real zeros of the following system

S x := g(x, y, z) = 0 h(x, y, z) = 0 d(x, y, z) = 0
As for the system in (3), this can be done by computing a univariate representation of the solutions followed by a numerical approximation. This provides a set of disjoint 3D boxes that isolate the real zeros of S x . Thenceforth, it suffices to refine the box B until it overlaps with only one box of this set to get a refined box that contains only the x-critical point (0, y 0 , z 0 ). Let denote to this new box by B .

The second computational step consists in refining the box B until all the intersections with C lie only on the desired sides. Let B be a box resulting from the refinement of B up to a certain precision. We denote by y min , y max (resp. z min , z max ), the rational endpoints of its interval on y (resp. z). Then, the box B admits intersections only on the sides that are orthogonal to the x-axis if and only if the following algebraic systems do not have real zeros.

g(x, y min , z) = 0 h(x, y min , z) = 0 g(x, y max , z) = 0 h(x, y max , z) = 0 g(x, y, z min ) = 0 h(x, y, z min ) = 0 g(x, y, z max ) = 0 h(x, y, z max ) = 0 Consequently, to compute the desired box, our algorithm consists in refining successively the initial box and checking the existence of real solutions for the above algebraic systems. The algorithm stops when no real solution is found.

CONCLUSION

We addressed in this paper some aspects related to the study of the stability of differential time-delay systems with commensurate delays. More precisely, following the approach in common use, based on the computation of the critical pairs, we proposed a new strategy for analyzing the asymptotic behavior at these points, that gets rid of Puiseux series computations. Instead, our strategy proceeds by locally approximating at the critical pairs, the quasi-polynomial by a polynomial, and then by analyzing the asymptotic behavior at the zeros of this polynomial seen as specific points of algebraic curve in R 3 . The implementation of this approach is the subject of an ongoing work. In future work, an interesting problem is to investigate more deeply the impact of numerical errors, due to the truncation as well as the local approximation at the critical pairs, on the certification of our result.

Fig. 1 .

 1 Fig. 1. Two zero-crossing points