Yacine Bouzidi
email: yacine.bouzidi@inria.fr

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

Alban Quadrat
email: alban.quadrat@inria.fr

Computer algebra methods for the stability analysis of differential systems with commensurate time-delays

. Explicit examples are given to illustrate our algorithms.

I. INTRODUCTION

This paper aims at studying the stability of the class of linear time-invariant (LTI) differential systems with commensurate time-delays by means of computer algebra methods recently developed by the symbolic computation community.

An example of such a system is defined by the statespace representation ẋ(t) = m k=0 A k x(t -k τ ), where τ ∈ R + := {τ ∈ R | τ ≥ 0}, m ∈ N + := {0, 1, . . .}, A 0 , . . . , A m , and B 1 , . . . , B m are constant matrices with entries in a field K (e.g., K = Q, R). The characteristic function of this system is defined by a quasi-polynomial f (s, τ ) = mn j=0 p j (s) e -j τ s , where the p j 's are polynomials in the complex variable s with coefficients in K.

In this paper, we investigate the asymptotic stability of LTI differential commensurate time-delay systems whose dynamics are defined by quasipolynomials. Recall that 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}. Our approach for analyzing the asymptotic stability of this class of systems is based on the so-called critical pairs of f (s, τ ), that is the pairs (ω ∈ R, τ ∈ R + ) such that f (i ω, τ ) = 0. For instance, see, e.g., [START_REF] Li | Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays[END_REF], [START_REF] Gu | Stability of Time-Delay Systems[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 asymptotic behavior of the coordinates s of these pairs, called critical imaginary roots of f (s, τ ), that is the way these critical imaginary roots behave under small variation of the time-delay τ . There are two independent steps to analyze the stability. First, the critical pairs of f (s, τ ) are computed and then for each critical pair (s 0 , τ 0 ), the asymptotic behavior of the root s 0 with respect to small variation of τ 0 has to be studied. See [START_REF] Li | Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays[END_REF], [START_REF] Gu | Stability of Time-Delay Systems[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.

There exist several methods for computing the critical pairs. For the study of the asymptotic behavior of critical imaginary roots, several methods were developed for the case of simple imaginary roots (see [START_REF] Li | Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays[END_REF], [START_REF] Gu | Stability of Time-Delay Systems[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, a recent method based on Puiseux series was developed in [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].

The contributions of this paper are twofold. We first present a new approach for the computation of the critical pairs. Applying some transformations, this approach roughly reduces the computation of the critical pairs of a quasipolynomial to the computation of the real solutions of a zero-dimensional polynomial system in two variables, i.e., a system admitting a finite number of complex solutions. This can be achieved using the so-called Rational Univariate Representation (RUR) [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF] which is, roughly speaking, a oneto-one mapping between the solutions of the polynomial system and the roots of a univariate polynomial.

Our second contribution concerns the study of the asymptotic behavior of the quasipolynomial f (s, τ ) at a critical pair (s 0 , τ 0 ). More precisely, using recent advances obtained in the computer algebra community, we present an efficient algorithm that allows us to compute the different terms of Puiseux series. Hence, we can study the variation ∆ s in terms of ∆ τ in a neighborhood of a critical pair (s 0 , τ 0 ).

The paper is organized as follows. In Section II, we present two efficient and certified algorithms which compute the critical pairs. An algorithm which computes Puiseux series expansions is given in Section III.

II. COMPUTATION OF THE CRITICAL PAIRS

In this section, we focus on the certified computation of the critical pairs of the quasipolynomial f (s, τ ), that is its zeros of the form (i ω, τ ), where (ω, τ ) ∈ R×R + . This problem is equivalent to computing the set of real zeros of f (i ω, τ ). Due to the presence of transcendental terms, this quasipolynomial usually admits an infinite number of zeros. By means of particular transformations, this problem can be reduced to the computation of the real solutions of a zero-dimensional polynomial system, i.e., a polynomial system admitting only a finite number of complex solutions. We present two kinds of transformations. The first one, the Rekasius transformation, appears in [START_REF] Rekasius | A stability test for systems with delays[END_REF] and was used in a series of works [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]. The second one is a Möbius transformation. Rekasius transformation. This transformation consists in replacing the terms e -τ i ω in the quasipolynomial f (i ω, τ ) by the rational fraction1-T i ω 1+T i ω , where T ∈ R. Cleaning the denominators, we obtain a polynomial of the form R(ω, T )+ i I(ω, T ). Computing the critical pairs then amounts to studying the real zeros of the following polynomial system: R(ω, T ) = 0, I(ω, T ) = 0.

(

) 1
Given a solution (ω, T ) ∈ R 2 of (1), the critical delays can be obtained by τ k = 2 ω (arctan(ω T ) + k π) for k ∈ Z. Möbius transformation. An alternative to Rekasius transformation consists in replacing the terms e -τ i ω by a new variable z. When the delay τ varies, e -τ i ω covers the complex torus T := {z ∈ C | |z| = 1}. The problem of studying the zeros of f (i ω, τ ) then amounts to studying the zeros of the polynomial f (i ω, z), where ω ∈ R and z ∈ T. For more details, see, e.g., [START_REF] Li | Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays[END_REF] and the references therein.

A first approach is to consider z = u+i v, where u, v ∈ R, so that f (i ω, u

+ i v) = R(ω, u, v) + i I(ω, u, v)
, where R and I are two polynomials with real coefficients. The problem then leads to the problem of computing the real solutions (ω, u, v) of the following polynomial system: R(ω, u, v) = 0, I(ω, u, v) = 0, u 2 + v 2 -1 = 0. Generically, the above system is zero-dimensional, i.e., it only admits a finite number of complex solutions. Standard computer algebra methods such that the Rational Univariate Representation (RUR) [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF] can then be used to obtain certified numerical approximations of the real solutions.

In order to reduce the number of indeterminates in the above equivalent polynomial system, we propose the Möbius transform x-i x+i which maps R ∪ {∞} to T. Similarly as above, substituting e -τ i ω by x-i x+i into f (i ω, τ ) and cleaning the denominators, we obtain a polynomial of the form R(ω, x) + i I(ω, x), where R and I are two polynomials with real coefficients. The critical pairs can then be computed by solving the following polynomial system

R(ω, x) = 0, I(ω, x) = 0, (2)
and then the critical delays are obtained as follows:

τ k = 1 ω arctan 2 x x 2 -1 + k π , k ∈ Z. ( 3
)
Remark 1: We can observe that z = 1 is sent to ∞ in the above Möbius transformation. This case, which corresponds to the delays of the form 2 k π, where k ∈ Z, can be independently studied. The corresponding quasipolynomial is then a polynomial f (i ω, 2 k π) in ω whose real roots can be isolated using, e.g., classical bisection algorithms [START_REF] Collins | Polynomial real root isolation using descarte's rule of signs[END_REF], [START_REF] Rouillier | Efficient isolation of polynomial real roots[END_REF].

Up to a sign, the Rekasius transformation is equivalent to the Möbius transformation by setting x = -T ω.

For the above transformations, the computation of the critical pairs can be reduced to solving [START_REF] Li | On computing puiseux series for multiple imaginary characteristic roots of LTI systems with commensurate delays[END_REF] or [START_REF] Li | Analytic Curve Frequency-Sweeping Stability Tests for Systems with Commensurate Delays[END_REF].

In what follows, we only consider the case of the Möbius transformation (the approach being similar for Rekasius transformation). Let us now consider the polynomial system (2) resulting from the above transformation. Without loss of generality, we can assume that (2) admits only a finite number of complex solutions, which means that the polynomials R and I do not have a non-trivial common factor. Our objective is to study the solutions of (2) by means of modern computer algebra methods.

A first important remark is that the ω-coordinates of the solutions of f (i ω, τ ) = 0 (called the imaginary roots of f ), are actually the roots of a resultant 1 of R and I with respect to the variable x [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF].

Theorem 1: The roots of Res x (R, I) are the projection onto the ω-axis of the common solutions of R and I and the common roots of a n and b m .

If we want to check whether or not (2) admits imaginary roots, we can first compute Res x (R, I) and then check whether or not Res x (R, I) admits real roots. A method to achieve this is, for instance, to use the Descartes rule of sign [START_REF] Rouillier | Efficient isolation of polynomial real roots[END_REF]. If f (i ω, τ ) does not admit solutions with real ωcoordinates, then it means that the stability of the system is independent from the delay τ . In that case, a method for checking the stability consists in substituting, for instance, τ = 0 in f (s, τ ) and then applying the standard Routh criterion to the univariate polynomial f (s, 0).

From the end-user point of view, if K = R or C, then a convenient way to express the solutions V ( R, I 2) is to use a RUR [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF], that is a one-to-one mapping

) := (ω, x) ∈ K 2 | R(ω, x) = I(ω, x) = 0 of (
V ( R, I ) -→ V ( f h ) (ω, x) -→ ξ, g h,ω (ξ) g h (ξ) , g h,x (ξ) g h (ξ) ←-ξ, (4)
between

V ( R, I ) and V ( f h ) := {t ∈ K | f h (t) = 0}
for a certain univariate polynomial f h . In order to achieve the one-to-one condition, the representation ( 4) is computed with respect to a separating linear form h

:= a 1 ω + a 2 x ∈ Q[ω, x],
for certain a 1 , a 2 ∈ Q, that takes different values when evaluated at the different points of V ( R, I ). Using (4), the solutions of V ( R, I ) are then defined by

f h (t) = 0, ω = g h,ω (t) g h (t) , x = g h,x (t) g h (t) , (5)
where f h , g h , g h,ω , g h,x ∈ Q[t] and f h and g h are coprime polynomials, i.e., gcd(f h , g h ) = 1.

Computing a RUR requires solving the two problems:

• Find a separating linear form h := a 1 ω + a 2 x.

• Given a linear form h ∈ Q[ω, x], compute a RURcandidate, that is the polynomials f h , g h , g h,ω and g h,x . Computation of the RUR-candidate. Given a polynomial h in Q[ω, x] (not necessarly separating), the RUR-candidate with respect to h can be computed using the algorithm given in [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF]. This algorithm requires the knowledge of a basis of the finite Q-vector space A := Q[ω, x]/ R, I and a reduction algorithm which computes normal forms modulo the ideal R, I . In order to explicitly characterize the polynomials appearing in the RUR-candidate, we first define the Q-endomorphism defined by the multiplication by h ∈ Q[ω, x] in A, i.e., we consider m h : A -→ A defined by m h (p) = h p, where p denotes the residue class of p ∈ Q[ω, x] in A (i.e., modulo R, I ). A representative of p is the normal form of p with respect to the Gröbner basis of R, I . For more details, see [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF].

Given a basis {e 1 , . . . , e n } of A (which can be deduced from a Gröbner basis for any monomial ordering [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF]), we can compute the (n × n)-matrix M h associated to m h . The polynomial f h of the RUR-candidate is defined as the characteristic polynomial of the matrix M h . Moreover, if

f h := f h gcd f h , df h dt = d i=0 v i t d-i ∈ Q[t] is the square- free part of f h , and if we note H j := j i=0 v i t j-i ∈ Q[t]
for j = 0, . . . , d -1, then, we can define:

                       g h := d-1 i=0 Trace(M i h ) H d-i-1 , g h,ω := d-1 i=0 Trace(M ω M i h ) H d-i-1 , g h,x := d-1 i=0 Trace(M x M i h ) H d-i-1 .
Finding of a separating linear form. For the computation of a separating polynomial of V ( R, I ), a critical remark is that the number of non separating elements is bounded by m := n (n -1)/2, where n denotes the cardinal of V ( R, I ), that is the number of lines passing by two distinct points of V ( R, I ). Thus, a separating form can always be found among the set {h := ω + a x | a = 0, . . . , m}. On the other hand, the Bezout theorem [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF], [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF], states that for polynomials R and I of total degree d, the cardinal of V ( R, I ) is bounded by d 2 . Hence, a strategy for computing a separating element for V ( R, I ) is to loop over m := d 2 (d 2 -1)/2 + 1 different integers a, for each a, compute the number of distinct roots of the polynomial f h (see above), i.e., the degree of its squarefree part f h , and finally select an a for which this number is maximal. This ensures that the degree of f h is equal to the cardinal of V ( R, I ), and thus that the roots of f h are in bijection with the points of V ( R, I ). This strategy is actually time-consuming since it requires the computation of m charateristic polynomials as well as their squarfree parts. In practice, noticing that an arbitrary chosen linear form is separating with high probability, one prefers a strategy that consists in choosing randomly a linear form and testing that the latter is separating a posteriori of the computation of a RUR-candidate. This test is based on the fact that a linear form h is separating for V ( R, I ) if and only if the polynomials u ω := g h ω -g h,ω and u x := g h x-g h,x vanish on all the points of V ( R, I ) (i.e., belong to the radical of R, I ). In computational terms, this condition can be translated into (Trace(M uω ei )) i=1,..,n = (0, . . . , 0) and (Trace(M ux ei )) i=1,..,n = (0, . . . , 0) [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF].

Efficient RUR algorithm. Actually, the above RUR algorithm is designed for general zero-dimensional systems and is thus not optimized for the specific case of two polynomials in two variables. The computation of a Gröbner basis of {R, C}, required for the computation of the RUR-candidate as well as for the separating form, is often time-consuming.

Alternatively, below, we propose a practically efficient method that computes a decomposition into univariate representations of the solutions of (2) and which avoids the computation of a Gröbner basis. This method consists in shearing the system, i.e., applying a change of variables and then using resultants and subresultants polynomials [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF] to compute parametrizations of the solutions. Providing that the change of variables transforms the system into a generic position, the computed parametrization then encodes the solutions of the system in a one-to-one correspondence with the roots of a univariate polynomial.

Let us start with the following result which shows that there exist h := ω + a x such that, up to a factor in Q, the polynomial f h is equal to the resultant of two polynomials resulting from R and I after a change of variables.

Theorem 2 ( [START_REF] Bouzidi | Separating linear forms and rational univariate representations of bivariate systems[END_REF]): Let R(ω, x),

I(ω, x) ∈ Q[ω, x]. De- fine R (t, x) := R(t -a x, x) and I (t, x) := I(t -a x, x),
where a ∈ Z is such that the leading coefficient of R and I with respect to x are coprime. Then, the resultant of R and I with respect to x, denoted Res x (R , I ), is equal to:

c (α1, α2)∈V ( R, I ) (t -α 1 -a α 2 ) µ(α1,α2) , c ∈ Q.
Given a linear form ω + a x, it is well-known that the latter is separating for V ( R, I ) if and only if the system {R , I } is in generic position, i.e., for each root α of Res x (R , I ) (where R and I are defined as in Theorem 2), the gcd of R (α, x) and I (α, x), denoted G(α, x), has exactly one root (i.e., there exists only one intersection point above each t = α).

To check the above genericity condition, we first perform a triangular decomposition of {R , I }. More precisely, we compute a set of triangular systems of the form {r k (t), Sres k (t, x)} k=1,...,l where l is a positive integer bounded by the minimum degree in x of R and I , and such that l k=1 r k (t) = Res x (R , I ), and Sres k (α, x) of degree k in x is equal to G(α, x) for any root α of r k (t). The polynomial Sres k (t, x) = k i=0 sr k,i x i is known as the k-th subresultant polynomial of R and I with respect to x (see [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]). We refer to [START_REF] Bouzidi | Separating linear forms and rational univariate representations of bivariate systems[END_REF] for more details.

Once a triangular decomposition is computed, the genericity condition is equivalent to the fact that Sres k (t, x) writes modulo r k as (a k (t) x -b k (t)) k with gcd(a k , b k ) = 1. To check the latter condition, we can use the following result of which the proof can be found in [START_REF] Daouda | On the computation of the topology of a non-reduced implicit space curve[END_REF].

Theorem 3: Let R(ω, x), I(ω, x) ∈ Q[ω, x].
Define the polynomials R (t, x), I (t, x) as in Theorem 2, and let {r k (t), Sres k (t, x)} k=1,...,l be the triangular decomposition of {R , I }. Then, ω + a x is a separating element for V ( R, I ) if and only if we have

k (k -i) sr k,i sr k,k -(i + 1) sr k,k-1 sr k,i+1 = 0 mod r k ,
for all k ∈ {1, . . . , l} and for all i ∈ {0, . . . , k -1}.

Provided that the conditions of Theorem 3 are satisfied (i.e., the system is in generic position), each

Sres k (t, x) is then of the form (a k (t) x -b k (t)) k with gcd(a k , b k ) = 1.
Moreover, we can obtain the explicit expression of a k and b k by computing the k-th derivative polynomial of Sres k (t, x), which yields that a k = k! sr k,k and b k = (k-1)! sr k,k-1 . Finally, the solutions of the system {R, I} can be parametrized using the following set of parametrizations:

               r k (t) = 0, x = -sr k,k-1 (t) k sr k,k (t) , ω = t -a sr k,k-1 (t) k sr k,k (t) , k = 1, . . . , l. (6)
Numerical approximations. Once a RUR of the solutions of ( 2) is computed, we can obtain certified numerical approximations of these solutions by first isolating the real roots of f h by means of intervals using, for instance, the algorithm in [START_REF] Rouillier | Efficient isolation of polynomial real roots[END_REF], and then substituting these intervals into the rational fractions

g h,ω (t)
g h (t) and g h,x (t) g h (t) in order to get isolating intervals for the coordinates ω and x of the solutions. Moreover, substituting these intervals into (3) yields intervals for the delays corresponding to each solution.

Efficient algorithms for computing the RUR (5) and the parametrization (6) are implemented in Maple (see the command Rational Univariate Representation of the package Groebner for the first one, and the external library RS 2 for the second). Moreover, the numerical approximations of the solutions, obtained by means of a RUR, can be computed using the Maple command Isolate of the package RootFinding with the option Method="RS" and Output="interval".

Example 1: We consider the quasipolynomial f (s, τ ) = e -3 τ s -3 e -2 τ s + 3 e -τ s + s 4 + 2 s 2 [2, example 4.6].

To compute the critical pairs, we consider the following system, which results from the substitutions s = i ω and e

-τ i ω = x-i x+i into f (s, τ ) R(ω, x) = x 3 -3 x ω 4 + -2 x 3 + 6 x ω 2 + x 3 -3 x = 0, I(ω, x) = 3 x 2 -1 ω 4 + -6 x 2 + 2 ω 2 + 3 x 2 + 7 = 0,
and the f (i ω, 2 k π) = (ω 2 -1) 2 (see Remark 1). The latter polynomial yields the critical pairs (i, 2 k π) and (-i, 2 k π) where k ∈ Z. Note that the critical imaginary roots have multiplicity two.

2 https://who.rocq.inria.fr/Fabrice.Rouillier/software.php Computing a RUR of the former system with respect to the separating linear form h = ω + x, we obtain: which, after substitution into the rational fractions of the RUR, yields x 1 = x 2 = 0 and thus ω j = t j for j = 1, 2.

         f h = t 4 -2 t 2 -
Note that the two ω-coordinates are roots of the polynomial t 4 -2 t 2 -7 which can be solved symbolically. We then get ω j = (-1) j 1 + 2 √ 2 for j = 1, 2. Finally, the critical delays τ j,k can then be obtained by [START_REF] Gu | Stability of Time-Delay Systems[END_REF]. In our case, the solutions (ω 1 , x 1 ) and (ω 2 , x 2 ) yields to τ j,k = k π ωj , where k ∈ Z and j = 1, 2. Alternatively, we can compute [START_REF] Olgac | An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems[END_REF]. Since all these solutions of the system are simple (i.e., their multiplicities are equal to 1), this yields the following single parametrization: The two algorithms, based respectively on the computation of ( 5) and ( 6), were implemented in Maple. Below, we report their running times (on a laptop with an Intel i-7 processor and L8 cache) for randomly generated quasipolynomials of total degree deg and density dens (sparse/dense).

     r 1 = t 4 -2 t 2 -
deg dens numsol timing for [START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF] In the table, timings are the running times in CPU seconds and numsol is the average number of critical pairs (ω, x).

III. COMPUTING PUISEUX SERIES

This section is organized as follows. We first provide definitions and algorithms for the case of polynomials. This situation is well-studied in the computer algebra literature. See [START_REF] Duval | Rational Puiseux expansions[END_REF], [START_REF] Poteaux | Calcul de développements de Puiseux et application au calcul de groupe de monodromie d'une courbe algébrique plane[END_REF], [START_REF] Poteaux | Complexity bounds for the rational newton-puiseux algorithm over finite fields[END_REF] for instance. We then explain how these results can be adapted to the case of quasipolynomials. As shown in [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 results can be used to effectively decide the stability of a LTI differential system with commensurate delays, which are multiples of τ , when τ changes. Definition (polynomial case). Let us consider a bivariate polynomial F ∈ K[X, Y ], where K is a field of characteristic 0 (e.g., K = Q, R, C). Let K denote the algebraic closure of K. If we consider the roots of F viewed as a univariate polynomial in Y , then we have the Puiseux theorem.

S ij (X) = ∞ k=ni α ik ζ jk ei (X -x 0 ) k e i , 1 ≤ i ≤ s,
where 0 ≤ j ≤ e i -1, n i ∈ Z, ζ ei is the e i -th root of unity, and

α ini = 0 if S ij = 0, such that F (X, S ij (X)) = 0.
Over most points x 0 , the polynomial F (x 0 , Y ) has d distinct roots. In such a case, the Implicit Function Theorem ensures that the series solutions are actually Taylor series. Such series can be quickly computed using, e.g., quadratic Newton iterations [START_REF] Kung | All algebraic functions can be computed fast[END_REF]. When dealing with multiple roots over a point x 0 , this may not be the case anymore.

Computing the first term of the Puiseux series. We will compute the Puiseux series using a variant of the well-known Newton-Puiseux algorithm [START_REF] Walker | Algebraic Curves[END_REF], [START_REF] Brieskorn | Plane Algebraic Curves[END_REF], namely the rational Newton-Puiseux algorithm [START_REF] Duval | Rational Puiseux expansions[END_REF]. Let us explain the main idea of this algorithm by means of an example, beginning with the computation of the first term of each Puiseux series.

Example 2: Let us compute the Puiseux series at

x 0 = 0 of F := Y 6 +Y 5 X+5 Y 4 X 3 -2 Y 4 X+4 Y 2 X 2 +X 5 -3 X 4 .
Note that we can always be in this situation via a change of the variable X ← X + x 0 . From Theorem 1, we know that the first term of any such series S(X) is of the form α X m q , α ∈ K and (m, q) ∈ N 2 . We then have:

F (X, α X m q + • • • ) = α 6 X 6 m q + α 5 X 5 m q +1 + 5 α 4 X 4 m q +3 -2 α 4 X 4 m q +1 + 4 α 2 X 2 m q +2 + X 5 -3 X 4 + • • •
To get F (X, Y (X)) = 0, at least two terms of the above sum must cancel, i.e. (m, q) must be chosen so that two of the above exponents coincide. This leads us to consider the support of a polynomial F .

Definition 1: Given F (X, Y ) = i,j a ij X j Y i , then Supp(F ) = {(i, j) ∈ N 2 | a ij = 0} is the support of F .
The condition on (m, q) can be translated as follows: two points of Supp(F ) belong to the same line m i + q j = l. In order to increase the X-order of the evaluation, there must be no point under this line.

Example 3: Considering again Example 2, we have Supp(F ) = {(0, 6), [START_REF] Li | On computing puiseux series for multiple imaginary characteristic roots of LTI systems with commensurate delays[END_REF][START_REF] Niculescu | Delay Effects on Stability: A Robust Control Approach[END_REF], [START_REF] Gu | Stability of Time-Delay Systems[END_REF][START_REF] Marshall | Time-Delay. Systems: Stability and Performance Criteria with Applications[END_REF], [START_REF] Li | On computing puiseux series for multiple imaginary characteristic roots of LTI systems with commensurate delays[END_REF][START_REF] Marshall | Time-Delay. Systems: Stability and Performance Criteria with Applications[END_REF], (2, 2), (5, 0), (4, 1)}, and we have two such lines i + 2 j = 6, and i + j = 4. They define the so-called Newton polygon of F . Definition 2: The Newton polygon N (F ) of F ∈ K[X, Y ] is the lower part of the convex hull of its support.

Example 4: Let us consider again Examples 2 and 3. We consider α corresponding to i + 2 j = 6 and we get F (T 2 , α T ) = (α 6 -2 α 4 + 4 α 2 ) T 6 -3 T 8 + α 5 T 7 + (5 α 4 + 1) T 10 + . . . meaning that α must be a root of P = Z 6 -2 Z 4 + 4 Z 2 . One can notice that P ∈ K[Z 2 ] and that 0 is a root of P . As we are not interested in such a root (we will get the first non-zero term of the corresponding Puiseux series by considering the other edge of the Newton polygon), we can see that it is more interesting to consider a root of the polynomial φ = T 2 -2 T + 4.

By construction, the polynomial P will always belong to K[Z q ]. The polynomial φ defined in Example 4 is a characteristic polynomial.

Definition 3:

If H = a ij X j Y i K[X, Y ], then the char- acteristic polynomial φ ∆ of ∆ ∈ N (H) is defined by φ ∆ (T ) = (i,j)∈∆ a ij T i-i 0 q
, where i 0 is the smallest value such that (i 0 , j 0 ) belongs to ∆ for some j 0 . Computing more terms. Up to now, we have described a way to compute the first non zero term of each Puiseux series: we compute the Newton polygon, the characteristic polynomial of each of its edges, and compute the roots of these polynomials. Note first that we assume the use of symbolic computation techniques: we have to decide whether or not coefficients are 0 to be able to certify the correctness of our computations. At the end of this section, we will briefly discuss a symbolic-numeric strategy developed by the second author to overcome possibly costly symbolic computations.

In our context (assuming the generalisation in the next paragraph), if one of the roots of the characteristic polynomial leads to a purely imaginary coefficient for one of the Puiseux series (in the sequel, such a root is denoted by ξ), we then have to compute its next term to decide on the stability of the system. This can actually be done by applying the above strategy to F X, Y + ξ

1 q X m q ∈ K ξ 1 q X 1 q , Y , or similarly to F (X q , Y +ξ 1 q X m ) ∈ K ξ 1 q [X, Y ].
We can also avoid taking a q-th root of ξ by considering polynomial

F (ξ v X q , X m (Y + ξ u ))/X l ∈ K(ξ)[X, Y ], where u, v ∈ Z are such that u q -m v = 1 [15, section 4].
Factor(K, φ) Input: K a field, and φ

∈ K[T ]. Output: A set {(φ i , M i )} i , where φ i ∈ K[T ] is monic and irreducible, and φ = c i φ Mi i for some c ∈ K. RNPuiseux(F, K) Input: K a field, F ∈ K[X, Y ] separable. Output: Truncated Puiseux series of F . 1. R ← {} 2. For ∆ in N (F ) do 3. Compute q, m, l, and φ ∆ associated to ∆ 4. Compute (u, v) ∈ Z 2 s. t. u q -m v = 1 5. For (φ, M ) in Factor(K, φ) do 6. If M = 1 then 7. R ← R ∪ [ξ v X q , ξ u X m ] 8. Else 9.
Let ξ be any root of φ 10.

G ← F (ξ v X q , X m (ξ u + Y ))/X l 11.

For (P, Q) in RNPuiseux(K(ξ), G) do 12.

R ← R ∪ {(ξ v P q , P m (ξ u + Q))} 13. Return R.

The quasipolynomial case. As mentioned in [START_REF] Li | On computing puiseux series for multiple imaginary characteristic roots of LTI systems with commensurate delays[END_REF], we can approximate the quasipolynomial by truncated power series.

In [START_REF] Poteaux | Calcul de développements de Puiseux et application au calcul de groupe de monodromie d'une courbe algébrique plane[END_REF], [START_REF] Poteaux | Complexity bounds for the rational newton-puiseux algorithm over finite fields[END_REF], it is shown that truncation bounds exist to ensure the exactness of the Puiseux series. A first answer is therefore to compute such an approximation and then apply the standard algorithm explained above.

For the computation of the first term another way, described in [START_REF] Li | On computing puiseux series for multiple imaginary characteristic roots of LTI systems with commensurate delays[END_REF], is to replace any coefficient access in the polynomial case by an evaluation of some derivatives of the quasipolynomial. For instance, if f (s, τ ) is a quasipolynomial, then 1 2! 1! ∂3 f ∂ 2 s ∂τ (τ 0 , s 0 ) corresponds to the coefficient of the term (s -s 0 ) 2 (τ -τ 0 ).

For higher order terms, we can apply a similar idea: instead of evaluating the derivatives at a point, we evaluate it at the truncated Puiseux series already computed. Details about this strategy will be developed in a forthcoming work.

Finally, we have to change the condition in line 6 in RNPuiseux. The aim of RNPuiseux is to desingularise the curve F (X, Y ) = 0, i.e. to stop when the multiplicity of the root ξ is 1. In our context, we might have to stop before such a condition is fulfilled or might have to consider more terms. We stop when the computed coefficient has a non-zero real part. A symbolic-numeric strategy. In [START_REF] Poteaux | Computing Monodromy Groups Defined by Plane Algebraic Curves[END_REF], [START_REF] Poteaux | Calcul de développements de Puiseux et application au calcul de groupe de monodromie d'une courbe algébrique plane[END_REF], a symbolicnumeric strategy was developed to compute a numerical approximation of the Puiseux series coefficients with correct exponents. Roughly speaking, the idea is to first compute the Puiseux series modulo a well-chosen prime number, which gives the structure (exponents. . . ) of the series, and then to use this structure to conduct purely numeric operations. Adapting this strategy to the quasi-polynomial case is an interesting challenge which will be studied in a future work. Some examples. We implemented a Maple prototype to compute the Puiseux series of a quasipolynomial as explained above. We illustrate the results with examples and we show the logs of some computations (including the timing).

Example 5: We consider again the quasipolynomial defined in Example 1 and compute its Puiseux series at (τ 0 , s 0 ) = (2 π, i). We get the following: 2 , ∆s = -0.0987 i (∆τ ) 1/3 -(0.03557 -0.0028 i) (∆τ ) 2/3 . In this last example, the coefficient growth is important. For instance, the root ξ that is not written above is:

Theorem 1 :

 1 Let F ∈ K[X, Y ] and d the degree of F in Y . Given x 0 ∈ K, there exist positive integers e 1 , . . . , e s with s i=1 e i = d and d distinct series

  Real part of the coefficient is 1/4 * 2ˆ(1/2) It took .88e-1 secondsHence, we obtain ∆s = (.3535 + .3535 i) (∆τ )

3 2 . 6 :

 326 Example We consider the quasipolynomial F (τ, s) = 4 k=0 a k (s) e -k s τ where a 4 (s) = s 3 + 2 s 2 + 2 s + 1, a 3 (s) = 3 s 3 + 9 s 2 + 9 s + 4, a 2 (s) = 5/4 π s 5 + 11/4 π s 4 + 3 s 3 -π s 3 + 1/2 s 2 π + 13 s 2 + 15 s -9/4 s π + 6 -9/4 π, a 1 (s) = 5/4 π s 5 + 11/2 π s 4 + s 3 + 7/2 π s 3 + s 2 π + 7 s 2 + 11 s + 9/4 s π + 4 -9/2 π, and a 0 (s) = 1 + 3 s + 9/2 π s 3 + s 2 -45/8 π 2 +9/2 s π -15/8 s 4 π 2 +1/2 s 2 π -75/8 s 2 π 2 -9/4 π + 11/4 π s 4 + 15/8 π 2 s 6 (this polynomial is given in http://web1.lss.supelec.fr/delsys2015/Site/ Program_files/Wed7_LI.pdf).Computing the Puiseux series at the point (5 π, i), we get:Newton polygon is [[0 2] [1 1] [4 0]] Edge 1, xi is -((1/2) * I)/Pi Real part is 0; Going for a recursive call Newton polygon is [[0 1] [1 0]] xi is (-3/208+(1/104) * I) * (117 * Pi+2-3 * I)/Piˆ2 Real part is -(27/16)/Pi, We're done. Edge 2, xi is ((6/5) * I)/(50 * Piˆ3-30 * Piˆ2-3 * Pi) Real part is 0; Going for a recursive call Newton polygon is [[0 1] [1 0]] xi is (see below) real part is non zero. We're done. It took 3.596 seconds This leads to the two following Puiseux series: ∆s = -0.1591 i ∆τ -(0.5371 -0.3644 i) (∆τ )

  7 t 8 -16 t 6 + 74 t 4 -32 t 2 + 25 , g h = t t 10 -15 t 8 + 66 t 6 -34 t 4 -143 t 2 + 29 , g h,ω = t 10 -9 t 8 -26 t 6 + 234 t 4 + 53 t 2 + 35, g h,x = 2 t 4 -2 t 2 -7 t 6 -10 t 4 + 17 t 2 -10 . Isolating the roots of f h , we obtain the two real roots

	t 1 ∈ [-268918098581 137438953472 , -67229524645 34359738368 ], t 1 ≈ -1.956636687,
	t 2 ∈ [ 67229524645 34359738368 , 268918098581 137438953472 ],	t 2 ≈ 1.956636687,

  7 t 8 -16 t 6 + 74 t 4 -32 t 2 + 25 , sr 1,0 = 65536 t 3 t 6 -18 t 4 + 7 t 2 + 18 , sr 1,1 = -327680 t 6 + 1179648 t 4 + 983040 t 2 + 655360.

Given an integral domain A and two polynomialsf = n i=0 a i x i and g = m i=0 b i x i in A[x], the resultant of f and g w.r.t x, denoted by Resx(f, g), is the determinant of the so-called Sylvester matrix. The latter belongs to A.

25π 2 200 iπ 2 -180 iπ -750 π 3 -21 i + 600 π 2 125000 π 6 -225000 π 5 + 112500 π

-6750 π 2 -810 π -27 .

This explains why the running time is longer on this example. This advocates for the development of an efficient symbolic-numeric strategy.

Supported by the ANR-13-BS03-0005.


Guillaume Rance
email: guillaume.rance@safrangroup.com

Yacine Bouzidi
email: yacine.bouzidi@inria.fr

Alban Quadrat
email: alban.quadrat@inria.fr

Arnaud Quadrat

A symbolic-numeric method for the parametric H ∞ loop-shaping design problem

In this paper, we present a symbolic-numeric method for solving the H∞ loop-shaping design problem for low order single-input single-output systems with parameters. Due to the system parameters, no purely numerical algorithm can indeed solve the problem. Using Gröbner basis techniques and the Rational Univariate Representation of zero-dimensional algebraic varieties, we first give a parametrization of all the solutions of the two Algebraic Riccati Equations associated with the H∞ control problem. Then, following the works of [1], [9] on the spectral factorization problem, a certified symbolic-numeric algorithm is obtained for the computation of the positive definite solutions of these two Algebraic Riccati Equations. Finally, we present a certified symbolic-numeric algorithm which solves the H∞ loop-shaping design problem for the above class of systems. This algorithm is illustrated with a standard example.

I. INTRODUCTION

Nowadays, automatic control usually consists in the solving of complex problems. As a consequence, this is usually done numerically for a particular system and hardware performances are required to handle the underlying numerical computations: the more a system or a problem is refined, i.e., takes into consideration finer dynamics of the physical system or constraints, the more time-computation is required.

An alternative to this approach is to solve these problems symbolically, that is to say to solve the problems corresponding to a certain class of systems. The goal is to compute the solution of a control problem for a set of systems depicted by some parameters. Its solution (e.g., a stabilizing controller) then depends on the system parameters. In what follows, this approach will be called a "symbolic method". One of the greatest benefits of this approach is to simply obtain the solution for a particular system belonging to the class of systems by numerically evaluating the system parameters in the closed-form solution. Since only evaluation operations are required, this task can be achieved at a cheap computational cost.

The purpose of this paper is to develop a symbolic method for the H ∞ loop-shaping design problem [START_REF] Glover | Robust stabilization of normalized coprime factor plant descriptions with H∞-bounded uncertainty[END_REF], [START_REF] Vinnicombe | Uncertainty and Feedback -H∞ Loop-shaping and the ν-gap metric[END_REF], [START_REF] Zhou | Robust and Optimal Control[END_REF] for linear systems containing parameters. This synthesis problem yields the computation of a stabilizing controller which takes into account robustness objectives in the frequency Paris, Massy, France. arnaud.quadrat@safrangroup.com domain (e.g., stability margins) and the modeling of different transfer functions via perturbations. Practically, H ∞ control provides a natural compromise between the performance and the robustness of the control-loop system. The desired performance properties and robustness conditions can be obtained by "weighting" the plant transfer function, which results in a "weighted plant".

On a practical point of view, such an approach is interesting in the design stage of a project, when a designer wants to know if the global design of its system has a chance to reach the specifications, or if a global rework is needed. As the global design is likely to change at some stage, the engineer has to quickly provide a range of performances and robustness that the closed-loop system could achieve in order to choose the appropriate global design. This approach is also interesting in adaptive control where the controller depends on the measurement of variables or parameters (e.g., a mechanical mode) and has to be re-computed in real time (i.e., quickly relative to the dynamics of the system). As the technologies of embedded processors are not powerful enough to compute H ∞ optimization problems in real time, the symbolic approach developed in this paper is justified.

Solving the H ∞ control problem for a system with parameters requires solving an Algebraic Riccati Equation (ARE), whose resolution was studied in [START_REF] Forsman | Solving the ARE symbolically[END_REF] using Gröbner basis techniques. A method to solve an H ∞ control problem was proposed in [START_REF] Kanno | Symbolic-numeric hybrid optimization for plant/controller integrated design in H∞ loop-shaping design[END_REF] based on properties of the spectral factorization of the Hamiltonian matrix associated with the ARE [START_REF] Kanno | Solution of algebraic Riccati equations using the sum of roots[END_REF]. In this article, a different method is proposed by means of the study of the polynomial system defined by an ARE (Section III). A Rational Univariate Representation of the solutions of this polynomial system is given (Section IV). In Section V, the equivalence between our method and the one developed in [START_REF] Kanno | Symbolic-numeric hybrid optimization for plant/controller integrated design in H∞ loop-shaping design[END_REF] is shown. In Section VI, an algorithm is given to compute the optimal solution to the H ∞ control problem by means of certified root isolation methods. Finally, in Section VII, a standard example illustrates our method.

II. THE STANDARD H ∞ -CONTROL PROBLEM

In this paper, we will consider a single-input single-output (SISO) finite-dimensional linear system (Figure 1) defined by y 1 = G e 1 , where the transfer function G is given by

G := c n s n + c n-1 s n-1 + . . . + c 1 s + c 0 a n s n + a n-1 s n-1 + . . . + a 1 s + a 0 , (1)
a i , c i ∈ R for i = 0, . . . , n and n is the order of G (i.e., a n = 0). In what follows, to get a strictly proper transfer

u 1 + + e 1 G e 2 u 2 + - y 2 y 1 K Fig. 1. Control scheme
function G, we set:

a n := 1, c n := 0. (2)
Let us consider its controller canonical form defined by the following state-space representation:

ẋ = A x + B e 1 , y 1 = C x, (3)
A :=           0 1 0 . . . 0 0 0 1 . . . . . . . . . . . . . . . . . . 0 0 0 . . . 0 1 -a 0 -a 1 . . . -a n-2 -a n-1           ∈ R n×n , B := (0 . . . 0 1) T ∈ R n×1 , C := (c 0 . . . c n-1 ) ∈ R 1×n .
(4) Let R(s) be the field of rational functions with real coefficients. Let K ∈ R(s) be a controller and let us consider the closed-loop system defined in Figure 1. If we denote the sensitivity transfer function by S := (1 + G K)

-1 , then:

e 1 y 1 = S K S G S G K S u 1 u 2 .
Let us consider the following standard problem.

Robust Control Problem (RCP): Given γ > 0, find a controller K which stabilizes G (i.e., such that the rational transfer functions S, K S and G S are proper and stable) and is such that:

S K S G S G K S ∞ < γ. (5)
For more details, see [START_REF] Glover | Robust stabilization of normalized coprime factor plant descriptions with H∞-bounded uncertainty[END_REF], [START_REF] Zhou | Robust and Optimal Control[END_REF], [START_REF] Duc | Commande H∞ et µ-analyse, des outils pour la robustesse[END_REF], [START_REF] Vinnicombe | Uncertainty and Feedback -H∞ Loop-shaping and the ν-gap metric[END_REF] and the references therein. The RCP yields a compromise between the performance of the closed-loop system and the robustness with respect to the perturbations u 1 and u 2 .

We briefly state a standard result of H ∞ -control theory. Theorem 1: [8, Cor. 5.1], [START_REF] Zhou | Robust and Optimal Control[END_REF]Ch. 18], [6, Th. 3.2] Let (A, B, C) be a controllable and observable state-space representation (3) of the transfer function G defined by [START_REF] Anai | Parametric polynomial factorization using the sum of roots and its application to a control design problem[END_REF]. Then, the optimal value of γ of ( 5) is given by

γ opt = 1 + λ max (Y X),
where X and Y are respectively the unique real positive definite solutions of the following Algebraic Riccati Equations (ARE) [START_REF] Forsman | Solving the ARE symbolically[END_REF] and λ max is the greatest eigenvalue of Y X (all the eigenvalues of Y X are real positive).

R := X A + A T X -X B B T X + C T C = 0, (6) R := Y A T + A Y -Y C T C Y + B B T = 0,
For γ > γ opt , a controller K γ satisfying the RCP is defined by

ż = A γ z + B γ e 2 , y 2 = C γ z, (8)
with the following notations:

         Z γ := I n + Y X -γ 2 I n -1 , A γ := A -B B T X + γ 2 Z γ Y C T C, B γ := -γ 2 Z γ Y C T , C γ := B T X.
We note that R T = R and R T = R. For SISO systems, the next lemma shows that ( 7) is a consequence of (6). Lemma 1: We have Y := QXQ, where Q = Q T is defined by

Q -1 := P = P 1 • • • P n , (9)
where:

P i T := C n-i j=0 a n-j A n-i-j , i = 1, . . . , n.
Corollary 1: Let ρ r (M ) be the set of all the real eigenvalues of M ∈ R n×n and λ (M ) := max λ∈ρr(M ) |λ|. We have:

γ opt = 1 + λ max (Y X) = 1 + λ 2 (Q X).
(10) In this paper, we focus on the computation of X in the case where the a i 's and c j 's are unknown parameters and not fixed numerical values. In particular, numerical algorithms for the computation of the positive definite solutions of ARE cannot be used. We have to consider symbolic methods to take into account the parameters a i 's and c j 's.

III. POLYNOMIAL SYSTEM DESCRIPTION OF AN ARE

We first study the equation R = 0. For 1 ≤ i ≤ j ≤ n, x i,j denotes the entries of X and we set x i,0 = x 0,j = 0. From direct computations, we get the following proposition.

Proposition 1: The equation R = 0, defined by [START_REF] Duc | Commande H∞ et µ-analyse, des outils pour la robustesse[END_REF], is equivalent to the following system of n(n+1) 2 polynomial equations in the n(n+1) 2 unknowns {x i,j } 1≤i≤j≤n :

-(x i,n + a i-1 ) (x j,n + a j-1 ) +x i,j-1 + x i-1,j + a i-1 a j-1 + c i-1 c j-1 = 0, 1 ≤ i ≤ j ≤ n. (11)
Let us introduce the following notations:

k < 0, b k := 0, 0 ≤ k ≤ n -1, b k := x k+1,n + a k , k = n, b k := 1, k > n, b k := 0. (12)
Theorem 2: For k = 1 . . . n, we set x k,0 = x 0,k := 0, and for (i, j) ∈ N 2 , we define:

N (i, j) := i -1, 2 ≤ i + j ≤ n + 1, N (i, j) := n -j + 1, n + 1 < i + j ≤ 2 n + 1.
The elements of X solution of (6) are determined only by the b k 's as follows

x k,n = b k-1 -a k-1 , x i,j-1 = N (i,j) k=0 (-1) k b i-1-k b j-1+k -θ N (i,j) , (13)
where 1 ≤ k ≤ n, 1 ≤ i < j ≤ n, and θ m is defined by:

θ m := m k=0 (-1) k (a i-1-k a j-1+k + c i-1-k c j-1+k ) . ( 14
)
Proof: The diagonal elements of the matrix R yield:

1 ≤ i ≤ n, 2 x i-1,i = 2 x i,i-1 = b 2 i-1 -a 2 i-1 + c 2 i-1 . (15
) Using ( 11) and ( 12), the i th row of X can be deduced from the (i -1) th row, which allows us to write x i,j-1 in terms of x i-1,j which itself can be expressed in terms of x i-2,j+1 , . . . After m substitutions, for 1 ≤ i < j ≤ n -m, we obtain

x i,j-1 = (-1) m+1 x i-1-m,j+m + m k=0 (-1) k b i-1-k b j-1+k -θ m ,
where θ m is defined by [START_REF] Vinnicombe | Uncertainty and Feedback -H∞ Loop-shaping and the ν-gap metric[END_REF].

We distinguish two halting conditions in the induction: 1) We reach the first row of X, i.e. an element of the form x i-1-m,j+m = x 0,j+m = 0 for m = i -1.

In this case, note that we have 1 ≤ j + m ≤ n, i.e. 2 ≤ i + j ≤ n + 1. By reaching this limit, one obtains the following expression:

∀(i, j) ∈ N 2 , i < j, 2 ≤ i + j ≤ n + 1, x i,j-1 = i-1 k=0 (-1) k b i-1-k b j-1+k -θ i-1 . (16)
2) We reach the last column of X, i.e. an element of the form

x i-1-m,j+m = x i-1-m,n for m = n -j.
In this case, note that we have

1 ≤ i -1 -m ≤ n, i.e. n + 1 < i + j ≤ 2 n + 1. By reaching this limit, using x i+j-n-1,n = b i+j-n-2 b n -a i+j-n-2 a n (see (2)
and ( 12)), we obtain the following expression:

∀(i, j) ∈ N 2 , i < j, n + 1 < i + j ≤ 2 n + 1, x i,j-1 = n-j+1 k=0 (-1) k b i-1-k b j-1+k -θ n-j+1 .
(17) Hence, the entries x i,j of X can be explicitly expressed in terms of the b k 's defined in [START_REF] Rouillier | Algorithmes pour l'étude des solutions réelles des systèmes polynomiaux[END_REF], i.e., up to constant a k-1 , in terms of the elements of the last column x k,n of X.

We note that for i = 1, using x 1,0 = 0, (15) then yields:

B 0 := b 2 0 -(a 2 0 + c 2 0 ) = 0.
For i = 1, x i-1,i appears twice in ( 15) and ( 16) or (17). Thus, combining ( 15) and ( 13) for i = j = k + 1, we obtain a polynomial system of n -1 equations in n -1 unknowns

{b i } 1≤i≤n-1 . If M is defined by:      M (k) := k, 1 ≤ k ≤ n -1 2 , M (k) := n -k, n -1 2 < k ≤ n -1, (18)
then this system, denoted by B, is defined by

B := {B k = 0} 0≤k≤n-1 , (19)
where the polynomials B k are given by

           B 0 := b 2 0 -d 0 , B k := b 2 k + 2 M (k) m=1 (-1) m b k-m b k+m -d 2k , 1 ≤ k ≤ n -1, (20)
and the constants d 2k are defined by

               d 0 := a 2 0 + c 2 0 , d 2k := 2 M (k) m=1 (-1) m (a k-m a k+m + c k-m c k+m ) + a 2 k + c 2 k , 1 ≤ k ≤ n -1, d 2n := 1.
(21) Remark 1: We point out that the parameters a i 's and c j 's appear only in the constant terms of the B k 's (compare with [START_REF] Rouillier | Solving zero-dimensional systems through the Rational Univariate Representation[END_REF]), which highly simplifies the study of (19).

If the d 2k 's are fixed, then the polynomial system (21) defines all the transfer functions G (see [START_REF] Anai | Parametric polynomial factorization using the sum of roots and its application to a control design problem[END_REF]) which give the same solution of [START_REF] Duc | Commande H∞ et µ-analyse, des outils pour la robustesse[END_REF]. The study of this algebraic variety will be studied in details in a future publication.

Let K := Q(d 0 , .., d 2n ) be the field for rational functions in d 0 , . . . , d 2n , K its algebraic closure and K := Q(a 0 , . . . , a n-1 , c 0 , . . . , c n-1 ). Using Gröbner basis techniques (see, e.g., [START_REF] Cox | Varieties, and Algorithms Springer[END_REF] and the references therein), we can now state some results on the polynomial system (19).

Theorem 3: The polynomial system B given by ( 19) and defined over

R := K[b 0 , .., b n-1 ] is zero-dimensional (i.e., A := R/ B , where B := { n i=1 r k B k | r k ∈ R} is the ideal of R generated by B, is a finite-dimensional K- vector space).
The dimension of the K-vector space A is 2 n . Moreover, B is a reduced Gröbner basis of B with respect to the degree reverse lexicographic (DRL) order b n-1 . . . b 0 .

Proof: According to the Buchberger criterion (see, e.g., [START_REF] Cox | Varieties, and Algorithms Springer[END_REF] where α := (α 0 , . . . , α n-1 ) ∈ N n , a basis of A is defined by {b α } α∈ 0,1 n . Finally, the dimension of A is equal to 2 n . Remark 2: Combining Theorems 2 and 3, we get that the polynomial system R = 0 associated with an ARE of a SISO system (1) is a finite-dimensional K-vector space of dimension 2 n with basis {b α } α∈ 0,1 n (see [START_REF] Rouillier | Algorithmes pour l'étude des solutions réelles des systèmes polynomiaux[END_REF]).

IV. PARAMETRIZATION OF THE SOLUTIONS

In Section III, we show that we can reduce the problem of computing the solutions of the ARE equation R = 0 to those of a zero-dimensional system B. This system is given by its Gröbner basis B with respect to the DRL order. From the end-user point of view, a convenient way to express the solutions of such a system is to use a Rational Univariate Representation (RUR) [START_REF] Rouillier | Solving zero-dimensional systems through the Rational Univariate Representation[END_REF], [START_REF] Rouillier | Algorithmes pour l'étude des solutions réelles des systèmes polynomiaux[END_REF], i.e., a one-to-one mapping

V ( B ) -→ V ( p ) b = (b 0 , • • • , b n-1 ) -→ ξ p 0 (ξ) q(ξ) , • • • , p n-1 (ξ) q(ξ) ←- ξ, (22)
between the solutions of the system B defined by

V ( B ) := b ∈ K n | r(b) = 0, ∀ r ∈ B ,
and the roots of a univariate polynomial p. In order to achieve the one-to-one condition, the above representation is computed with respect to a separating element t, that is a linear combination t :=

n-1 i=0 u i b i of the b i 's, that takes different values when evaluated at different points of V ( B ).
Generically, i.e. for almost any value of the parameters, the system B admits t := b n-1 as a separating element. In that case, the above representation has the following form

                       p(b n-1 ) = 0, b n-1 = b n-1 , b n-2 = p n-2 (b n-1 ) q(b n-1 ) , • • • b 0 = p 0 (b n-1 ) q(b n-1 ) , (23)
where p, p 0 , . . . , p n-2 , q are in K[b n-1 ] and p and q coprime polynomials, i.e., gcd(p, q) = 1.

For instance, given explicit values for the parameters, using (23), we can compute certified numerical approximations of the real solutions by first isolating the real roots of p by means of intervals [START_REF] Rouillier | Efficient isolation of polynomial's real roots[END_REF], and then substituting these intervals into the expressions pi(bn-1) q(bn-1) for i = 0, . . . , n -2 to get the corresponding isolating intervals for the b i 's.

From the computation point of view, a RUR of B can be obtained using the algorithm given in [START_REF] Rouillier | Algorithmes pour l'étude des solutions réelles des systèmes polynomiaux[END_REF]. This algorithm requires the knowledge of a basis of A := R/ B (e.g., {b α } α∈ 0,1 n ; see the proof of Theorem 3) and a reduction algorithm which computes normal forms modulo the ideal B . In order to explicitly characterize the polynomials appearing in the RUR, we first introduce the multiplication K-endomorphism of A associated with r ∈ R defined by

m r : A -→ A a -→ r a,
where a denotes the residue class of a ∈ R in A (i.e., modulo B ). A representative of a is the normal form of a with respect to the Gröbner basis B. Using the basis {b α } α∈ 0,1 n , the K-endomorphism m r can be defined by means of a matrix M r ∈ K 2 n ×2 n . Then, p is the characteristic polynomial of M bn-1 . Moreover, if

p(b n-1 ) := p gcd p, dp dbn-1 = d i=0 v i b d-i n-1 ∈ K[b n-1 ]
denotes the square-free part of p and if we note

H j (b n-1 ) := j i=0 v i b j-i n-1 ∈ K[b n-1 ], j = 0, . . . , d -1,
then, for k = 0, . . . , n -2, we can define:

             q(b n-1 ) := d-1 i=0 Trace(M i bn-1 ) H d-i-1 (b n-1 ), p k (b n-1 ) := d-1 i=0 Trace(M b k M i bn-1 ) H d-i-1 (b n-1 ).
For more details, see [START_REF] Rouillier | Solving zero-dimensional systems through the Rational Univariate Representation[END_REF], [START_REF] Rouillier | Algorithmes pour l'étude des solutions réelles des systèmes polynomiaux[END_REF]. Note that the polynomials p, p 0 , . . . , p n-2 , q defining a RUR ( 22) are not unique. The above formulas give us a way to compute the RUR. But the polynomials appearing are usually not the shortest ones. For small orders (i.e., n < 5), a RUR can easily be obtained by direct computations yielding simple expressions. Finally, if b n-1 is a separating element and the ideal B is radical (i.e., B := {r ∈ R | ∃ l ∈ N : r l ∈ B } = B ), a similar representation can be computed by performing a change of order for the Gröbner basis (passing from the DRL order to the lexicographic order b 0 . . . b n-1 ) [START_REF] Rouillier | Algorithmes pour l'étude des solutions réelles des systèmes polynomiaux[END_REF]. However, the size of the outputs is then much more larger than the one given by a RUR [START_REF] Rouillier | Algorithmes pour l'étude des solutions réelles des systèmes polynomiaux[END_REF], a fact increasing the computational cost. For instance, on a regular machine (Intel Core i7 2.60 GHz processor and 16 GB of RAM), we can compute the RUR of B for n = 6 with all the parameters a i 's and c j 's in 20 minutes Maple 18 CPU time. By comparison, a similar computation based on a change of order for the Gröbner basis B did not finish after 3 days (due to a lack of memory).

Remark 3: If p is square-free (i.e., the discriminant Disc(p) ∈ K of p is not zero), then the ideal B is radical and b n-1 is a separating element. In a future work, the conditions of non degeneracy will be further studied.

Finally, combining (23) and ( 13), all the entries x i,j of X can be expressed in terms of b n-1 = x n,n -a n-1 (see [START_REF] Rouillier | Algorithmes pour l'étude des solutions réelles des systèmes polynomiaux[END_REF]).

Remark 4: Note that the above effective method yields the explicit resolution of an ARE for a general SISO system with parameters. It can be used to study many standard problems (optimal control and estimation problems such as H ∞ or H 2 optimal control, LQ control, model order reduction, Kalman filtering, . . . ) for SISO system with parameters.

Once an explicit expression of the solutions of R = 0 is known (see ( 23)), we then have to characterize the positive definite one, i.e. X of Theorem 1, among the 2 n solutions.

V. FROM ARE TO SPECTRAL FACTORIZATION

Another way to obtain a solution of R is by means of the invariant subspaces of the Hamiltonian matrix defined by:

H := A -B B T -C T C -A T ∈ K 2n×2n .
Using (4), the computation of the characteristic polynomial of H , denoted by f (λ) := det(H -λ I n ), gives:

f (λ) = (-1) n n i=0 n j=0
δ i,j λ i+j , δ i,j := (-1) i (a i a j + c i c j ) .

(24) We have:

f (λ) := λ 2n +(-1) n n-1 k=0     0≤i,j≤n i+j=2 k δ i,j λ 2k + 0≤i,j≤n i+j=2 k+1 δ i,j λ 2k+1     .
(25) Note that the coefficient of f of degree 2 k + 1 in λ vanishes for k = 0, . . . n -1. Moreover, some coefficients of λ 2k appear twice and after tedious computation, we obtain

f (λ) := (-1) n n k=0 (-1) k d 2k λ 2k , (26)
where d 2k is defined in (21). As a consequence, f is an even polynomial, i.e. f (-λ) = f (λ).

The even property of f leads to the spectral factorization problem which consists in finding a polynomial

g(λ) := n k=0 β k λ k such that f (λ) = (-1) n g(λ) g(-λ), (27)
where all the roots of g have negative real parts. We first concentrate on the condition (27). We have:

f (λ) = (-1) n n i=0 n j=0 (-1) i β i β j λ i+j . (28)
After similar tedious computation as the ones done to get (26) , using the change of variable i := k-m and j := k+m, we obtain

f (λ) = n k=0 (-1) n+k   β 2 k + 2 M (k) m=1 (-1) m β k-m β k+m   λ 2k , (29)
where M (k) is defined by (18). Combining (26) and (29), we finally obtain

   β 2 0 = d 0 , β k 2 + 2 M (k) m=1 (-1) m β k-m β k+m = d 2k , for 1 ≤ k ≤ n -1, which coincides with (20) with β k = b k for k = 0, . . . , n -1.
Let us now characterize b n-1 = β n-1 so that all the roots of g have negative real parts, or equivalently characterize b n-1 = β n-1 which defines the unique positive definite solution X of ( 6). In the sequel, the parameters are fixed to explicit values. Since f is even, if λ is a solution of f , then so is -λ. Let {α i } i=1,...,n be the roots of f with negative real parts (see ( 26)). We want to determine b n-1 = β n-1 , and thus all the β i = b i 's by ( 23) such that all the α i 's are roots of g. Given g satisfying ( 27), then α i is either a root of g(λ) = 0 or of g(-λ) = 0. Since β n = 1, -β n-1 is a sum of the roots of g, we get β n-1 = -n i=1 i α i where i ∈ {-1, 1}. If {ν k } k=1,...,2 n is the set of complex solutions of B, where ν k := (ν k,0 , . . . , ν k,n-1 ), and ν k,n-1 is the n th coordinate of ν k , using

β n-1 = b n-1 , then we get ν k,n-1 = - n i=1 k,i α i for = 1, . . . , 2 n , where k,i ∈ {-1, 1}. Moreover, we have σ := - n i=1 α i = max k∈{1,...,2 n } ν k,n-1 , (30)
For more details, see [START_REF] Kanno | Solution of algebraic Riccati equations using the sum of roots[END_REF]. We obtain the following theorem. Theorem 4 ( [START_REF] Kanno | Solution of algebraic Riccati equations using the sum of roots[END_REF]): The positive definite solution X of R = 0 corresponds to the real solution

ν k of B such that b n-1 (ν k ) is maximal.
Proposition 2: We have b 0 = a 2 0 + c 2 0 . Proof: Using ( 19) and ( 21), we have b 2 0 = a 2 0 + c 2 0 . We also note that we have

b 0 = β 0 = (-1) n n i=1 λ i ,
where λ i 's are roots of g with negative real parts. Let us distinguish the following two cases:

• If n is even then b 0 = (-1) n n/2 i=1 λ i λ i = n/2 i=1 |λ i | 2 > 0. • If n is odd then we have b 0 = (-1) n λ 0 (n-1)/2 i=1 λ i λ i = (-1) λ 0 (n-1)/2 i=1 |λ i | 2 ,
where λ 0 < 0 so that b 0 > 0 and concludes the proof.

VI. CERTIFIED COMPUTATION OF γ opt

As stated in Theorem 1, an optimal solution of the RCP is given by [START_REF] Kanno | Solution of algebraic Riccati equations using the sum of roots[END_REF]. Using the RUR (23) of B and ( 13), we can express the entries x i,j (resp., y i,j ) of X (resp., Y ) in terms of b n-1 = x n,n -a n-1 , and then we compute the characteristic polynomial of QX:

h(λ, b n-1 ) := det(QX -λ I n ).
The numerator of h, denoted by h , belongs to K [λ, b n-1 ]. Hence, the problem of computing γ opt is equivalent to that of the maximal real solution of h (λ, σ) = 0, where σ is defined by (30) (i.e., the maximal real root of p). We note again that σ defines the solutions X and Y .

When explicit values of the parameters are given, this problem can be tackled using purely numerical methods.

For instance, one approach consists in first isolating the root σ of p by means of an interval, substituting this interval into h (λ, b n-1 ) and isolating again in order to find λ (see Corollary 1). But such a numerical method fails to provide a certified output (see [START_REF] Collins | Polynomial real roots isolation using Descartes' rule of signs[END_REF], [START_REF] Rouillier | Efficient isolation of polynomial's real roots[END_REF]), especially when h (λ, σ) is not square-free in λ. In order to achieve a certification, we can instead consider the following triangular system

p(b n-1 ) = 0, h (λ, b n-1 ) = 0, (31)
and compute a certified numerical approximation of the solution (σ, λ ) [START_REF] Collins | Polynomial real roots isolation using Descartes' rule of signs[END_REF], [START_REF] Rouillier | Efficient isolation of polynomial's real roots[END_REF].

There exist several exact methods for solving (31) and thus for deducing the desired solution (σ, λ ). We can mention the algorithm of [START_REF] Cheng | Complete numerical isolation of real roots in zero-dimensional triangular systems[END_REF] which uses some separation bounds in order to isolate the roots of the polynomial h (λ, b n-1 ) over a root of p, or the algorithm given in [START_REF] Boulier | Real root isolation of regular chains[END_REF] which requires the computation of a decomposition of (31) into a set of square-free triangular systems so that for any root α of p the polynomial h (λ, α) is square-free.

We now use again the idea developed in Section IV and compute a RUR of the solutions of (31). For a well chosen separating element τ := u 1 b n-1 + u 2 λ, these solutions can be expressed by:

             v(τ ) = 0, b n-1 = v 1 (τ ) w(τ ) , λ = v 2 (τ ) w(τ ) . (32)
Finally, in order to obtain the solution (σ, λ ), it is sufficient to isolate the real roots of v, evaluate the resulting intervals in the expressions v1(τ ) w(τ ) and v2(τ ) w(τ ) , which yields boxes isolating the solutions, and finally choose the one that satisfies the requirement (i.e., the maximal λ obtained from the maximal b n-1 ).

Getting λ , we deduce γ opt := 1 + λ 2 , which allows us to compute an optimal H ∞ controller by Theorem 1.

Finally, we recapitulate the above approach for solving the RCP by outlining an algorithm.

Algorithm 1:

1) Using ( 20) and ( 21), compute B.

2) Compute a RUR (23) of the ideal B .

3) Using Theorem 2, compute X in terms of b n-1 . 4) Compute Q defined by (9). 5) Compute Q X and the numerator of its characteristic polynomial h . 6) Compute a RUR of (31) to get (32). 7) Isolate the solutions of (32) by means of 2D-boxes and select the desired box corresponding to (σ, λ ). 8) Compute γ opt := 1 + λ 2 . 9) Compute Y := Q X Q. 10) For γ > γ opt , compute the controller K γ , defined by [START_REF] Glover | Robust stabilization of normalized coprime factor plant descriptions with H∞-bounded uncertainty[END_REF], which satisfies the RCP.

