Y Bouzidi
email: yacine.bouzidi@inria.fr.

A Quadrat

F Rouillier

Certified Non-conservative Tests for the Structural Stability of Discrete Multidimensional Systems

Keywords: Multidimensional systems, structural stability, stability analysis, computer algebra

In this paper, we present new computer algebra based methods for testing the structural stability of n-D discrete linear systems (with n ≥ 2). More precisely, we show that the standard characterization of the structural stability of a multivariate rational transfer function (namely, the denominator of the transfer function does not have solutions in the unit polydisc of C n ) is equivalent to the fact that a certain system of polynomials does not have real solutions. We then use state-of-the-art computer algebra algorithms to check this last condition, and thus the structural stability of multidimensional systems.

Introduction

Multidimensional systems is a class of systems for which the information propagates in more than just one dimension as for the classical dynamical systems (this dimension being the continuous/discrete time). The latter class of systems is usually referred as 1-D systems whereas multidimensional systems are also called n-D systems, where n denotes the number of dimensions in which the information propagates. Within the frequency domain approach, such a system is defined by means of a transfer function G of the following form

G(z 1 , . . . , z n ) = N (z 1 , . . . , z n ) D(z 1 , . . . , z n ) , (1)
where D and N are two polynomials in the complex variables z 1 , . . . , z n with real coefficients, i.e., D, N ∈ R[z 1 , . . . , z n ], which are supposed to be factor prime i.e., gcd(D, N ) = 1, where gcd stands for the greatest common divisor of D and N . As for 1-D systems, a fundamental issue in the multidimensional systems theory is stability analysis. In this paper, we are interested in testing the structural stability of discrete linear multidimensional systems defined by multivariate rational transfer function. A system such as (1) is said to be structurally stable if the denominator D of G is devoid of zeros in the complex unit polydisc U n defined by:

U n := n k=1 {z k ∈ C | |z k | ≤ 1}.
In other words, if V C (D) = {z = (z 1 , . . . , z n ) ∈ C n | D(z) = 0} denotes the hypersurface formed by the complex zeros of D, then (1) is structurally stable if the following condition holds:

V C (D) ∩ U n = ∅. (2)
In the context of discrete linear multidimensional systems, structural stability implies bounded-input bounded-output (BIBO) stability and, under certain conditions, other types of stability such as asymptotic or exponential stability. For further details on the connections between the different concepts of stability of multidimensional systems, see Bachelier et al. (2016a) and the references therein.

The simplicity of (2) significantly contrasts with the difficulty to develop effective algorithms and efficient implementations for testing it. One important first step toward this objective was the formulation of new conditions that are equivalent to the above condition (see ( 2)) but easier to handle. The following theorems, due to [START_REF] Strintzis | Tests of stability of multidimensional filters[END_REF] and DeCarlo et al. [START_REF] Decarlo | Multivariable Nyquist theory[END_REF], are two good representatives of these reformulations.

Theorem 1 [START_REF] Strintzis | Tests of stability of multidimensional filters[END_REF]) Condition ( 2) is equivalent to: Theorem 2 [START_REF] Decarlo | Multivariable Nyquist theory[END_REF]) Condition ( 2) is equivalent to:

                 D(
                 D(z 1 , 1, . . . , 1) = 0, |z 1 | ≤ 1, D(1, z 2 , 1, . . . , 1) = 0, |z 2 | ≤ 1, . . . . . . D(1, . . . , 1, z n ) = 0, |z n | ≤ 1, D(z 1 , . . . , z n ) = 0, |z 1 | = . . . = |z n | = 1.
(1994); [START_REF] Hu | On two-dimensional filter stability test[END_REF]; [START_REF] Bistritz | Testing stability of 2-D discrete systems by a set of real 1-D stability tests[END_REF]; [START_REF] Xu | Further improvements on Bose's 2-D stability test[END_REF]; [START_REF] Fu | Author's reply to a counterexample to generalized eigenvalue-based stability tests for 2-d linear systems: Necessary and sufficient conditions[END_REF] and the references therein. Common to all these tests is the fact that they proceed recursively on the number of variables, reducing the computations with polynomials in two variables to computations with a set of univariate polynomials with algebraic coefficients by means of symbolic computation tools such as resultant and subresultant polynomials (see, e.g., [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]). Such a recursive approach, which shows its relevance for 2-D systems, becomes rather involved when it comes to n-D systems with n > 2, mainly due to the exponential increase of the degree of the intermediate polynomials. This fact prevents these 2-D tests from being efficiently generalized to n-D systems.

For n-D systems with n > 2, very few implementations for the stability analysis exist. Among the recent work on this problem, one can mention the work of [START_REF] Serban | Multidimensional systems: Bibo stability test based on functional schur coefficients[END_REF] where, using an extension of the 1-D Schur-Cohn criterion, a new stability condition is proposed as an alternative to the conditions of Theorems 1 and 2. As a result, the stability is expressed as a positivity condition of n -1 polynomials on the unit polycircle:

T n-1 := n-1 k=1 {z k ∈ C | |z k | = 1} .
Unfortunately, such a condition becomes considerably hard to effectively test when the involved systems are not of low degree in few variables. To achieve practical efficiency, [START_REF] Dumitrescu | Stability test of multidimensional discrete-time systems via sum-ofsquares decomposition[END_REF][START_REF] Dumitrescu | LMI stability tests for the Fornasini-Marchesini model[END_REF] proposes a sum of squares approach to test the last Decarlo's condition (Theorem 2). The proposed method is however conservative, i.e., it provides only a sufficient stability condition. Finally, LMI test also exist for n = 2 (Bachelier et al. (2016b,a)). To sum up, the existing stability tests for n-D systems with n > 2 are either nonconservative but inefficient or efficient (polynomial time) but conservative.

The contribution of this paper is threefold. First, a new algebraic approach for testing the stability of n-D systems is presented. Our approach, which starts with the stability conditions given by [START_REF] Decarlo | Multivariable Nyquist theory[END_REF] (see Theorem 2), transforms the problem of testing these conditions into that of deciding the existence of real zero in some algebraic or semi-algebraic sets. Hence, state-of-the-art real algebraic geometry techniques can then be used for this purpose. Unlike the existing counterparts, this new approach is not conservative. Moreover, our approach shows good practical performances for relatively small dimensions n.

Secondly, we address the specific case of 2-D systems with the main objective of achieving practical efficiency. Following the same approach as for n-D systems but taking advantage from the recent developments in solving bivariate algebraic systems of equations (see Bouzidi (2014b)), we propose a stability test based on the existence of real solutions of bivariate algebraic systems which is efficient in practice. Namely, this test makes use of the software RS3 [START_REF] Rouillier | RS3: Real roots isolation for bivariate algebraic systems of equations with rational coefficients with a finite number of complex roots[END_REF]) which provides very efficient tools for the symbolic solving of bivariate systems of equations.

Finally, the above 2-D stability test is extended in order to handle system parameters. More precisely, using the concept of discriminant variety developed in the computer algebra community [START_REF] Lazard | Solving parametric polynomial systems[END_REF]), we provide a new method which, given a 2-D system depending on an arbitrary set of real parameters, decomposes the parameter space into regions inside which the stability is guaranteed.

The plan of the paper is the following. In Section 2.1, we first reformulate the last condition of Theorem 1 as the emptiness of a certain semi-algebraic set. We then present state-of-part computer algebra techniques (namely, Cylindrical Algebraic Decomposition, critical point methods and Rational Univariate Representation) which can be used to effectively study this last problem. In Section 2.2, we present a new approach for testing the last condition of Theorem 1 based on the Möbius transform and a critical point method. Algorithms are presented. They are then illustrated with explicit examples in Section 2.3 and their implementations in Maple are discussed and timings are shown. A new stability test for 2-D systems is presented in Section 3 based on a recent approach, developed by two of the authors, for the efficient computation of real solutions of systems of two polynomial equations in two variables. This approach is based on the efficient computation of Rational Univariate Representations using resultants, subresultant sequences and triangular polynomial systems. Finally, Section 4 shows how to use the mathematical concept of a characteristic variety, developed by one of the author, to study the stability of 2-D systems with parameters. We show how the parameter space can be explicitly decomposed into cells so that the 2-D system is either stable or unstable for all values of the parameters belonging to each cell.

Stability of n-D discrete linear systems

In this section, Subsection 2.1 overviews in broad lines computer algebra methods for computing the real zeros of semi-algebraic sets (namely, unions of sets defined by a finite number of polynomial equations and inequalities) and recall the basic ideas behind these methods. Then, Subsection 2.2 shows how we can obtain new stability conditions that can be tested efficiently using the previously introduced computer algebra methods. Finally, Subsection 2.3 illustrates our stability test on non-trivial examples and show its practical efficiency thanks to experimental tests.

Computer algebra methods

Recall that the transfer function G defined by ( 1) is said to be structurally unstable if the set

E := V C (D) ∩ U n = {(z 1 , . . . , z n ) ∈ C n | D(z 1 , . . . , z n ) = 0, |z 1 | ≤ 1, . . . , |z n | ≤ 1}
is not empty. The set E is a semi-algebraic subset of R 2n . Indeed, if we note z k := x k + i y k , where x k (resp., y k ) is the real part (resp., the imaginary part) of z k and i is the imaginary unit, then the polynomial D can be rewritten as follows

D(z 1 , . . . , z n ) = R(x 1 , . . . , x n , y 1 , . . . , y n ) + i I(x 1 , . . . , x n , y 1 , . . . , y n ),
where R, I ∈ R[x 1 , . . . , x n , y 1 , . . . , y n ], and the inequalities |z k | ≤ 1 are equivalent to x 2 k + y 2 k ≤ 1 for k = 1, . . . , n, which shows that E is equivalently defined by the following semi-algebraic set:

