# Modified Schur-Cohn Criterion for Stability of Delayed Systems

Authors: Juan Ignacio Mulero-Martínez
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

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 ofs 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.


---

## Body

## 1. 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 time-delayed systems (LTI-TDS) with commensurate delays described by the state-space equation [2]:(1)x˙t=A0xt+∑k=1m‍Akxt-kτ,m∈Z+,wherex∈Rn stands for the state vector, A0,Ak∈Rn×n are given system matrices, and τ>0 is a delay time. For these systems the characteristic quasi-polynomial pcs possesses the form of(2)pcs,τ=∑k=0q‍akse-kτs,whereq is the degree of commensuracy in the dynamics and ak∈Rs are polynomials of degree at most n-k.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 ofpc 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 pc as τ varies [3–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–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 2q2, while (ii) requires finding the roots of a polynomial of degree q2q. 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 q2, 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 degreeq2 whose real roots determine the pure imaginary roots of the characteristic equation. The original Schur-Cohn procedure computes the determinant of a partitioned matrix A=Λ1sΛ2sΛ2*sΛ1*s where Λ1,Λ2∈Cs are appropriate q×q matrices over the ring of polynomials in the unknown s and Λ* denotes the Hermitian (or conjugated complex) of Λ. The classical Schur-Cohn procedure for high dimensions (2q2>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 s=jω+ϵ, ϵ≪1 and j=-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 A 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 Section2, 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.


## 2. Delayed Systems. Schur-Cohn Modified Method