VII. EXAMPLE: A TWO-MASS-SPRING SYSTEM

We illustrate the above approach with the two-massspring problem (Figure 2) considered in [START_REF] Wie | A benchmark problem for robust control design[END_REF], [START_REF] Vinnicombe | Uncertainty and Feedback -H∞ Loop-shaping and the ν-gap metric[END_REF], which mathematical model is similar to the one of a 2-dimensional gyro-stabilized sight 1 . Fig. 2. Two-mass-spring system Two masses m 1 and m 2 are linked by a spring of stiffness k. With the notations of Figure 1, we study the displacement of m 2 , denoted by y 1 , while m 1 is excited by a force e 1 . We consider the transfer function of the physical plant

P := c s 2 (s 2 + a 2 ) , c := k m 1 m 2 , a 2 := m 1 + m 2 m 1 m 2 k,
between the input e 1 and the output y 1 . As in [14, §2.6, §4], we consider a weight w > 0 and define the fictive plant by

G := w P = c 0 s 2 (s 2 + a 2 )
, c 0 := w c > 0, a 2 > 0, which will be used in Algorithm 1. As seen above, we can design the loop-shaping controller K γ stabilizing G.

From this controller K γ , we then get the robust stabilizing controller C γ := w K γ of the physical plant P .

For the system G, the polynomials B k 's defined in Theorem 3, are given by

         B 0 := b 2 0 -d 0 = 0, B 1 := b 2 1 -2 b 0 b 2 -d 2 = 0, B 2 := b 2 2 -2 b 1 b 3 + 2 b 0 -d 4 = 0, B 3 := b 2 3 -2 b 2 -d 6 = 0, (33)