E ∼ = {(x 1 , . . . , x n , y 1 , . . . , y n ) ∈ R 2n | R(x 1 , . . . , x n , y 1 , . . . , y n ) = 0, I(x 1 , . . . , x n , y 1 , . . . , y n ) = 0, x 2 k + y 2 k ≤ 1, k = 1, . . . , n}. (3)
Testing ( 2) is thus equivalent to testing that the above semi-algebraic set is empty. This test can be performed using classical computer algebra methods for computing the real zeros of semi-algebraic systems which will be briefly overviewed in the next section.

Cylindrical Algebraic Decomposition & Critical point methods

To study the real zeros of semi-algebraic sets, two classes of symbolic algorithms are available: the algorithms based on Cylindrical Algebraic Decompositions (CAD) and those based on the study of the critical points of well-chosen functions (see, e.g., [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]).

The Cylindrical Algebraic Decomposition (CAD): Introduced originally by [START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition[END_REF]) in the seventies, the CAD has become a standard tool for the study of real zeros of semi-algebraic sets. CAD refers to both an object and an algorithm for computing this object. In short, a CAD associated to a finite set of polynomial F = {P 1 , . . . , P s } ⊂ R[x 1 , . . . , x n ] is a partition of R n into connected semi-algebraic sets, called cells, on which each polynomial P i has a constant sign (i.e., either +, -or 0). For instance, the CAD of a set of univariate polynomials in R[x] is an union of points and open intervals that form a partition of R. Such a partition is called F -invariant. Let Π k : R n -→ R n-k denote the projection onto the first n -k components of R n . The CAD is called cylindrical since for every two cells c 1 and c 2 , we either have

Π k (c 1 ) = Π k (c 2 ) or Π k (c 1 ) ∩ Π k (c 2 ) = ∅.
This implies that the images of the cells by Π k define a cylindrical decomposition of R n-k .

Example 1 A CAD associated to

P = x 2 1 + x 2 2 -1 ∈ R[x 1 , x 2
] is a partition of R 2 into the following algebraic sets (cells) in each of which the polynomial P has a constant sign:

-C 1 = {(x 1 , x 2 ) ∈ R 2 | x 1 < -1}, -C 2 = {(x 1 , x 2 ) ∈ R 2 | x 1 = -1, x 2 < 0}, -C 3 = {(x 1 , x 2 ) ∈ R 2 | x 1 = -1, x 2 = 0}, -C 4 = {(x 1 , x 2 ) ∈ R 2 | x 1 = -1, x 2 > 0}, -C 5 = {(x 1 , x 2 ) ∈ R 2 | -1 < x 1 < 1, x 2 1 -x 2 2 -1 > 0, x 2 < 0}, -C 6 = {(x 1 , x 2 ) ∈ R 2 | -1 < x 1 < 1, x 2 1 -x 2 2 -1 = 0, x 2 < 0}, -C 7 = {(x 1 , x 2 ) ∈ R 2 | -1 < x 1 < 1, x 2 1 -x 2 2 -1 < 0}, -C 8 = {(x 1 , x 2 ) ∈ R 2 | -1 < x 1 < 1, x 2 1 -x 2 2 -1 = 0, x 2 > 0}, -C 9 = {(x 1 , x 2 ) ∈ R 2 | -1 < x 1 < 1, x 2 1 -x 2 2 -1 > 0, x 2 > 0}, -C 10 = {(x 1 , x 2 ) ∈ R 2 | x 1 = 1, x 2 < 0}, -C 11 = {(x 1 , x 2 ) ∈ R 2 | x 1 = 1, x 2 = 0}, -C 12 = {(x 1 , x 2 ) ∈ R 2 | x 1 = 1, x 2 > 0}, -C 13 = {(x 1 , x 2 ) ∈ R 2 | x 1 > 1}.
The CAD algorithm mainly consists of two distinct phases: a projection phase and a lifting phase. During the projection phase, one proceeds recursively on the number of variables. Starting from the initial set of polynomials F 1 = {P 1 , . . . , P s }, a first set of polynomials F 2 ⊂ R[x 1 , . . . , x n-1 ] is computed which has the property that a partition of R n-1 that is F 2 -invariant, naturally lifts to a partition of R n that is F 1 -invariant. Then, from F 2 , another set F 3 ⊂ R[x 1 , . . . , x n-2 ] with the same property is computed and so on, until obtaining a set of univariate polynomials F n . At each step, the projection essentially consists in computing resultants (and subresultants) for all possible pairs of polynomials as well as their discriminants (see, e.g., [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]).

Example 2 The projection phase for Example 1 consists in only one projection and yields the discriminant -4 x 2 1 + 4 of P = x 2 1 + x 2 2 -1 with respect to x 2 .

Starting from the set of univariate polynomials, the lifting step, which also proceeds recursively, consists in isolating the real roots of univariate polynomials with real algebraic numbers as coefficients, which can be viewed as solving a socalled triangular zero-dimensional system, namely, a system with a finite number of complex solutions having a triangular form.

Example 3 Continuing Example 2, the real roots of -4 x 2 1 + 4 are first isolated, which yields a partition of R that consists of the cells ] -∞, -1[, {-1}, ] -1, 1[, {1} and ]1, +∞[. Then, the real roots of the polynomial x 2 1 + x 2 2 -1 are isolated for each x 1 in these cells, which yields the partition of R 2 given in Example 1.

Although the computation of the CAD answers our problem i.e,. deciding the existence of real zeros, the partition of R n given by the CAD gives more information than required. Moreover, its computation requires a number of arithmetic operations which is doubly exponential in the number of variables of the polynomial ring, i.e., n, due to, at least, the iterative computation of the resultants (and subresultants). Despite of this, it is worth mentioning that, to some extent, most of the existing algorithms for testing the stability of multidimensional systems can be viewed as particular variants of the CAD.

Critical point methods: Instead of a complete partition of the real space, critical points based methods basically compute at least one point in each real connected component of a given semi-algebraic set, which is sufficient to answer our question. Roughly speaking, the basic principle of these methods consists in computing the critical points of a well-chosen function restricted on an algebraic set. Under certain conditions, the set of critical points is defined by a zero-dimensional polynomial system (i.e., which admits a finite number complex solutions) and it meets each real connected component of the original semi-algebraic set.

Example 4 We consider again the polynomial given in Example 1, i.e., P = x 2 1 + x 2 2 -1. Our goal is to compute at least one point in the single real connected component of the algebraic set

V R (P ) = {α = (α 1 , α 2 ) ∈ R 2 | P (α) = 0}. Let Π : R 2 -→ R denote the projection onto the first component x 1 of (x 1 , x 2 ) ∈ R 2 .
The critical points of Π restricted to V R (P ) are the points of V R (P ) on which the derivative ∂P ∂x 2 vanishes. They are given as the real zeros of the following system

   P = x 2 1 + x 2 2 -1 = 0, ∂P ∂x 2 = 2 x 2 = 0,
which yields (1, 0) and (-1, 0). These points belong to the real connected component of P = 0.

Let us consider the real hypersurface

V R (P ) := {α = (α 1 , . . . , α n ) ∈ R n | P (α) = 0} defined by P ∈ Q[x 1 , . . . , x n ].
Let us suppose that V R (P ) is smooth, i.e., it does not contain singular points, namely, points where the rank of the Jacobian matrix is not maximal, and V R (P ) is compact (i.e., closed and bounded) for the so-called Zariski topology (see, e.g., [START_REF] Cox | Varieties, and Algorithms[END_REF]). The set of critical points of the projection with respect to some variable x i restricted to V R (P ) defined by

C Π (V R (P )) := α ∈ V R (P ) | ∂P ∂x k (α) = 0, ∀ k ∈ {1, . . . , i -1, i + 1, . . . , n} , (4
) is finite and meets each real connected component of V R (P ) (see [START_REF] Bank | Polar varieties and efficient real elimination[END_REF]; Safey El Din and Schost ( 2003)).

There are several ways to circumvent the hypotheses (i.e., compactness and smoothness). For instance, we can use a distance function to a fixed point instead of the projection function to get rid of the compactness assumption [START_REF] Seidenberg | A new decision method for elementary algebra[END_REF]), deform the variety to get a compact and smooth one [START_REF] Rouillier | Finding at least one point in each connected component of a real algebraic set defined by a single equation[END_REF]), introduce more general notions of critical points [START_REF] Bank | Generalized polar varieties: Geometry and algorithms[END_REF]), or separately study the subsets of singular points of the variety [START_REF] Aubry | Real solving for positive dimensional systems[END_REF]).

For systems of polynomial equations and, more generally systems of polynomial inequalities, several strategies have been proposed by different authors (see, e.g., [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]; [START_REF] Bank | Polar varieties and efficient real elimination[END_REF]). Some are based on the use of sums of squares to reduce the problem of studying an algebraic set to the problem of studying an hypersurface [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]), on the use of infinitesimal deformations by adding some variables to avoid singularities and to deal with inequalities [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]), or on the introduction of a special kind of critical points (generalized critical values) to circumvent the compactness hypothesis. But the basic ideas stay the same.

As already said, critical points methods compute less information than the CAD but they are sufficient in our case since we just have to decide if a semialgebraic set is empty. Moreover, a key advantage of these methods is that they transform the problem into solving a zero-dimensional polynomial system and this transformation is performed within a number of arithmetic operations that is single exponential in the number n of variables.

Remark 1 In practice, for systems depending on strictly less than 3 variables, the use of CAD is usually preferred since, it provides more information than the critical point methods for a negligible additional cost.

Symbolic resolution of univariate polynomials and zero-dimensional systems

The methods described in Section 2.1.1 are based on the resolution of univariate polynomials and, more generally, on zero-dimensional polynomial systems. For stability analysis of multidimensional systems, we mainly have to decide whether or not a polynomial system admits real solutions. For polynomial systems with a finite number of solutions, we can use an additional processing that turns this last problem into computing a univariate parameterization of all the solutions.