The stability of the system (1) is determined by its characteristic quasi-polynomial in s [21]:(3)pcs;e-sτ,…,e-msτ=det⁡sI-A0-∑k=1m‍Ake-skτ=∑k=0q‍akse-kτs,(4)a0s=sn+∑i=0n-1‍a0isi,aks=∑i=0n-1‍akisi,k=1,…,q.The characteristic quasi-polynomial contains time delays of commensurate nature with degreeq; that is, pc depends on e-kτs for k=0,…,q, with aqs≠0 and rankAk≤q≤n=rankA0. The starting-point of all the stability criteria is the determination of the roots of pc 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 pcs;e-sτ,…,e-msτ, and to this end it is convenient to introduce the variable z=e-sτ and write pc as a polynomial in two unknowns over the field of reals, pcs,z=∑i=0q‍akszk∈Rs,z.We define the bivariate polynomialr∈Cω,z as rω,z=∑i=0q‍biωzi such that biω=aijω∈Cω. With every bivariate polynomial over the complex field, rω,z, we associate two triangular matrices over the polynomial ring Cω:(5)Uqrω=b¯qωb¯q-1ωb¯q-2ω⋯b¯1ω0b¯qωb¯q-1ω⋯b¯2ω⋮⋮0⋯0b¯qωb¯q-1ω0⋯00b¯qω,Lqrω=0⋯00b0ω0⋯0b0ωb1ω⋮0b0ωb1ω⋯bq-2ωb0ωb1ωb2ω⋯bq-1ω.For the sake of simplicity we writeUq, Lq instead of Uqrω, Lqrω. The Schur-Cohn criterion solves the problem of computing the values of ω such that pcjω;z,…,zq=0 by multiplying pcjω;z,…,zq by z-k for k=0,…,q-1, and pcjω;z,…,zq¯ by zk for k=1,…,q. This leads to a system of 2q homogeneous equations in the unknowns zk, k=-q,…,0,…,q-1 which has the following matrix representation:(6)LqU¯qJJUqJL¯qJe1e2=00,whereJ is the counteridentity matrix (i.e. an antidiagonal matrix with ones), e1=Jlq-1Tz, e2=Jlq-1Tz-1z-1, and lq-1Tz=1,z,z2,…,zq-1; the homogeneous system in (7) is the starting point of the Schur-Cohn method. If we make the change of unknowns ei↦Jei and we take into account the idempotence of J, (6) can be transformed into(7)LqU¯qJUqJL¯qJe1Je2=00.It is easily verified thatUq*L¯q is of the same structure as Lq, and in particular it is symmetric. Therefore Uq*L¯q=Uq*L¯qT=Lq*U¯q. As a consequence we obtain(8)Lq*-Uq*OIIOOJLqU¯qJUqJL¯q=Lq*-Uq*OILqU¯qUqL¯q=Lq*Lq-Uq*UqOUqL¯q.The homogeneous system (7) has a nontrivial solution so long as det⁡LqU¯qJUqJL¯q=0. We use the fact that(9)det⁡Lq*-Uq*OIIOOJ=det⁡Lq*-Uq*JOJ=det⁡Lq*J=det⁡L¯q-1q,and taking determinants in (8) yields(10)det⁡L¯q-1qdet⁡LqU¯qJUqJL¯q=det⁡L¯qdet⁡Lq*Lq-Uq*Uq.As a result,(11)det⁡LqU¯qJUqJL¯q=-1qdet⁡Lq*Lq-Uq*Uq.Therefore, the problem has been reduced to find the pure imaginary roots ofdet⁡Lq*Lq-Uq*Uq which is just the determinant of the Jury matrix Jqrω over Cω defined as(12)Jqrω=LqUqU¯qL¯q∈C2q×2qω.This can be seen by a similar argument as above in (8),(13)Lq*-Uq*OILqU¯qUqL¯q=Lq*Lq-Uq*UqLq*U¯q-Uq*L¯qUqL¯q=Lq*Lq-Uq*UqOUqL¯q,anddet⁡Jqrω=det⁡Lq*Lq-Uq*Uq.In the sequel we define the polynomialζω≔det⁡Jqrω and the set of zeros of a polynomial p is denoted as Zp.Proposition 1.
The real polynomialζω is even.Proof.
We prove thatζλ=0 if and only if ζ-λ=0. First of all note that bi-ω=biω¯ and, then, r-ω,z=r¯ω,z. As a consequence,(14)Uqr¯ω=U¯qrω=Uqr-ω,Lqr¯ω=L¯qrω=Lqr-ω.
Thus,(15)Lq*r-ωLqr-ω-Uq*r-ωUqr-ω=LqTrωL¯qrω-UqTrωU¯qrω,
andζλ=ζ-λ.From the above development we obtain the following criterion.Criterion 1.
Letω0∈R+ be a real root of ζs, then ±jω0 is a pair of pure imaginary roots of the LTI-TDS system in (1).Proof.
It is immediate that±jω is a pair of pure imaginary roots of pcs;z,…,zq if and only if ω is a positive real root of pcjω;z,…,zq. According to the Schur-Cohn method the roots of pcjω;z,…,zq 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 ζs, ±ω0 is a pair of real roots of ζs and ±jω0 is a pair of pure imaginary roots of pcs;z,…,zq.Note that this criterion is advantageous since firstly we only need to compute real roots and secondlydeg⁡ζ·=q. In the following the set of crossing frequencies is defined as(16)Ω≔ω∈R+:ω2∈Zζ·∩R+.Example 2.
Let us consider the LTI-TDS systemx˙t=A0xt+A1xt-τ where the system matrices are defined as(17)A0=-113.5-1-3-1-2-2-1-4,A1=-5.97.1-70.32-15206.
The quasi-polynomial characteristic is(18)pcs;e-sτ,…,e-msτ=a3se-3τs+a2se-2τs+a1se-τs+a0s
with(19)a0s=s3+6s2+45.5s+111a1s=0.9s2-116.8s-22.1a2s=90.9s-185.1a3s=119.4.
In the frequency domain (s=jω), we define the bivariate polynomial rω,z=∑i=03‍aijωzi, and from (19) we build the triangular matrices U3r and L3r as follows:(20)U3=119.4-185.1-0.9ω2-22.10119.4-185.100119.4+jω0-90.9116.800-90.9000L3=00111.0-6.0ω20111.0-6.0ω2-0.9ω2-22.1111.0-6.0ω2-0.9ω2-22.1-185.1+jω0045.5-ω2045.5-ω2-116.845.5-ω2-116.890.9.
Then we compute the polynomialζs according to Criterion 1:(21)ζs=s9-165.81s8-18682s7+2.4869×105s6+35185s5-4.1197×109s4+9.4110×1010s3-7.1137×1011s2+1.8933×1012s-1.0144×1012,
whose positive real roots areZζ·∩R+=8.4802,4.4564,9.2141,0.70634,240.35 and then(22)±jΩ=±j2.9121,±j2.111,±j3.0355,±j0.84044,±j15.503
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 thatZζ·∩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-polynomialpcs,e-sτ and to compute the delay margin for LTI delay systems. It is common practice to make the following assumption on pc.Assumption 3.
The characteristic quasi-polynomialpc is stable at τ=0, that is, Zpc·;0⊂C-, where C-≔s∈C:Res<0.Recall that the delay marginτ¯ is the smallest deviation of τ from τ=0 such that the system becomes unstable:(23)τ¯≔min⁡τ≥0:pcjω,e-jωτ=0forsomeω∈R+.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 pc from the stable region to the unstable one. The determination of the delay margin requires solving the bivariate polynomial pc∈Cω,z in two phases: (i) computing the positive real roots of ζs so as to obtain the crossing frequencies ω∈Ω of pc, (ii) finding the roots of pcjω∈Cz 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 pc is stable at τ=0, the two-step stability method is summarized in Table 1.Table 1
Stability analysis of the quasipolynomialpc in terms of the delay margin τ-.
Set of squares of crossing frequencies
Solutions ofpc(ω)∈C[z] on ∂D, ω∈Ω
Stability of the quasi-polynomialpc∈C[ω,z]
Z
(
ζ
(
·
)
)
∩
R
+
=
⌀
Stable, independent of delay (τ-=∞)
Z
(
ζ
(
·
)
)
∩
R
+
=
(
0
)
Z
(
ζ
(
·
)
)
∩
R
+
≠
⌀, Z(q)∩R+≠(0)
Z
(
p
c
(
j
ω
)
)
∩
∂
D
=
⌀
Z
(
p
c
(
j
ω
)
)
∩
∂
D
≠
⌀
Stable, dependent of delay (τ-<∞)Example 4.
Taking into account the positive crossing frequenciesωk∈Ω computed for the LTI-TDS system in Example 2, we can derive the set of roots of pcjωk∈Cz onto the unit circle ∂D. Let zk∈∂D∩Zpcjωk, from the identity zk=e-jωkτ=ejargzk it follows that τk=2π-argzk/ωk (with argzk∈0,2π) and the system can cross the stability boundary ∂C=jR from the stable half-plane C- to the unstable half-plane C+≔C∖C-∪∂C or vice versa at τk. If the root crosses ∂C from C+ to C- at τk, there exists a τ<τk such that makes the root to cross from C- to C+. As a consequence the stability margin is defined as(24)τ-=min⁡1≤k≤q2π-argzkωk:ωk∈Ω,zk∈∂D∩Zpcjωk.Under Assumption3 the system is unstable at τ=τ- and at least a root crosses ∂C from C- to C+ at τ-. If we choose τ∈0,τ-, then τωk≠θk and pcjωk;e-jωτ≠0 for all ω∈R+. Therefore, pcs;e-sτ is stable for all τ∈0,τ-. For each crossing frequency ωk a delay τk is computed as shown in Table 2. The delay margin is then computed as the minimum value of these delays, that is, τ-=min⁡τk:k=1,…,5=0.16230.Table 2
Delays corresponding to the crossing frequencies of the LTI-TDS system of the Example2.
Crossing frequency (ωk)
Roots ofpc(jωk) on ∂D (zk=e-jωkτ)
Phase ofzk(θk)
Delay (τk)
2.9121
0.85692
-
j
0.51548
5.7416
0.18597
2.111
-
0.26768
-
j
0.96352
4.4414
0.87247
3.0355
0.88110
-
j
0.47297
5.7905
0.16230
0.84044
0.97521
+
j
0.22129
0.22314
7.2106
15.503
-
0.95537
+
j
0.29541
2.8417
0.22199The 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 (25)RTω,τ=sgnRe∂s∂τjω0,τ0,wheresgn· denotes the sign function. In the following lemma we propose a simple method to compute the root tendency based on the implicit function theorem.Lemma 5 (root tendency).
Letpc be a characteristic quasi-polynomial as defined in (3) and let τ0 be the delay for which the system crosses the stability boundary at jω0, for some ω0∈Ω, that is, pcjω0;τ0=0. Let us assume that ∂pcs;τ/∂sjω0,τ0=∑k=0q‍ak′jω0-k-1τakω0e-jω0kτ0≠0. Then, the root tendency of jω0 is given by(26)RTω0,τ0=-sgne-jkω0τ0∑k=1qk-1-1Im∑k=1qkakjω0e-jkω0τ0HHHHHHHHHH×∑k=0qe-jkω0τ0ak′jω0-k-1τ0akjω0HHHHHHHHHHHH×e-jkω0τ0∑k=1qk-1-1.Proof.
In virtue of the implicit function theorem, there exist aδ>0, a neighborhood U0⊂C of jω0, and a continuous function s:τ0-δ,τ0+δ→U0 such that pcsτ;τ=0 for all τ∈τ0-δ,τ0+δ; this is also obvious from the continuity of the roots s· with respect to the delay τ. Furthermore, since pc is differentiable and ∂pcs;τ/∂sjω0,τ0≠0, it is clear that(27)∂s∂τjω0,τ0=-∂pcs;τ/∂τjω0,τ0∂pcs;τ/∂sjω0,τ0,
where(28)∂pcs;τ∂sjω0,τ0=∑k=0q‍ak′jω0-k-1τ0akjω0e-jkω0τ0s+pcjω0;τ0,∂pcs;τ∂τjω0,τ0=-jω0∑k=1q‍kakjω0e-jkω0τ0.
Owing topcjω0;τ0=0 and from the chain rule of derivation, the result in (26) follows immediately.Example 6.
Following with Example2 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+.Table 3
Root tendency for the crossing frequencies of the LTI-TDS system of the Example2.
Crossing frequency (ωk)
Delay (τk)
∑
k
=
0
q
‍
(
a
k
′
(
j
ω
0
)
-
(
k
-
1
)
τ
0
a
k
(
j
ω
0
)
)
e
-
j
ω
0
k
τ
0
∑
k
=
1
q
‍
k
a
k
(
j
ω
0
)
e
-
j
k
ω
0
τ
0
(
R
T
)
ω
,
τ
2.9121
0.18597
-
48.520
+
j
30.258
74.136
-
j
58.701
-
1
2.111
0.87247
-
140.62
+
j
342.86
148.78
-
j
182.37
+
1
3.0355
0.16230
-
47.511
+
j
27.405
94.101
-
j
41.172
+
1
0.84044
7.2106
872.31
-
j
692.90
-
119.43
+
j
99.556
-
1
15.503
0.22199
-
902.35
-
j
919.93
1825.5
+
j
4475.4
+
1We finish this analysis by noting that it is not possible in general to separateζs with a factor containing exactly all the real roots of ζs as stated in the following lemma.Lemma 7.
Given a polynomialp∈Rns there does not exist a general procedure of factorization of p in polynomials pR∈Rs and pC∈Rs, ps=pRspCs, such that pR includes exactly all the real roots of p (and pC all the pairs of complex conjugated roots of p not in R).Proof.
We prove this assertion by a counterexample: without loss of generality we assume thatp is monic. Let us assume that p has only real roots which are ordered as α1<α2<⋯<αm, m≤n. Always we can find a value δ∈R such that α1<δ<α2. We make the change of variable s↦s+δ resulting into a polynomial pδs=ps+δ. The roots of pδ are located at αk′=αk-δ and it is easy to see that α1′<0 and αk′>0 for k=2,…,m. Then we apply the map s↦s2 which produces a polynomial qs=pδs2 with a pair of pure conjugated complex roots at ±iδ-a1 and pairs of symmetric real roots at ±ak-δ for k=2,…,m. Under the hypothesis of separability, qs=qRsqCs and, undoing in reverse order the changes of variables, s↦s↦s+δ, we obtain that qRs-δ=s-α1 and qCs-δ=ps/s-α1. Therefore, we might isolate roots of p by successively repeating the procedure onto qC. This is a contradiction since from the Galois theory it is well known that it is not always possible to solve a polynomial equation by radicals.


## 3. Conclusions

We overcome the drawbacks of the original Schur-Cohn criterion for delay systems, namely, the high degree of the polynomial that emerges from the Schur-Cohn matrix and the difficulty of detection of complex roots, especially pure imaginary roots, due to the propagation of numerical errors. In particular, we decrease the degree of the generated polynomial to a half and we reduce the problem numerically to that of seeking the real roots of this polynomial, not the complex ones. These properties are exactly those that make the Rekasius substitution criterion so attractive in the stability theory of delay systems. Therefore, from a practical viewpoint, the modified Schur-Cohn criterion is comparable to other widely used criteria in the literature. Our method is based on making transformations on the Schur-Cohn matrix via triangular matrices over a polynomial ring in one unknown. This leads to the determinant of the Jury matrix over a polynomial ring, which is proved to be an even polynomial. The pure imaginary roots of the Schur-Cohn polynomial are those of the even polynomial, and in virtue of the change of variableω↦ω, we only need to seek real roots, which is more precise from a numerical viewpoint, while the degree of the polynomial is halved.



---

*Source: 846124-2015-03-19.xml*