where the constants d 2k 's are defined by:

d 0 := c 2 0 , d 2 := 0, d 4 := a 2 2 , d 6 := -2 a 2 .
Using Proposition 2, we get b 0 = c 0 . Moreover, if b 3 = 0, then the last equation of (33) yields b 2 = a 2 , and thus the last but one gives b 0 = c 0 = 0, i.e. G = 0. Since in the sequel we suppose that G = 0, and thus, b 3 = 0, a RUR of the system {B 1 = 0, B 2 = 0, B 3 = 0} for the separating element b 3 can easily be obtained by direct computations

         b 1 = b 4 3 + 4 a 2 b 2 3 + 8 c 0 8 b 3 , b 2 = 1 2 b 2 3 + a 2 , b 3 = b 3 , (34)
where b 3 satisfies the following polynomial equation:

p(b 3 ) := b 8 3 +8 a 2 b 6 3 +16 (a 2 2 -3 c 0 ) b 4 3 -64 a 2 c 0 b 2 3 +64 c 2 0 = 0. (35)
Noting τ = b 2 3 / √ 8 c 0 , we obtain:

c 0 τ 4 + √ 8 c 0 a 2 τ 3 + 2 (a 2 2 -3 c 0 ) τ 2 - √ 8 c 0 a 2 τ + c 0 = 0.
(36) The latter polynomial is anti-palindromic, so we can apply the change of variable ς := τ -τ -1 in order to get:

p s (ς) := c 0 ς 2 + 2 √ 2 c 0 a 2 ς + 2 a 2 2 -4 c 0 = 0. (37)
The discriminant of (37) is 16 c 2 0 = 0. Thus, p s has two distinct real roots, which are defined by:

ς + := - 2 c 0 a 2 + 2, ς -:= - 2 c 0 a 2 -2.
Since a 2 > 0, we have

ς -< 0, ς + > ς -, |ς -| > |ς + |. Then, τ verifies τ 2 -τ ς -1 = 0 which discriminant is ς 2 + 4 > 0.
Its solutions are given by:

       τ + (ς ± ) := 1 2 ς ± + ς 2 ± + 4 , τ -(ς ± ) := 1 2 ς ± -ς 2 ± + 4 , ς ± ∈ {ς + , ς -}.
Then all the roots of (36) are real. It is clear that τ + (ς + ) is the maximal real root of (36), which leads to the following expression of σ, maximal real root of p:

σ := √ 2 √ 2 c 0 -a 2 + √ 2 c 0 -a 2 2 + 2 c 0 . ( 38
)
Using Theorem 2 (see ( 13)), the matrix X is defined by

X :=        c 0 b 1 c 0 b 2 c 0 b 3 c 0 c 0 b 2 b 1 b 2 -c 0 b 3 b 1 b 3 -c 0 b 1 c 0 b 3 b 1 b 3 -c 0 b 2 b 3 -b 1 b 2 -a 2 c 0 b 1 b 2 -a 2 b 3       
, where b 1 and b 2 are defined by (34). By (9), we have:

Q :=        0 0 0 c -1 0 0 0 c -1 0 0 0 c -1 0 0 -a 2 c -1 0 c -1 0 0 -a 2 c -1 0 0       
.

The numerator h of the characteristic polynomial of QX divided by c 2 0 , i.e., h s := h c -2 0 , is defined by: h

s (λ, b 3 ) = 4 c 0 λ 4 -α(b 3 ) λ 3 + β(b 3 ) λ 2 + α(b 3 ) λ + 4 c 0 , ( 39
) α := b 4 3 + 8 c 0 , β := -2 (b 4
3 + 4 c 0 ). Since h s is anti-palindromic, using t := λ -λ -1 (λ = 0 is not a root of h s since its constant term 4 c 0 is non-zero), h s can then be rewritten as:

h s (t, b 3 ) = 4 c 0 t 2 -α(b 3 ) t -2 b 4 3 . ( 40
)
The real roots of (40) are then given by:

       t + (b 3 ) := 1 8 c 0 α(b 3 ) + α(b 3 ) 2 + 32 c 0 b 4 3 , t -(b 3 ) := 1 8 c 0 α(b 3 ) -α(b 3 ) 2 + 32 c 0 b 4 3 .
The roots of h are then of the form:

       λ + (t ± ) := 1 2 t ± + t 2 ± + 4 , λ -(t ± ) := 1 2 t ± -t 2 ± + 4 , t ± ∈ {t + , t -}.