Given a zero-dimensional polynomial system and I ⊂ Q[x 1 , . . . , x n ] the ideal generated by the corresponding polynomials, a Rational Univariate Representation (RUR) of V C (I) is defined by a separating linear form t := a 1 x 1 + . . . + a n x n and univariate polynomials f, f 1 , f x 1 , . . . , f x n ∈ Q[t] such that we have the following one-to-one correspondence between V C (I) and V C (f ):

φ t : V C (I) ≈ V C (f ) := {β ∈ C | f (β) = 0} α -→ t(α) = a 1 α 1 + . . . + a n α n f x 1 (β) f 1 (β) , . . . , f x n (β) f 1 (β) ←-β.
(5)

If ( 5) is a RUR of V C (I), then the bijection φ t between the zeros of I and those of f preserves the multiplicities and the real zeros [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF]).

According to its definition, the computation of a Rational Univariate Representation can be divided in two independent parts. The first one consists in computing a separating linear form for V (I), that is a linear combination of the variables a 1 x 1 + . . . + a n x n so that the following map

V (I) -→ C (α 1 , . . . , α n ) -→ a 1 α 1 + . . . + a n α n
is injective. The second part consists in computing the univariate polynomials f, f 1 , f x 1 , . . . , f x n that define the one-to-one mapping given in (5).

If we suppose that a separating linear form t := a 1 x 1 + . . . + a n x n is given, the polynomials of the RUR can be computed by means of simple linear algebra operations. Indeed, provided that I is a zero-dimensional ideal, the quotient algebra A := Q[x 1 , . . . , x n ]/I is then a finite-dimensional Q-vector space which dimension is equal to the number of complex zeros of

V C (I) = {α = (α 1 , . . . , α n ) ∈ C n | ∀ P ∈ I : P (α) = 0}
counted with multiplicities, and using Gröbner bases (see, e.g., [START_REF] Cox | Varieties, and Algorithms[END_REF]), one can compute a basis of the Q-vector space A as well as the matrix associated with the Q-endomorphism of A defined by multiplying any element of A by an element P ∈ Q[x 1 , . . . , x n ] as follows

M P : A -→ A π(a) -→ π(P a), (6)
where π : Q[x 1 , . . . , x n ] -→ A denotes the Q-algebra epimorphism which sends a to its residue class π(a) in A. Then, the first polynomial of the RUR f is the characteristic polynomial of M t where t = a 1 x 1 + . . . + a n x n is the separating form for V (I). Moreover, if

f (t) := f gcd(f, ∂f /∂t) = d i=0 a i t d-i ∈ Q[t]
denotes the square-free part of f and if we note

H j (t) := j i=0 a i t j-i , j = 0, . . . , d -1,
then, the other polynomials of the RUR are defined as follow:

           f 0 := d-1 i=0 Trace(M i t ) H d-i-1 (t), f x k := d-1 i=0 M x k M i t H d-i-1 (t), k = 1, . . . , n.
As mentioned above, the computation of the RUR requires to find a separating linear form. Such a separating form can be sought in a set of infinite number of linear forms that contains only a finite number of bad-separating forms. For instance, the set {x

1 + i x 2 + • • • + i n-1 x n | i ∈ Q} contains at most (n -1) d(d-1)
2 bad separating forms (see (Basu et al., 2006, Lemma 4.9)). A random choice of a linear form then yields a form which is separating with probability close to one. Therefore, a standard strategy consists in choosing an arbitrary linear form, computing a RUR candidate using this form, and then checking whether or not this RUR candidate is a RUR. A well-known fact is that checking that a RUR candidate is a RUR, and thus that the chosen linear form is separating, is closely related to the problem of computing the number of distinct complex solutions of the system. Indeed, given the number of distinct solutions of the system, one can check that a linear form t is separating by checking whether or not the squarefree part of the first polynomial of the candidate RUR, i.e., f , has a degree equal to this number. In practice, several strategies can be used to compute the number of distinct complex solutions (e.g., computing the radical ideal of the polynomial ideal defined by the system, computing the rank of a well-chosen quadratic form [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF])). To avoid costly computations resulting from many trials, the separation test is usually performed modulo a prime number, which allows one to obtain a good prediction [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF]).

Finally, note that alternative algorithms exist for the computation of univariate representations which do not require the pre-computation of a Gröbner basis (see, e.g., [START_REF] Giusti | A Gröbner free alternative for solving polynomial systems[END_REF]). In addition, for the specific case of polynomial systems with only two variables, univariate representations as well as separating forms can be efficiently obtained using algorithms based on resultants and subresultants sequence (see Bouzidi (2014a)). For more details, see Section 3.

Once a RUR ( 5) is known, computing the real solutions of

I, namely comput- ing V R (I) = {(x 1 , . . . , x n ) ∈ R n | ∀ P ∈ I : P (x 1 , . . . , x n ) = 0} (resp.
, deciding whether or not the polynomial system defined by I has real solutions) reduces to computing the real roots of the univariate polynomial f 1 ∈ Q[t] (resp., deciding whether or not f 1 has real roots). This can be done using classical bisection algorithms such as Sturm's sequences or methods based on Descartes' rule of signs (see, e.g., [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]) which gives a set of intervals which isolate the real roots of f 1 .

Example 5 Let us illustrate the concept of univariate representation with the polynomial system which encodes the critical points of P = x 2 1 + x 2 2 -1 = 0 as defined in Example 4, namely:

   x 2 1 + x 2 2 -1 = 0, ∂P ∂x 2 = 2 x 2 = 0. ( 7
)
We note that (7) forms a Gröbner basis of the ideal I = P, ∂P ∂x 2 for the graded reverse lexicographic order (see, e.g., [START_REF] Cox | Varieties, and Algorithms[END_REF]) and a basis of the Q-vector space A = Q[x 1 , x 2 ]/I is given by the monomials {1, x 1 }, which implies that the number of complex solutions counted with their multiplicities is equals to two. Since the univariate polynomial in x 2 i.e., 2 x 2 , has degree one, t = x 2 is not a separating form. We can obtain a RUR of ( 7) by computing a Gröbner basis of I for a monomial ordering that eliminates x 2 such as:

x 2 1 -1 = 0, x 2 = 0.
Example 6 Another example is given by the following polynomial system:

P 1 = 36 x 2 1 x 2 -10 x 1 -6 x 2 = 0, P 2 = 12 x 2 1 + 30 x 1 x 2 -2 = 0. (8)
Computing a Gröbner basis of the ideal I = P 1 , P 2 generated by ( 8) for the graded reverse lexicographic order, we obtain that I = 9 x 2 2 + 1, 6 x 2 1 + 15 x 1 x 2 -1 , which shows that (8) is equivalently defined by these two polynomials. On the above Gröbner basis of I, we can then read that the dimension of the Q-vector space A = Q[x 1 , x 2 ]/I is 4: the set of monomials {1, x 1 , x 2 , x 1 x 2 } does not satisfy Qlinear relations which are algebraic consequences of (8). As a consequence, (8) admits 4 complex solutions counted with their multiplicities. Let us solve (8). Computing a Gröbner basis of I for an order which eliminates x 2 [START_REF] Cox | Varieties, and Algorithms[END_REF]), we get that the ideal I = 36 x 4 1 + 13 x 2 1 + 1, 36 x 3 1 + 19 x 1 + 15 x 2 . Hence, (8) can be parametrized as follows:

   x 2 = - 1 15 x 1 (36 x 2 1 + 19), 36 x 4 1 + 13 x 2 1 + 1 = 0.
Solving the last univariate equation, we obtain that the four solutions of ( 8) are of the form of:

x 1 , x 2 = - 1 15 x 1 (36 x 1 2 + 19) , x 1 ∈ ± i 2 , ± i 3 .
In particular, (8) does not admit real solutions, a fact that could be directly checked by applying Descartes' rule of signs on the univariate polynomial 36 x 4 1 + 13 x 2 1 + 1.

Real algebraic stability conditions

As shown at the beginning of Section 2.1, (2) can be reduced to checking the emptiness of the semi-algebraic set of R 2n defined by (3). This problem can be achieved using the methods described in Section 2.1.2. However, such an approach presents an important drawback: it doubles the number of variables (see (3)), which yields an overhead of computations in practice due to the exponential cost of the methods described in Sections 2.1.1 and 2.1.2.

To avoid this computational issue, we can start directly with the DeCarlo et al.'s stability conditions (see Theorem 2). From the computational point of view, the first n conditions of Theorem 2 can easily be checked using classical univariate stability tests (see, e.g., [START_REF] Marden | The geometry of the zeros of a polynomial in the complex plane[END_REF]; [START_REF] Jury | Theory and Application of the z-Transform Method[END_REF]; [START_REF] Bistritz | Zero location with respect to the unit circle of discrete-time linear system polynomials[END_REF][START_REF] Bistritz | Zero location of polynomials with respect to the unit-circle unhampered by nonessential singularities[END_REF]). We are then left with the last condition of Theorem 2, i.e.:

V C (D) ∩ T n = ∅. ( 9
)
This condition replaces the search for zeros of D in the unit polydisc U n (see (2)) by the search for zeros over the unit polycircle T n . Now, our approach to test the above condition (9) consists in applying a transformation that maps the unit poly-circle T n to the real space R n . More precisely, for each complex variable z k , we perform a change of variable

z k := φ(x k ) such that z k ∈ T if and only if x k ∈ R.
In particular, such a transformation allows us to keep unchanged the number of variables. A classical transformation that satisfies the above requirement is the so-called Möbius transformation which definition is recalled in the next definition.

Definition 1 Denoting the extended complex plane by C, namely, C := C ∪ {∞}, a Möbius transformation is the following rational function

φ : C -→ C z -→ a z + b c z + d ,
where a, b, c, d ∈ C are fixed and satisfy a d -b c = 0. We formally write:

φ - d c = ∞, φ(∞) = a c .
Denoting by H the class of circles of arbitrary radius in C (this class includes lines which can be considered as circles of infinite radius). Then, the set of Möbius transformations have the property of mapping H to itself, i.e., each circle in C is mapped to another circle in C. In particular, the following transformation φ(z) = z -i z + i , which corresponds to the Möbius transformation with a = 1, b = -i, c = 1 and d = i, maps the real line R to the complex unit circle deprived from the point 1, that is, T \ {1}.

Remark 2 Different transformations such as the classical parametrization of T \ {-i} defined by

∀ t k ∈ R, x k = 2 t k 1 + t 2 k , y k = 1 -t 2 k 1 + t 2 k , k = 1, . . . , n,
with the notation z k = x k + i y k , also fulfill the above requirement but usually yield a polynomial with higher degree than the one obtained by a Möbius transformation.

Accordingly, the following result holds.

Proposition 1 Let D ∈ R[z 1 , . . . , z n ].
Two polynomials R and I can be obtained such that:

V C (D) ∩ [T \ {1}] n = ∅ ⇐⇒ V R ( R, I ) = ∅.
Proof Given a polynomial D(z 1 , . . . , z n ), we can handle the following substitution

z k = x k -i x k + i = x 2 k -1 x 2 k + 1 -i 2 x k x 2 k + 1 , k = 1, . . . , n, (10)
which yields a rational function in C(x 1 , . . . , x n ) whose numerator can be written as

R(x 1 , . . . , x n ) + i I(x 1 , . . . , x n ),
where

R, I ∈ R[x 1 , . . . , x n ]. If (z 1 , . . . , z n ) is a zero of D(z 1 , . . . , z n ) that belongs to [T \ {1}] n , then, by construction, (x 1 , . . . , x n ) = i (z 1 + 1) 1 -z 1 , . . . , i (z n + 1) 1 -z n is a real zero of R(x 1 , . . . , x n ) + i I(x 1 , . . . , x n ). Conversely, if (x 1 , . . . , x n ) is a real zero of R(x 1 , . . . , x n ) + i I(x 1 , . . . , x n ), then (z 1 , . . . , z n ) = x 1 -i x 1 + i , . . . , x n -i x n + i is also a zero of D(z 1 , . . . , z n ) that belongs to [T \ {1}] n .
Remark 3 If we denote by d k the degree of D with respect to the variable z k , then we can easily check that the transformation used in the proof of Proposition 1 (see ( 10)) yields two polynomials R and I of total degrees at most n k=1 d k .

Example 7 Let us consider D(z 1 , z 2 ) = (2 z 2 1 + 10 z 1 + 12) + (z 2 1 + 5 z 1 + 6) z 2 . Applying the transformation (10) for k = 1, 2, we obtain the following polynomial system of total degree 3

R = 36 x 2 1 x 2 -10 x 1 -6 x 2 = 0, I = 12 x 2 1 + 30 x 1 x 2 -2 = 0,
which was considered in (8).

We can also consider

D(z 1 , z 2 ) = -3 z 2 1 z 2 2 + 2 z 2 1 z 2 + 2 z 1 z 2 2 -3 z 2 1 + 4 z 1 z 2 -3 z 2 2 + 2 z 1 + 2 z 2 -3
which, after transformation (10) for k = 1, 2, yields the following two polynomials of total degree 2: R = 0,

I = x 2 1 + x 2 2 -1 = 0.
According to Proposition 1, we can test that a polynomial D ∈ R[z 1 , . . . , z n ] does not have complex zeros in [T \ {1}] n by first computing the polynomials R(x 1 , . . . , x n ) and I(x 1 , . . . , x n ) and then checking that the following polynomial system R(x 1 , . . . , x n ) = 0,

I(x 1 , . . . , x n ) = 0, (11)
does not have any solution in R n by means of the techniques described in Section 2.1. Note that to check the last condition of Theorem 2, the above test is not sufficient since it excludes the points of the poly-circle that have at least one coordinate equal to 1. Hence, we also have to check that the polynomial D does not vanish at any of these points. Let us explain how this can be done in a systematic manner. Starting with D, we first compute the following polynomials:

D i (z 1 , . . . , z k-1 , z k+1 , . . . , z n ) := D(z 1 , . . . , z k-1 , 1, z k+1 , . . . , z n ), k = 1, . . . , n.
To each D k , we then apply the Möbius transformation (10) for z j with j = 1, . . . , k -1, k + 1, . . . , n, followed by the test of Proposition 1. Similarly as above, this test allows us to check whether or not each D k does not have complex zeros on [T \ {1}] n-1 . But we still need to check that D k does not vanish at the excluded points, namely, points having at least one coordinate in {z 1 , . . . , z k-1 , z k+1 , . . . , z n } equal to 1. This can then be done in the same way as above by considering the polynomials D k l obtained by substituting the variable z l by in the D k 's. Proceeding recursively until obtaining univariate polynomials of the form D(1, . . . , 1, z k , 1, . . . , 1), we can then check that D does not vanish on the unit poly-circles T.

Note that at the step m of the above process, the set of polynomials we have to consider are exactly the polynomials obtained from D after substituting m of the n variables z i 's by 1. From this observation, we obtain the following algorithm to check (10) based on Theorem 2.

Algorithm 1 1: procedure IntersectionEmpty(D(z 1 , . . . , z n )) return true if D satisfies (9)

2:

for k = 0 to n -1 do 3:

Compute the set S k of polynomials obtained by substituting k of the variables z i 's by 1 in D(z 1 , . . . , z n ) 4: return IntersectionEmpty(D(z 1 , . . . , z n )) 8: end procedure Remark 4 In Algorithm 1, the polynomials considered at the last step are D(1, . . . , z k , . . . , 1) for k = 1, . . . , n. Since the stability of these polynomials are checked at step 1 of Algorithm 2, we can skip this test in Algorithm 2 by stopping Algorithm 1 at step n -2.

for each D k in S k do 5: {R, I} = Möbius transform(D k ) 6: if V R ({R, I}) = ∅
Remark 5 From the computational cost viewpoint, it should be stressed that the most dominant part of Algorithm 1 is the first iteration of the outer loop which consists in checking that the polynomial D(z 1 , . . . , z n ) is devoid from zero in the poly-circle T n . Indeed, for any iteration k of the outer loop, the algorithm performs n k calls to the routine for checking the existence of real zeros of an algebraic system with n -k variables. Since this checking step requires a cost that is at least single exponential in the number of variables (see [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]), this implies that the cost of the outer loop, and thus of the whole algorithm, is dominated by the cost of the first iteration.

Examples and experiments

Let us illustrate Algorithms 1 and 2 with some explicit examples.

Example 8 We consider D(z 1 , z 2 ) = (2 z 2 1 + 10 z 1 + 12) + (z 2 1 + 5 z 1 + 6) z 2 which appears in several articles on the stability analysis [START_REF] Xu | Further improvements on Bose's 2-D stability test[END_REF] [START_REF] Li | Stability and stabilisation of linear multidimensional discrete systems in the frequency domain[END_REF]). It is known that D is structurally stable. Let us check again this result.

The first step of our procedure consists in checking that the two polynomials D(z 1 , 1) = 3 z 2 1 + 15 z 1 + 18 and D(1, z 2 ) = 12 z 2 + 24 are stable, which can be directly checked by, e.g., inspecting their solutions (i.e., {-3, -2} and {-2}).

In a second step, we apply Algorithm 1 to D(z 1 , z 2 ). As we have already checked that D(z 1 , 1) and D(1, z 2 ) are stable, we only have to consider D(z 1 , z 2 ) itself. Using the Möbius transformation (10), this polynomial yields the polynomial system defined by ( 8). In Example 6, we proved that (8) does not admit real solutions.

Example 9 If we consider

D(z 1 , z 2 ) = -3 z 2 1 z 2 2 + 2 z 2 1 z 2 + 2 z 1 z 2 2 -3 z 2 1 + 4 z 1 z 2 -3 z 2 2 + 2 z 1 + 2 z 2 -3,
then, the Möbius transformation (10) yields only one polynomial x 2 1 + x 2 2 -1 (see Example 7) that admits an infinite number of zeros. Checking for the existence of real zeros of this polynomial can be done by checking for the existence of real solutions for the system of its critical points (see ( 7)).

Example 10 We consider D(z 1 , z 2 , z 3 ) = (z 2 1 + z 2 2 + 4) (z 1 + z 2 + z 3 + 5)
which is known to be structurally stable [START_REF] Li | Stability and stabilisation of linear multidimensional discrete systems in the frequency domain[END_REF]. Our procedure first checks that the polynomials D(z 1 , 1, 1) = (z2 1 + 5) (z 1 + 7), D(1, z 2 , 1) = (z 2 2 + 5) (z 2 + 7) and D(1, 1, z 3 ) = 6 z 3 + 42 are stable. Then, applying Algorithm 1 to D, we have to test the existence of zeros on the polycircle of the following polynomials

         D(z 1 , z 2 , z 3 ) = (z 2 1 + z 2 2 + 4) (z 1 + z 2 + z 3 + 5), D(z 1 , 1, z 3 ) = (z 2 1 + 5) (z 1 + z 3 + 6), D(1, z 2 , z 3 ) = (z 2 2 + 5) (z 2 + z 3 + 6), D(z 1 , z 2 , 1) = (z 2 1 + z 2 2 + 4) (z 1 + z 2 + 6), (12)
by considering the set of polynomial systems obtained by applying the Möbius transformation (10) to each of them. The main difficult computation is to decide whether or not the following polynomial system, which corresponds to D(z 1 , z 2 , z 3 ), has real solutions:

                     48 x 3 1 x 3 2 x 3 -72 x 3 1 x 2 2 -96 x 3 1 x 2 x 3 -72 x 2 1 x 3 2 -184 x 2 1 x 2 2 x 3 -96 x 1 x 3 2 x 3 +24x 3 1 + 120 x 2 1 x 2 + 72 x 2 1 x 3 + 120 x 1 x 2 2 + 176 x 1 x 2 x 3 + 24 x 3 2 + 72 x 2 2 x 3 -40 x 1 -40 x 2 -24 x 3 = 0, 36 x 3 1 x 3 2 + 100 x 3 1 x 2 2 x 3 + 100 x 2 1 x 3 2 x 3 -68 x 3 1 x 2 -36 x 3 1 x 3 -124 x 2 1 x 2 2 -180 x 2 1 x 2 x 3 -68 x 1 x 3 2 -180 x 1 x 2 2 x 3 -36 x 3 2 x 3 + 44 x 2 1 + 116 x 1 x 2 +68 x 1 x 3 + 44 x 2 2 + 68 x 2 x 3 -12 = 0.
