Modified Schur-Cohn Criterion for Stability of Delayed Systems

A modified Schur-Cohn criterion for time-delay linear time-invariant systems is derived. The classical Schur-Cohn criterion has two main drawbacks; namely, (i) the dimension of the Schur-Cohn matrix generates some round-off errors eventually resulting in a polynomial of s with erroneous coefficients and (ii) imaginary roots are very hard to detect when numerical errors creep in. In contrast to the classical Schur-Cohn criterion an alternative approach is proposed in this paper which is based on the application of triangular matrices over a polynomial ring in a similar way as in the Jury test of stability for discrete systems.The advantages of the proposed approach are that it halves the dimension of the polynomial and it only requires seeking real roots, making this modified criterion comparable to the Rekasius substitution criterion.


Introduction
Stability theory for linear time invariant systems with delays have long been in existence, but it still remains an active area of research, with applications in power system control, biology, human-in-the-loop control, neural networks, control in the context of automotive engines, haptic control, and synthetic biology, [1].
In this work we treat with linear time invariant timedelayed systems (LTI-TDS) with commensurate delays described by the state-space equation [2]: where  ∈ R  stands for the state vector,  0 ,   ∈ R × are given system matrices, and  > 0 is a delay time.For these systems the characteristic quasi-polynomial   () possesses the form of where  is the degree of commensuracy in the dynamics and   ∈ R[] are polynomials of degree at most  − .
It is well known that the analysis of the stability of LTI-TDS lies on the root continuity argument; that is, the location of the poles of   varies continuously with respect to delay.This means that any root crossing from the left half plane to the right half-plane will need to pass through the imaginary axis; as a result the computation of crossing frequencies is crucial when analyzing the stability and also for the plot of the root locus of   as  varies [3][4][5][6][7][8][9].
There are five criteria in the literature to compute delay margins of general time-delayed systems: (i) Schur-Cohn (Hermite matrix formulation) [2,[10][11][12]; (ii) elimination of transcendental terms in the characteristic equation [13]; (iii) matrix pencil, Kronecker sum method [2,11,12,14]; (iv) Kronecker multiplication and elementary transformation [15]; (v) Rekasius substitution, [4,16,17].These criteria are based on the determination of all the imaginary roots of the characteristic equation and require numerical procedures which result in vast range of precisions in finding pure imaginary roots of polynomials.Some criteria show subtle inapplicabilities mainly from a numerical implementation viewpoint [18].In particular methods (i), (iii), and (iv) lead to polynomials of degree 2 2 , while (ii) requires finding the roots of a polynomial of degree 2  .The Kronecker multiplication method generates the results with the smallest number of significant digits among the above five methods.In the literature, the Rekasius substitution approach is the most attractive since the generated polynomial has degree  2 , much smaller than all the other criteria.Another advantage of this method is that only real roots are searched for, instead of complex ones, which is very crucial from a practical view point since complex roots are difficult to compute when numerical errors are dragged and propagated.For a total comparison of these methods, demonstrating their strengths and weakness, the reader is referred to [19].
This paper is focused on improving the Schur-Cohn procedure so as to obtain the advantageous properties of the Rekasius substitution criterion.Specifically, we modify the Schur-Cohn criterion to generate a real polynomial of degree  2 whose real roots determine the pure imaginary roots of the characteristic equation.The original Schur-Cohn procedure computes the determinant of a partitioned matrix  = ( are appropriate  ×  matrices over the ring of polynomials in the unknown  and Λ * denotes the Hermitian (or conjugated complex) of Λ.The classical Schur-Cohn procedure for high dimensions (2 2 > 10) becomes numerically unreliable due to repeated round-off errors in the determinant expansion procedure unless the operation is performed using very large number of significant digits.This particular point alone brings a weakness from the numerical deployment viewpoint, since the operation ultimately yields poor precision in determining the desired imaginary roots: because of the accumulated numerical errors the imaginary roots shape up in the form of  =  + ,  ≪ 1 and  = √ −1, which leads the numerical error in real parts to be at least in the order of 10 magnitudes larger than the computational precision.Therefore, it is problematic to decide whether these roots, which are very close to the imaginary axis, are really imaginary roots or not.As an attempt to overcome the problems of the Schur-Cohn method some authors have proposed a modification based on the determination of the real eigenvalues of a constant matrix, [11].By contrast, our modification stems from the applications of some transformations on  by using triangular matrices over the polynomial ring in a similar way as in the Jury test of stability for discrete systems [20].
The remainder of the paper is structured as follows: in Section 2, the basic idea of the modified stability criterion is explained and a numerical example is provided.Finally, a summary and an outlook are given in Section 3.

Delayed Systems. Schur-Cohn Modified Method
The stability of the system (1) is determined by its characteristic quasi-polynomial in  [21]: The characteristic quasi-polynomial contains time delays of commensurate nature with degree ; that is,   depends on  − for  = 0, . . ., , with   () ̸ = 0 and rank(  ) ≤  ≤  = rank( 0 ).The starting-point of all the stability criteria is the determination of the roots of   on the imaginary axis.In fact, it is well known from the stability theory of LTI-TDS that the pure imaginary characteristic roots are the only possible transition points from stable to unstable behavior, and vice versa.We want to determine all the imaginary roots of   (;  − , . . .,  − ), and to this end it is convenient to introduce the variable  =  − and write   as a polynomial in two unknowns over the field of reals, We define the bivariate polynomial With every bivariate polynomial over the complex field, (, ), we associate two triangular matrices over the polynomial ring C[]: ) . (5) For the sake of simplicity we write   ,   instead of   (()),   (()).The Schur-Cohn criterion solves the problem of computing the values of  such that   (; , . . .,   ) = 0 by multiplying   (; , . . .,   ) by  − for  = 0, . . ., −1, and   (; , . . .,   ) by   for  = 1, . . ., .This leads to a system of 2 homogeneous equations in the unknowns   ,  = −, . . ., 0, . . ., −1 which has the following matrix representation: where  is the counteridentity matrix (i.e. an antidiagonal matrix with ones), e 1 =   −1 (), e 2 =   −1 ( −1 ) −1 , and . .,  −1 ); the homogeneous system in (7) is the starting point of the Schur-Cohn method.If we make the change of unknowns e   → e  and we take into account the idempotence of , (6) can be transformed into It is easily verified that  *    is of the same structure as   , and in particular it is symmetric.Therefore  *    = ( *    )  =  *    .As a consequence we obtain The homogeneous system (7) has a nontrivial solution so long as det ( and taking determinants in (8) yields As a result, Therefore, the problem has been reduced to find the pure imaginary roots of det( *    −  *    ) which is just the determinant of the Jury matrix   (()) over C[] defined as This can be seen by a similar argument as above in (8), and det(  (())) = det( *    −  *    ).
In the sequel we define the polynomial () := det(  (())) and the set of zeros of a polynomial  is denoted as Z().
From the above development we obtain the following criterion.
Criterion 1.Let  0 ∈ R + be a real root of (√), then ±√ 0 is a pair of pure imaginary roots of the LTI-TDS system in (1).
Proof.It is immediate that ± is a pair of pure imaginary roots of   (; , . . .,   ) if and only if  is a positive real root of   (; , . . .,   ).According to the Schur-Cohn method the roots of   (; , . . .,   ) can be determined by searching for the real roots in the system (7).Additionally, it was shown above that the solution of (7) implies the computing of the roots of the even polynomial (); henceforth, if  0 ∈ R + is a root of (√), ±√ 0 is a pair of real roots of () and ± 0 is a pair of pure imaginary roots of   (; , . . .,   ).
Note that this criterion is advantageous since firstly we only need to compute real roots and secondly deg((√⋅)) = .In the following the set of crossing frequencies is defined as Example 2. Let us consider the LTI-TDS system ẋ () =  0 () +  1 ( − ) where the system matrices are defined as ) . ( The quasi-polynomial characteristic is In the frequency domain ( = ), we define the bivariate polynomial (, ) = ∑ 3 =0   ()  , and from (19) we build the triangular matrices  3 () and  3 () as follows: ) . ( Then we compute the polynomial (√) according to Criterion 1: is the set of pure imaginary roots of the delayed system.Note that from the fundamental theorem of algebra there are only a finite number of crossing frequencies and that Z((√⋅)) ∩ R + represents the set of squares of crossing frequencies.
The method presented above is the first step to test the stability of the characteristic quasi-polynomial   (,  − ) and to compute the delay margin for LTI delay systems.It is common practice to make the following assumption on   .Assumption 3. The characteristic quasi-polynomial   is stable at  = 0, that is, Z(  (⋅; 0)) ⊂ C − , where C − := { ∈ C : Re{} < 0}.
Recall that the delay margin  is the smallest deviation of  from  = 0 such that the system becomes unstable: For any  ∈ [0, ), the system is stable, and whenever  = ∞, the system is said to be stable independent of delay.For each finite delay  we can associate a crossing frequency which represents the first contact or crossing of the roots of   from the stable region to the unstable one.The determination of the delay margin requires solving the bivariate polynomial   ∈ C[, ] in two phases: (i) computing the positive real roots of (√) so as to obtain the crossing frequencies  ∈ Ω of   , (ii) finding the roots of   () ∈ C[] on the boundary of the unit disk D (denoted as D), at each crossing frequency .A detailed analysis of the delay margin  can be found in [11] and in Theorem 2.11, pp.59-60 in [2].Under the assumption that   is stable at  = 0, the two-step stability method is summarized in Table 1.( Under Assumption 3 the system is unstable at  =  and at least a root crosses C from C − to C + at .If we choose  ∈ [0, ), then   ̸ =   and   (  ;  − ) ̸ = 0 for all  ∈ R + .Therefore,   (;  − ) is stable for all  ∈ [0, ).For each crossing frequency   a delay   is computed as shown in Table 2.The delay margin is then computed as the minimum value of these delays, that is,  = min{  :  = 1, . . ., 5} = 0.16230.
The crossing direction of those roots onto the imaginary axis (either from stable to unstable complex plane or vice versa) as  increases is a quantity called the root tendency (RT), [7,8].For a cross frequency  0 and a delay  0 , the root tendency is defined as where sgn(⋅) denotes the sign function.In the following lemma we propose a simple method to compute the root tendency based on the implicit function theorem.
Owing to   ( 0 ;  0 ) = 0 and from the chain rule of derivation, the result in (26) follows immediately.
Example 6.Following with Example 2 and using the root tendency lemma we can compute the crossing direction for each crossing frequency as shown in the fifth column in Table 3; RT = −1 indicates that the crossing is from C + to C − and RT = +1 implies a crossing from C − to C + .
We finish this analysis by noting that it is not possible in general to separate (√) with a factor containing exactly all the real roots of (√) as stated in the following lemma.Lemma 7. Given a polynomial  ∈ R  [] there does not exist a general procedure of factorization of  in polynomials   ∈ R[] and   ∈ R[], () =   ()  (), such that   includes exactly all the real roots of  (and   all the pairs of complex conjugated roots of  not in R).

Example 4 .
Taking into account the positive crossing frequencies   ∈ Ω computed for the LTI-TDS system in Example 2, we can derive the set of roots of   (  ) ∈ C[] onto the unit circle D.Let   ∈ D ∩ Z(  (  )), from the identity   =  −   =   arg(  ) it follows that   = (2 − arg(  ))/  (with arg(  ) ∈ [0, 2]) and the system can cross the stability boundary C = R from the stable halfplane C − to the unstable half-plane C + := C \ (C − ∪ C) or vice versa at   .If the root crosses C from C + to C − at   , there exists a  <   such that makes the root to cross from C − to C + .As a consequence the stability margin is defined as  = min 1≤≤ { 2 − arg (  )   :   ∈ Ω,   ∈ D ∩ Z (  (  ))} .

Table 1 :
Stability analysis of the quasipolynomial   in terms of the delay margin .Set of squares of crossing frequenciesSolutions of   () ∈ C[] on D,  ∈ Ω Stability of the quasi-polynomial

Table 2 :
Delays corresponding to the crossing frequencies of the LTI-TDS system of the Example 2.   (  ) on D (  =  −   )

Table 3 :
Root tendency for the crossing frequencies of the LTI-TDS system of the Example 2.   ( 0 )) − 0  0 ∑