Let us determine the root of h of greatest absolute value. Since t + t -= -2 b 4 3 4 c0 < 0 and t + > t -, t -is necessarily negative while t + is positive. Furthermore, since t + + t -= α(b3)

4 c0 > 0, then we get t + = |t + | > -t -, i.e., |t + | > |t -|.
The root of h of greatest absolute value is given by:

λ := λ + (t + ) = 1 2 t + + t 2 + + 4 . (41)
Hence, using (38), the minimal γ is given by:

γ opt (σ) := 1 +   t + (σ) 2 + 1 + t + (σ) 2 2   2 . ( 42
)
Remark 5: (42) is also valid for a 2 < 0.

Using the explicit formula [START_REF] Glover | Robust stabilization of normalized coprime factor plant descriptions with H∞-bounded uncertainty[END_REF] given in Theorem 1, we can compute a robust stabilizing controller K γ of G and, since K γ G = C γ P , we deduce a stabilizing controller C γ of P satisfying S ∞ < γ and G K γ S ∞ < γ [14, Cor. 5.1]. Figure 3 shows the plot of the function γ opt in terms of the values of the parameters a 2 and c 0 . We note that we have: The weight w can be seen as a scalar tuning parameter for γ opt . As a consequence, for γ > γ opt , we can obtain an auto-tuned controller C γ (c, a 2 ) = K γ (c, a 2 )/w. By on-line identification of the modal frequency a 2 and of the gain c, the controller C γ (c, a 2 ) can be computed in real-time by a simple embedded calculator (which does not have to contain optimization routines).