Using one of the methods described in Section 2.1.1, we can check that the above polynomial system does not have real zeros. Similarly, the second, third and fourth polynomials of ( 12) can be proved to be devoid from zeros in the corresponding polycircle, and we find again that D is stable.

Our stability test was implemented in a Maple routine named IsStable. It takes a polynomial defining the denominator of a transfer function in input and returns true if this polynomial satisfies (9) and false otherwise. For testing the first n conditions of Theorem 2, we use the classical 1-D Bistritz test (see [START_REF] Bistritz | Zero location with respect to the unit circle of discrete-time linear system polynomials[END_REF]) that was also implemented in Maple. To test the emptiness of a real algebraic set, which is the main critical step in Algorithm 1, we have implemented the two presented methods1 . The first one uses the classical cylindrical algebraic decomposition. Such a decomposition is provided by the Maple routine CylindricalAlgebraicDecompose which can be found in the native package RegularChains [SemiAlgebraicSetTools]. The second method is based on the computation of the set of critical points of a given function restricted to the real algebraic set under consideration. An efficient implementation of this method has been done is the external Maple library RAGlib (Safey El Din ( 2007)) (see the command HasRealSolutions). Finally, we use the routine Isolate of the Maple package RootFinding in order to compute numerical approximations of the solutions through the computation of a univariate representation.

In Table 1, we show the average running times in CPU seconds of the IsStable routine for random (sparse or dense) polynomials in 2 and 3 variables with rational coefficients 2 . The two running time columns correspond to the two variants IsStableCAD and IsStableCRIT (depending on the method used for testing the emptiness of a real algebraic set) of the routine IsStable. Remark 6 Note that Algorithm 2 can check the structural stability of polynomials in 4 variables with degree up to 12 in less than 20 minutes. However, when the polynomials have more variables (i.e., larger than 4) or have larger degrees, these methods do not return a result in a reasonable time.

Data

A stability test for 2-D systems

In Section 2, a general framework was proposed for the stability analysis of n-D systems with n ≥ 2. In this section, we restrict the study to the particular case of n = 2 and we show that substantial improvements with respect to practical efficiency can be obtained by using state-of-the-art algorithms developed in Bouzidi (2014b) for the computation of the solutions of bivariate algebraic systems of equations.

Recall that testing the stability of a 2-D system can be reduced to deciding whether an algebraic system of the form of {R(x 1 , x 2 ) = I(x 1 , x 2 ) = 0} admits real solutions. In the present case, without loss of generality, we can assume that the ideal I generated by the two polynomials R and I is zero-dimensional, i.e., gcd(R, I) = 1 3 . Our contribution in this section is a dedicated method for deciding

3 If gcd(R, I) is non-trivial, then it is sufficient to compute their gcd G in Q[x 1 , x 2 ] and
to consider the two zero-dimensional systems: the system R G , I G and the system G, ∂G ∂x 1 which encodes the critical points of G with respect to x 1 .

if a system of two polynomial equations in two variables, having a finite number of complex solutions, admits real solutions.

Since the ideal I is zero-dimensional, we can use the univariate representation techniques described in Section 2.1.2 to reduce the problem of deciding the existence of real solutions of I to that of deciding the existence of real roots of a univariate polynomial. Indeed, under the hypothesis that I is zero-dimensional, the quotient algebra A := Q[x 1 , x 2 ]/I inherits a finite-dimensional Q-vector space structure and, for any polynomial P ∈ Q[x 1 , x 2 ], we can define the Q-endomorphism M P of A defined by ( 6) whose the characteristic polynomial can be written as follows

C P (t) = det(t I d -M P ) = (α 1 , α 2 )∈V (I) (t -P (α 1 , α 2 )) µ(α 1 , α 2 ) , (13)
where I d denotes the identity matrix of size d = dim Q (A) and µ(α 1 , α 2 ) the multiplicity of (α 1 , α 2 ) as a zero of I. Furthermore, if P = x 1 +a x 2 is a separating form of V (I) 4 , then the polynomial C P coincides with the polynomial f of the Rational Univariate Representation of I computed with respect to P (see ( 5) and Section 2.1.2) which yields an important property regarding to the existence of real solutions of V (I) and of C P . The following result can be proved considering a univariate representation of the solutions (see, e.g., [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF]).

Theorem 3 Let P = x 1 +a x 2 be a separating form for V (I). Then, the univariate polynomial C P has real roots if and only if V (I) has real solutions.

Consequently, the computation of a separating form of V (I) and the corresponding polynomial C P reduces the problem of searching for real solutions of V (I) to the problem of searching for real roots of C P . In the following, instead of the classical strategy which requires computations in the quotient algebra A := Q[x 1 , x 2 ]/I, we propose an alternative approach based on the computation of the so-called resultant polynomial as well as the computation of a generic position. Thus, before going further, let us introduce the concept of resultant and subresultant sequences, as well as some of their basic properties which are useful for the description of our algorithm.

Resultant and subresultant sequence

Let A be a unique factorization domain [START_REF] Cox | Varieties, and Algorithms[END_REF]), e.g.,

A := K[y], where K is a field. Let f = p i=0 a i x i ∈ A[x] and g = q j=0 b j x j ∈ A[x]
, that is, the a i 's and b j 's belong to A. Let us suppose that a p = 0 and b q = 0 so that deg x f = p and deg x g = q, and p ≥ q. Let A[x] n = {P ∈ A[x] | deg x P ≤ n} be the set of polynomials with degree at most n and {x i } i=0,...,n the standard basis of the free 4 It has been shown that a separating form can be sought in the following set of linear forms d (d-1) 2

x 1 + a x 2 | a = 0, . . . ,
, where d denotes the cardinal of V (I), i.e., the dimension dim Q (A) of A as a Q-vector space. Indeed, such a set contains at least one separating form since the number of non-separating forms is bounded by the number

d (d-1) 2
of lines passing through two distinct points among d points in the plane.

A-module A[x] n of rank n. We set A[x] n = 0 for negative integer n. For 0 ≤ k ≤ q, we can consider the following homomorphism of free A-modules:

ϕ k : A[x] q-k-1 × A[x] p-k-1 -→ A[x] p+q-k-1 (U, V ) -→ U f + V g.

Using the standard basis of

A[x] q-k-1 (resp., A[x] p-k-1 , A[x] p+q-k-1
) and identifying the polynomial q-k-1 i=0 u i x i ∈ A[x] q-k-1 with the row vector (u 0 , . . . , u q-k-1 ) ∈ A 1×(p-k) , we obtain that ϕ k (u 0 , . . . , u q-k-1 , v 0 , . . . , v p-k-1 ) = (u 0 , . . . , u q-k-1 , v 0 , . . . , v p-k-1 ) S k , where the matrix S k is the matrix defined by: p+q-k) .

S k = U k V k ∈ A (q-k+p-k)×(p+q-k) , U k =      a 0 a 1 . . . a q-k . . . a p 0 . . . 0 0 a 0 . . . a q-k-1 . . . a p-
     ∈ A (q-k)×(p+q-k) , V k =      b 0 b 1 . . . b p-k . . . b q 0 . . . 0 0 b 0 . . . b q-k-1 . . . b q-1 b q . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 . . . 0 b 0 . . . . . . . . . . . . b q      ∈ A (p-k)×(
To be coherent with the degree of polynomials, we attach index i -1 to the i th column of S k so that the index of the columns goes from 0 to p + q -k -1.

Definition 2 For 0 ≤ j ≤ p+q -k -1 and 0 ≤ k ≤ q, let sr k,j be the determinant of the submatrix of S k formed by the last p+q-2 k-1 columns, the column of index j and all the p + q -2 k rows. The polynomial Sres k (f, g) = sr k,k x k + . . . + sr k,0 is the k th subresultant of f and g, and its leading term sr k,k is the k th principal subresultant of f and g. The matrix S 0 ∈ A (p+q)×(p+q) is the Sylvester matrix associated with f and g, and Res x (f, g) = det S 0 is the resultant of f and g.

Remark 7 For k < j ≤ p + q -2 k -1, we note that sr k,j = 0 since sr k,j is the determinant of a matrix having twice the same column. Moreover, we can check that we have: sr q,q = b p-q q , ∀ q < p, Sres q (f, g) = b p-q-1 q g. Since A is an integral domain, we can consider its field of fractions which we denotes by F, namely,

F := { n d | 0 = d, n ∈ A}, and the Euclidean domain F[x]. Since f, g ∈ F[x],
we can define the greatest common factor gcd(f, g), which is defined up to a non-zero element of F, so that we can suppose that gcd(f, g) ∈ A.

Theorem 4 [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]) The first Sres k (f, g) such that sr k,k = 0 is equal to gcd(f, g). Now, if we consider two polynomials in two variables f = p i=0 a i (x 1 ) x i 2 and g = q j=0 b j (x 1 ) x j 2 , and x = x 2 so that A = K[x 1 ] and A[x] = K[x 1 , x 2 ], then we have the important results.

Theorem 5 [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]

) Let f (x 1 , x 2 ), g(x 1 , x 2 ) ∈ K[x 1 , x 2 ] be two bivariate polynomials.
-The roots of Res x 2 (f, g)(x 1 ) are the projection onto the x 1 -axis of the common solutions of f and g and of the common roots of a p and b q . -For any α root of Res x 2 (f, g) such that a p (α) and b q (α) do not both vanish, the first polynomial Sres k (α, x 2 ), for increasing k, that does not identically vanish is of degree k and is the gcd of f (α, x 2 ) and g(α, x 2 ), up to a nonzero constant in the field of fractions of A(α).

The subresultant sequences can be computed either by means of determinant computations or by applying a variant of the classical Euclidean algorithm (see [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]). The latter method, combined with evaluation/interpolation strategies, turns out to be much more efficient in practice, especially for the case of univariate or bivariate polynomials.

Computation of a separating form of bivariate polynomial systems

Given a linear form x 1 + a x 2 (not necessarily separating), the following theorem shows that, up to a non-zero factor in Q, the univariate polynomial C x 1 +a x 2 (see ( 13)) is equal to the resultant of the two polynomials obtained by applying a change of variables to R and I.

Theorem 6 (Bouzidi et al. (2015b)) Let R, I ∈ Q[x 1 , x 2 ] and let us define the polynomials

R (t, x 2 ) := R(t -a x 2 , x 2 ), I (t, x 2 ) := I(t -a x 2 , x 2 ), ( 14
)
where a ∈ Z is such that the leading coefficient of R and I with respect to x 2 are coprime. Then, the resultant of R and I with respect to x 2 , denoted by

Res x 2 (R , C ), is equal to: c (α 1 , α 2 )∈V (I) (t -α 1 -a α 2 ) µ(α 1 ,α 2 ) , c ∈ Q \ {0}.
In practice, the computation of C x 1 +a x 2 as a resultant (see Theorem 6) is much more efficient than computing the characteristic polynomial of the Q-endomorphism M x 1 +a x 2 (see ( 6)) since the computation of the matrix M x 1 +a x 2 usually requires the costly pre-computation of a Gröbner basis of the ideal I = R, I for the graded reverse lexicographic order.

Let us now focus on the computation of a separating form for V (I). Below, we propose a method that consists in applying a change of variables to the initial system and then, using resultant and subresultants, to check whether or not the resulting system is in generic position as defined below.

Definition 3 Let f (x 1 , x 2 ), g(x 1 , x 2 ) ∈ Q[x 1 , x 2 ].
If S denotes the cardinality of a finite set S, then the system {f, g} is said to be in generic position with respect to x 1 if we have:

∀ α ∈ C, {β ∈ C | f (α, β) = g(α, β) = 0} ≤ 1.
Let us first illustrate our approach with an example.

Example 11 Consider the ideal

I = {f = x 2 -x 2 1 , g = x 2 1 + x 2 2 -2}
whose set of solutions V (I) consists in four points of C 2 . The resultant Res x 2 (f, g) of f and g with respect to x 2 is equal to x 4 1 + x 2 1 -2. The roots of Res x 2 (f, g) correspond to the projections of the four solutions of V (I) onto the x 1 -axis. Since all these roots are distinct, x 1 is a separating form (see Figure 1 for the real solutions of V (I)). The fact that the solutions of V (I) project distinctly onto x 1 can be algebraically described by the fact that for each root α of Res x 2 (f, g), the gcd

x 2 + α 2 of f (α, x 2 ) = x 2 -α 2 and g(α, x 2 ) = x 2
2 + α 2 -1 has only one root. Let us now consider the ideal

I = {f = (x 1 -2) 2 +x 2 2 -2, g = x 2 1 +x 2 2 -2}.
The resultant Res x 2 (f , g ) of f and g with respect to x 2 , namely, 16 (x 1 -1) 2 , has a single (real) root 1 of multiplicity 2, and gcd(f (1, x 2 ), g (1, x 2 )) = x 2 2 -1 admits two distinct roots -1 and 1 which correspond to two different solutions of V (I ). This means that the system is not in generic position, and thus that x 1 is not a separating form (see Figure 2). In order to compute a separating form for V (I ), we can apply a change of variables to f and g , for instance t = x 1 + x 2 , and then compute the resultant of these new polynomials f (t -x 2 , x 2 ) and g (t -x 2 , x 2 ) with respect to x 2 . This yields the polynomial t (t -2) whose two distinct roots {0, 2} are the projections of the solutions onto the t-axis (see Figure 3). For α 1 = 0 (resp., α 2 = 2), the gcd of f (-x 2 , x 2 ) and g (-x 2 , x 2 ) (resp., of f (2 -x 2 , x 2 ) and g (2 -x 2 , x 2 )) is x 2 + 1 (resp., x 2 -1). Since both gcds admit only one root, then the system {f (t -x 2 , x 2 ), g (t -x 2 , x 2 )} is in generic position with respect to t and thus x 1 + x 2 is a separating form for V (I ).

(x 1 = t -x 2 , x 2 )
Given a linear form t = x 1 + a x 2 , it can be shown that it is separating for V ( R, I ) if and only if the system {R , I } is in generic position (see Definition 3). Algebraically, this means that for each root α of Res x (R , I ) (where R and I are defined as in Theorem 6), the gcd of R (α, x) and I (α, x), denoted G(α, x), has exactly one distinct root.

To check the above genericity condition, we can consider a triangular description of the solutions of {R , I } given by a finite union of triangular systems:

V ( R , I ) = l k=1 (α, β) ∈ C 2 | r k (α) = G k (α, β) = 0 .
Such a triangular description can be obtained via a triangular decomposition algorithm based on the resultant and subresultant polynomials (see Algorithm 1 of Bouzidi et al. (2015b) for more details). Such a triangular decomposition yields a set of triangular systems of the form {r k (t), Sres k (t, x 2 )} k=1,...,l , where l = min{deg x 2 R , deg

x 2 I }, Res x 2 (R , I ) = l k=1 r k (t), r k ∈ K[t]
is the factor of Res x 2 (R , I ) (possibly equal to 1) whose roots α's satisfy the property that the degree of G(α, x 2 ) (i.e., gcd(R (α, x 2 ), I (α, x 2 ))) in x 2 is equal to k and Sres k (t, x 2 ) = k i=0 sr k,i (t) x i 2 is the k th subresultant of R and I . Once a triangular decomposition {r k (t), Sres k (t, x 2 )} k=1,...,l of {R , I } is computed, the genericity condition is equivalent to the fact that Sres k (t, x 2 ) can be written as (a k (t) x 2 -b k (t)) k modulo r k (t), with gcd(a k , b k ) = 1. The next theorem checks this last condition.

Theorem 7 [START_REF] Daouda | On the computation of the topology of a non-reduced implicit space curve[END_REF]

) Let R(x 1 , x 2 ), I(x 1 , x 2 ) ∈ Q[x 1 , x 2 ]. De- fine the polynomials R (t, x 2 ), I (t,
x 2 ) as in Theorem 6, and let {r k (t), Sres k (t, x 2 )} k=1,...,l be the triangular decomposition of {R , I }. Then, {R , I } is in generic position 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 , (15)
for all k ∈ {1, . . . , l} and for all i ∈ {0, . . . , k -1}.

Finally, our algorithm for checking whether the system {R, I} admits real solutions consists in computing the above triangular decomposition for the system {R , I } obtained after applying an arbitrary linear change of variable t = x 1 + a x 2 . If the condition of Theorem 7 is satisfied, then x 1 + a x 2 is a separating form. It then remains to check if the resultant Res x 2 (R , I ) of R and I with respect to x 2 has real roots, a fact which can be done using, for instance, Sturm sequences [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]).

Remark 8 In practice, several strategies are used in order to reduce the computational time of the above algorithm. For instance, the computation is stopped when the resultant, computed for some linear form x 1 + a x 2 , that is the resultant of R(t -a x 2 , x 2 ) and I(t -a x 2 , x 2 ) with respect to x 2 , is devoid from multiple factors, which implies that the form x 1 + a x 2 is separating by Theorem 6. The computation is also stopped when the computed resultant does not have real zeros, since it implies that the system does not have real zeros as well. Another improvement is the way we can choose the form x 1 + a x 2 candidate to be a separating form. Indeed, in order to increase the probability that a form is separating, a first computation is performed modulo a prime number p (coefficients are then considered in the finite field F p = Z/Z p). Such a computation turns out to be very fast since it avoids coefficient swells in the algorithm. Providing that a linear form is separating modulo a prime p, then, with high probability, it is also separating over Z and we can choose it as a candidate for the algorithm over Z.

Experiments

In order to measure the gain of our algorithm with respect to the general algorithm described in Section 2, we compare it with the general method Isolate partially developed by the same authors and available in the package RootFinding of the Maple computer algebra system. This function first computes a Rational Univariate Representation [START_REF] Rouillier | Solving zero-dimensional systems through the rational univariate representation[END_REF]) from a Gröbner basis computed with the F 4 algorithm [START_REF] Faugère | A new efficient algorithm for computing Gröbner bases (F 4 )[END_REF]), and then uses of a variant of Descartes algorithm [START_REF] Rouillier | Efficient isolation of polynomial real roots[END_REF]) as well as multi-precision interval arithmetic [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the mpfi library[END_REF]) to isolate the real roots of the system.