γ inf = inf a2>0, c0>0 γ opt = 4 + 2 √ 2 ≈ 2.6131.
According to [14, Th.2.10], we have the following gain and phase margins:

   ∆ G (P, C γ ) = ∆ G (G, K γ ) ≥ 1 + γ -1 1 -γ -1 , ∆ Φ (P, C γ ) = ∆ Φ (G, K γ ) ≥ 2 arcsin γ -1 .
For γ opt = 3, with γ opt > γ inf , we have

∆ G (P, C γopt ) = ∆ G (G, K γopt ) ≥ 6 dB, ∆ Φ (P, C γopt ) = ∆ Φ (G, K γopt ) ≥ 39 • ,
and to achieve γ opt = 3, w can be chosen as follows:

w = 952 6561 √ 2 -8 √ 7 223074 √ 2 -129472 72048449 -49968576 √ 2 a 2 2 c 0 .
Remark 6: Since we usually cannot solve by radicals univariate polynomials of degree larger than or equal to 4, we have to use the symbolic-numeric method described in Algorithm 1. Let us show how to compute the second RUR at Step 6 of Algorithm 1 for the above example. We can first compute a RUR of the polynomial system defined by p and h s given by ( 35) and ( 40 (43) After some eliminations, we obtain T = δ(t) θ -1 (t), where

     δ(t) = -c 0 t 4 +
Substituting t = λ -λ -1 into (44) and (45), we obtain a RUR for the polynomial system defined by ( 35) and (39).

For explicit values of a 2 and c 0 , root isolation techniques can then be used to compute the H ∞ controller. For instance, for a 2 = 10 and c 0 = 1, the above computations can be directly obtained by applying Maple command RootFinding[Isolate] to (43) with options method = "RS", output=interval. In our example, we obtain:

VIII. CONCLUSION

In this paper, we have initiated a new symbolic-numeric method for solving the H ∞ loop-shaping design problem for low order SISO systems with parameters. To do that, we first showed that solving an ARE of size n × n is equivalent to solving a certain zero-dimensional polynomial system of n equations in n unknowns, and that a RUR can be used to parametrize all its real solutions. Then, using a result of spectral factorization [START_REF] Kanno | Solution of algebraic Riccati equations using the sum of roots[END_REF], we showed how to solve an H 2 control problem. The computation of the controller satisfying the robust 4-block problem requires a second RUR computation and symbolic/numeric root isolation techniques.

When the univariate polynomials of the two RUR can be reduced to polynomials of degrees less than or equal to 4, then an explicit formula of the H ∞ controller can be obtained. Then, the dependence of this controller to the system parameters can then be analytically studied, which is useful, e.g., in the applications to gyro-stabilized sights. For higher order systems, a certified numerical root isolation technique has to be used to compute the H ∞ controller.

Further works include the search for tractable RURs in order to deal with highest order systems. Furthermore, our approach will be applied to small order models of gyrostabilized sights to develop an adaptive control scheme.

  and the references therein), B is a Gröbner basis with respect to the DRL order b n-1 . . . b 0 . Moreover, this Gröbner basis is reduced according to [4, ch. 2, §7, Definition 4]. Since the leading monomial of each polynomial in the Gröbner basis B is of the form b 2 k , by [4, ch. 5, §6, Theorem 4], A is a finite-dimensional Kvector space. Using the notation b α := n-1 i=0 b αi i