For the present experiments, we re-use black boxes developed for the algorithms described in Bouzidi et al. (2015a[START_REF] Bouzidi | Solving bivariate systems: Efficient randomized algorithm for computing rational univariate representations and application to the topology of plane curves[END_REF]; Bouzidi (2014a) which use exactly the same technical base to design the component of the algorithm that computes the univariate polynomial C x 1 +a x 2 and performs the separation check. All the other components are shared with the Maple RootFinding[Isolate] function.

For dense polynomials with coefficients that can be encoded on 23 bits (such as if there were coming from floating point numbers), the results -obtained on a core i7 3.5 Ghz with 32 GB of memory -are summarized in the following table in which Degree denotes the total degree of the polynomial D(z 1 , z 2 ) to be studied, V (I) the number of complex roots of the bivariate system to be solved to decide stability, RootF inding the computational time of the function RootFinding[Isolate] and Dedicated the computational time of our new dedicated algorithm.

Degree

V For these examples, note that we did not report the computation times required for the two 1-D stability tests (i.e., the stability test for D(1, z 2 ) and D(z 1 , 1)) since they are small compared with the resolution of the bivariate polynomial system.

Finally, we point out that our naive implementation of the Möbius transform in Maple is the main bottleneck of our dedicated algorithm compare with the extremely efficient algorithm for the real solution computation of systems of two polynomials in two variables.

A stability test for 2-D systems with parameters

In what follows, let u = {u 1 , . . . , u m } denote a set of parameters. Throughout this section, these parameters are assumed to be real (for complex parameters, see Remark 10 at the end of the present section). In this last section, we study the structural stability of 2-D systems given by a transfer function of the form of (1) that depends on the set of parameters u, i.e., where n = 2 and D ∈ Q[z 1 , z 2 , u]5 . In other words, our goal is to study (2) in terms of the parameters u k 's. Roughly speaking, our approach consists in computing a set of polynomials {h 1 , . . . , h s }, where h k ∈ Q[u 1 , . . . , u r ] for k = 1, . . . , s, satisfying the property that the stability of (1) does not change provided that the sign of the sequence {h 1 , . . . , h s } does not change. Then, R m can be decomposed into cells in which the signs of {h 1 , . . . , h s } remain the same and the cells for which the system is structurally stable can then be selected by testing the structural stability of the system for a single parameter's value in each cell.

Considering D(z 1 , z 2 , u) as a polynomial in the variables z 1 and z 2 with coefficients in Q[u 1 , . . . , u r ], we can apply the transformation given in Section 2.2, which yields the following set of conditions:

     D(z 1 , 1, u) = 0, |z 1 | ≤ 1, D(1, z 2 , u) = 0, |z 2 | ≤ 1, V ( R(x 1 , x 2 , u), I(x 1 , x 2 , u) ) ∩ R 2 = ∅. (16)
We start with the study of the first two conditions involving univariate polynomials with parameters. We first transform these conditions so that continuous stability tests can be applied. More precisely, we can apply the following change of variables

s 1 = 1 -z 1 1 + z 1 , s 2 = 1 -z 2 1 + z 2
to the polynomials D(z 1 , 1, u) and D(1, z 2 , u). We denote by D 1 (s 1 , u) (resp.

D 2 (s 2 , u)), the numerator of D - 1 -s 1 1 + s 1 , 1, u (resp. D 1, - 1 -s 2 1 + s 2 , u ).
The first two conditions of ( 16) then become:

D 1 (s 1 , u) = 0, ∀ s 1 ∈ C : Re(s 1 ) ≥ 0, D 2 (s 2 , u) = 0, ∀ s 2 ∈ C : Re(s 2 ) ≥ 0.
Then, we can use a classical result of Liénard and Chipart ((Basu et al., 2006, Thm. 9.30)) that expresses the stability condition of a continuous polynomial D(s) as a positivity condition of its coefficients as well as a certain signed principal subresultant sequence of two polynomials F (s) and G(s) satisfying D(s) = F (s 2 )+ s G(s 2 ) (see (Basu et al., 2006, Thm. 9.30)). Using the specialization property of subresultants (see [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]), we can generalize this result to the case of univariate polynomials depending on parameters. In particular, applying this test to the polynomials D 1 (s 1 , u) and D 2 (s 2 , u) yields a set of polynomials depending only on the parameters u, and the stability of D 1 (s 1 , u) and D 2 (s 2 , u), and thus of D(z 1 , 1, u) and D(1, z 2 , u), is then satisfied providing that these polynomials {h i (u)} i=1,...,t are all positive.

The next problem is to decide whether or not the following system R(x 1 , x 2 , u) = 0,

I(x 1 , x 2 , u) = 0, (17)
admits real solutions. In what follows, we can assume that ( 17) is generically zerodimensional, that is, for almost all values of the parameters u ∈ R m , (17) admits a finite number of complex solutions. The main tool we use to solve this problem is the so-called discriminant variety, first introduced in [START_REF] Lazard | Solving parametric polynomial systems[END_REF], and recalled in the next section.

Discriminant varieties: definition and properties

Before recalling the definition of a discriminant variety of an algebraic set, let us start with some useful notations. For a set of polynomials p 1 , . . . , p l ∈ Q[x 1 , . . . , x n-m , u 1 , . . . , u m ], the corresponding algebraic set is defined as:

S = {α ∈ C n | p 1 (α) = 0, . . . , p l (α) = 0}. ( 18
)
We consider the canonical projection onto the parameter space C m , namely, the following map

Π u : C n -→ C m (x 1 , . . . , x n-m , u 1 , . . . , u m ) -→ (u 1 , . . . , u m ),
and we denote by Π u (S) the so-called Zariski closure of the projection of S onto the parameter space C m . For more details, see [START_REF] Cox | Varieties, and Algorithms[END_REF].

Definition 4 [START_REF] Lazard | Solving parametric polynomial systems[END_REF]) With the above notations, an algebraic set V of C m is called a discriminant variety of S if and only if:

1. V is contained in Π u (S).
2. The connected components U 1 , . . . , U s of Π u (S) \ V are analytic sub-manifolds (note that if Π u (S) is connected, there is only one component). 3. For j = 1, . . . , s, (Π -1 u (U j ) ∩ S, Π u ) is an analytic covering of U j .

In broad terms, a discriminant variety yields a partition of the parameter's space C m into cells U, such that for each cell, the cardinal of Π -1 u (µ) ∩ S, where µ ∈ C m , is locally constant on U, and Π -1 u (U) ∩ S consists of a finite collection of sheets which are all locally homeomorphic to U.

A consequence of Definition 4, stated in the following theorem, is a fundamental property of the discriminant variety regarding to the number of solutions. In this theorem, we assume that the polynomial system S defined by ( 18) is generically zero-dimensional, i.e., for almost all values of the parameters µ ∈ C m , the polynomial system S u=µ , obtained by substituting the parameters u to µ admits a finite number of complex solutions.

Theorem 8 [START_REF] Lazard | Solving parametric polynomial systems[END_REF]) Let S be an algebraic system and U 1 , . . . , U s defined as in Definition 4. Then, for two vectors of parameters µ, ν ∈ U j , the specialized polynomial systems S u=µ and S u=ν have exactly the same number of distinct zeros.

For a given set S, the smallest algebraic variety V that satisfies the conditions of Definition 4 is called the minimal discriminant variety (see [START_REF] Lazard | Solving parametric polynomial systems[END_REF]).

Example 12 A classical exemple of a discriminant variety is the zeros of the discriminant of a quadratic univariate polynomial f := a x 2 +b x+c whose coefficients are given as parameters. This discriminant is given as b 2 -4 a c and satisfies that for all a 0 , b 0 and c 0 such that b 2 0 -4 a 0 c 0 = 0, the polynomial a 0 x 2 + b 0 x + c 0 has exactly two distinct roots.

In the sequel, we simplify say "discriminant variety" for "minimal discriminant variety".

Computation of discriminant varieties

For a system F defined by a set {p 1 , . . . , p l } ⊂ Q[x 1 , . . . , x n-m , u 1 , . . . , u m ], by means of variable eliminations using, e.g. standard Gröbner bases computations (see, e.g., [START_REF] Lazard | Solving parametric polynomial systems[END_REF]), we can compute a sequence of polynomials {h 1 , . . . , h s } ⊂ Q[u 1 , . . . , u m ] whose zeros define the discriminant variety of F. For instance, in the case of the quadratic polynomial given in Example 12, the discriminant is computed by eliminating the variable x in the system defined by f and its derivative ∂f ∂x with respect to x, which can be done, e.g. by computing the resultant of f and ∂f ∂x with respect to x. In our setting, namely a system of two polynomial equations in two variables S = V ( R, I ), the discriminant variety, denoted by V D , consists in the union of the following two subsets (see [START_REF] Lazard | Solving parametric polynomial systems[END_REF] for details):

-The set O ∞ of α ∈ C m such that Π -1 u (U) ∩ S is not compact for any compact neighborhood U of α in C m .
-The set O c of the union of the critical values of Π u and of the projection of the singular points of V (S) onto C m .

Intuitively, O ∞ represents parameter values such that there exist either vertical leafs of solutions or leafs that go to infinity above some of their neighborhoods, while O c represents parameter values such that above some of their neighborhoods, the number of leafs varies. Thus, the minimal discriminant variety V D roughly represents parameter values over which the number of solutions of (17) changes. Furthermore, an important remark for the computation of the discriminant variety of S is that both O ∞ and O c are algebraic sets. V D can thus be described as the union of two algebraic sets that can be computed independently.

Both O ∞ and O c are projections of algebraic set. Computing these varieties remains to eliminating variables in the systems of equations corresponding to these varieties, which corresponds to the following problem: given

I = f 1 , . . . , f l ⊂ K[x 1 , x 2 , u], compute Π u (V (I)) = V (I u ), where I u ⊂ K[u] is defined by I u = I ∩ K[u].
Algorithmically, I u can be obtained by means of a Gröbner basis for any elimination ordering < satisfying u < x 1 , x 2 . More precisely, it suffices to compute a Gröbner basis for such an ordering and to keep only the polynomials that belong to K[u].

In [START_REF] Lazard | Solving parametric polynomial systems[END_REF], it was shown that the set O c is equal to Π u (V ( R, I, Jac x 1 ,x 2 (R, I) )), where Jac x 1 ,x 2 (R, I) denotes the determinant of the Jacobian matrix with respect to the variables x 1 and x 2 . Hence, computing an ideal I c such that O c = V (I c ) remains to computing the determinant Jac x 1 ,x 2 (R, I) and a Gröbner basis of the ideal R, I, Jac x 1 ,x 2 (R, I) for any elimination monomial ordering < satisfying u < x 1 , x 2 . It was also noticed that such an elimination ordering allows us to compute an ideal

I ∞ ⊂ Q[u] such that O ∞ = V (I ∞ ).
Precisely, suppose that G is a reduced Gröbner basis of R, I for a monomial ordering < u,x 1 ,x 2 , that is, the product of two degree reverse lexicographic orderings < u for the parameters and < x 1 ,x 2 for the variables. For more details, see [START_REF] Cox | Varieties, and Algorithms[END_REF]. Let us define the following ideal

I i ∞ = {LM < x 1 ,x 2 (g) | g ∈ G, ∃ m ≥ 0, 1 ≤ i ≤ 2, LM < x 1 ,x 2 (g) = x m i },
where LM < denotes the leading monomial of a polynomial for an admissible monomial ordering < (see [START_REF] Cox | Varieties, and Algorithms[END_REF]). Then, we have:

-I i ∞ ⊂ K[u] is a Gröbner basis for < u . -O ∞ = V (I 1 ∞ ) V (I 2 ∞ ) = V (I 1 ∞ I 2 ∞ ).

Discussing the number of real solutions

Once a discriminant variety V D of S = V ( R, I ), represented by a set of polynomials {h 1 , . . . , h s }, is computed, we can compute a CAD adapted to these polynomials (see Section 2.1) in order to obtain a partition of the parameter space R m defined by the discriminant variety V D and the connected components of its complementary R m \V D (which has the property that over any cell U that does not meet W D , Π -1 u (U) is an analytic covering of U). In particular, the number of zeros of S is constant over any connected set that does not intersect the discriminant variety.

Also, for computing the constant number of solutions over a connected component that does not meet the discriminant variety, it suffices to take a particular value of parameter values µ in this component and then solve the corresponding zero-dimensional polynomial system S u=µ .

Remark 9 Note that the structure of the solutions is not known above the discriminant variety itself. Since the discriminant variety is a set of null measure, it is useless here to study what is going on for such parameter values. However, the discriminant variety is defined by a polynomial system which can be added to the original system in order to follow the study recursively.

The discriminant variety is defined for the complex solutions of (17). For real solutions, only two cases may occur:

1. Π u (S ∩ R m+2 ) ⊂ V D . We then need to study S ∩ Π -1 u (V D ) instead of S. 2. Π u (S ∩ R m+2 ) V D . Then, we have V D ∩ R m is a discriminant variety for S ∩ R m+2 , which is the usual situation.

In the second case, note that if V D is minimal for S, then V D ∩ R m is not necessarily minimal for S ∩ R m+2 .

Computing stability regions

Let us now come back to the original problem, namely the computation of the regions of the parameter space such that (16) are satisfied, and thus define stable systems. As mentioned at the beginning of Section 4, we can compute a set of polynomials {p i (u)} i=1,...,t such that the first two conditions of ( 16) are satisfied if and only if p i (u) > 0 for i = 1, . . . , t. We can also compute a set of polynomials {q j (u)} j=1,...,s that defines a partition of the parameter space in which the number of real solutions of ( 17) is constant. Now, considering the global set of polynomials F := {p 1 (u), . . . , p t (u), q 1 (u), . . . , q s (u)}, we can compute a CAD adapted to F (see Section 2.1). This yields a disjoint union of cells in R m in which the signs of all the polynomials of F (both p i 's and q i 's) are constant. In particular, inside each of these cells, both the sign of the polynomials p i 's and the number of real solutions of (17) are constant. This implies that the system is either stable or unstable. To determine the cells for which the system is stable, it suffices to select a simple point u = µ in each cell and to test (16) after the evaluation of the parameters.

Finally, in practice, to reduce the running time computation, we only compute the cells that have maximal dimension during the CAD.

Remark 10 In this section, the parameters were assumed to be real. The main reason is that the coefficients of a transfer function is usually assumed to be real in control theory. Moreover, the univariate stability test we use holds only for real coefficients and the Maple implementation of the CAD provides a partition of the space R m . But our method can be extended to the case of complex parameters: we can first use a univariate stability test that handles complex coefficients (see, e.g., [START_REF] Bistritz | A circular stability test for general polynomials[END_REF]) and then carry out a complex variant of the CAD, that is a partition of the space C m . For the latter, it is sufficient to consider the discriminant variety as an algebraic variety over the space R 2m (considering the real and imaginary parts of each parameter) and to perform the CAD associated to this Decomposing the parameter space cylindrically with respect to these 14 curves gives 1161 cells shown in Figure 4.

For parameters which belong to these cells, the system is either stable or unstable. To test the stability of the corresponding system, it is sufficient to test the stability of the system obtained by evaluating the parameters u to a numerical value µ in this cell and to count the number of real solutions of the (non paramet-ric) zero-dimensional polynomial system (17) defined with u = µ, and to perform the stability test of D(z 1 , 1) and D(1, z 2 ). It should be noticed that in some regions of the parameter space, some cells are very small.

Finally, it turns out that 118 of these regions correspond to unstable systems. For instance, the cell containing the pair (u 1 = -.4717912847, u 2 = -.5389591122) defines unstable systems while the cell containing the pair (u 1 = -.6152602220, u 2 = -.5389591122) as well as the cell containing the pair (u 1 = -.3942379536, u 2 = -.5389591122) define stable systems (see Figure 6).

Conclusion

The main goal of this paper was to point out some advantages of using classical techniques from the computer algebra community in the context of the stability analysis of multidimensional systems. Indeed, using state-of-the-art algorithms for solving algebraic systems of equations, several methods for the study of structural stability of these systems have been developed. The novelty of these methods compared to the existing ones is that they are both non-conservative and show promising results in practice especially for 2D and 3D systems. Moreover, despite of their own interests for testing the stability, the obtained algorithms can also be used for solving similar problems such as the computation of stabilizing feedback control for 1D linear systems or for the stabilization of n-D systems. From the computational point of view, we would try to improve the practical behavior of these methods in the case of n-D systems (n > 3) by investigating the use of numerical routines while keeping the exactness aspect of the approach since it is critical in our problems. This investigation will be the subject of further works. In addition, other classes of linear systems such as time-delay systems share the same type of representation, and can thus be addressed using the same computer algebra techniques.

  0, . . . , 0, z n ) = 0,|z n | ≤ 1, D(0, . . . , 0, z n-1 , z n ) = 0, |z n-1 | ≤ 1, |z n | = 1, z 2 , . . . , z n ) = 0, |z 2 | ≤ 1, |z i | = 1, i = 3, . . . , n, D(z 1 , z 2 , . . . , z n ) = 0, |z 1 | ≤ 1, |z i | = 1, i = 2, . . . , n.

Fig. 1

 1 Fig. 1 Intersection between a circle and a parabola

Fig. 2

 2 Fig. 2 Intersection between two circlesFig.3After the change of variables (x 1 = t -x 2 , x 2 )

Fig. 4

 4 Fig. 4 Global view -Parameter space decomposition

Fig. 6

 6 Fig. 6 Zoom around an unstable region : u 1 = -0.4 . . . -0.6, u 2 = -0.4 . . . -0.6 -Parameter space decomposition

  thenLet us now state our n-D stability test.if D(1, . . . , z k , . . . , 1) is not stable then

	3:
	4:	return False
	5:	end if
	6:	end for
	7:
	7:	return False
	8:	end if
	9:	end for
	10:	end for
	11:	return True
	12: end procedure
	Algorithm 2 1: procedure IsStable(D(z 1 , . . . , z n ))	return true if D
		satisfies (2)
	2:	for k = 1 to n do

Table 1

 1 CPU times in seconds of IsStableCAD and IsStableCRIT run on random polynomials in 2 and 3 variables with rational coefficients. t/o: time out

				Running time
	nbvar	degree	density	IsStableCAD IsStableCRIT
		3	sparse dense	0.039 0.047	0.074 0.078
	2	5 8	sparse dense sparse dense	0.055 0.17 0.29 2.74	0.087 0.13 0.21 0.61
		10	sparse dense	1.91 8.59	0.38 1.82
		3	sparse dense	t/o t/o	0.31 0.36
	3	5 8	sparse dense sparse dense	t/o t/o t/o t/o	0.51 1.05 2.31 9.77
		10	sparse dense	t/o t/o	4.71 38.70

The user can choose one of these two methods by means of an option in the routine IsStable.

The experiments were conducted on 1.90 GHz

3-Core Intel i3-3227U with 3MB of L3 cache under a Linux platform.

For particular values of u , gcd(D(z, u ), N (z, u )) can become non-trivial. In this case, to test structural stability, we have to consider D(z, u )/gcd(D(z, u ), N (z, u )) instead of D(z, u ).

This work was supported by the ANR-13-BS03-0005

variety. Of course, doubling the number of variables substantially increases the running time of this step.

Remark 11 Another interesting problem related to the stability of parameterdepending systems concerns the robustness analysis. One typical question is the following: given a nominal value of the parameters that ensures the stability, how to compute the minimal distance from this nominal value that can be tolerated while preserving stability. This is the subject of an ongoing work by the authors where the use of critical methods is investigated in order to compute this minimal distance as a critical value of a distance function restricted to some semi-algebraic set.

An illustrating example

We consider a 2D system defined by a transfer function G (see ( 1)) depending on the parameters u = {u 1 , u 2 } and whose denominator D is defined by:

Applying the algebraic transformation defined in Section 2.2 to D, the bivariate polynomial system ( 17) is defined by:

The minimal discriminant variety of this bivariate system with respect to the projection onto (u 1 , u 2 ) can be obtained by means of the Maple function RootFinding [Parametric][Discriminant Variety]. The discriminant variety is the union of 8 lines, 2 quadrics and 1 curve of degree 6 defined by:

Now, computing the conditions on the parameters u that discriminate the situations where D(z 1 , 1) (resp., D(1, z 2 )) has (or not) roots in the complex unit disc lead to following 5 lines

where [7 u 1 -3 u 2 + 7] is already in the discriminant variety.