Fig. 3 .

 3 Fig. 3. γopt = Γ(a 2 , c 0 )

  ) with T = b 2 3 , i.e.:T 4 + 8 a 2 T 3 + 16 (a 2 2 -3 c 0 ) T 2 -64 a 2 c 0 T + 64 c 2 0 = 0, 4 c 0 t 2 -T 2 + 8 c 0 t -2 T 2 = 0.

σ

  = [0.481024117223, 0.481024117224], λ = [2.436943741649, 2.436943741650].

  -4 a 2 2 + 16 c 0 t 3 -8 c 0 t 2 +16 a 2 2 -4 c 0 t -16 c 0 , θ(t) = 2 a 2 t 3 -2 t 2 -12 t -8 ,

	(44)
	and t satisfies the following polynomial equation:
	c 2 0 t 8 -8 c 0 (a 2 2 + 4 c 0 ) t 7 + 16 (a 4 2 + 17 c 2 0 ) t 6
	-32 c 0 (a 2 2 + 4 c 0 ) t 5 -32 (4 a 4 2 + 61 c 2 0 ) t 4
	+128 c 0 (a 2 2 + 4 c 0 ) t 3 + 256 (a 4 2 + 17 c 2 0 ) t 2
	+512 c 0 (a 2 2 + 4 c 0 ) t + 256 c 2 0 = 0.

A gyro-stabilized sight is a system of cameras controlled in rotation by motors and gyroscopes.