Efficient Numerical Analysis of Stability
of High-Order Systems With a Time Delay
George Armanious* and Rick Lind†
University of Florida, Gainesville, FL 32612

arXiv:1710.10289v1 [] 27 Oct 2017

Time delays are a common perturbation in systems with many states, such as networked, distributed,
or decentralized systems. Current methods analyzing the stability of large systems with time delay
typically produce very conservative results. While more exact methods exist, these become inefficient for large systems. This paper provides a methodology for analyzing the stability of time-delayed
systems that is derived from exact methods but is efficient for high-order systems. The computational and memory cost of this new technique is compared to the costs of existing techniques, and its
efficiency is shown using a distributed system with over four hundred states

Nomenclature
A0

=

Undelayed state matrix

A1

=

Delayed state matrix

x

=

System state vector

t

=

Time

τ

=

Time delay

ω

=

Natural frequency at root crossing

T

=

Rekasius substitution variable

n

=

Number of system states

s

=

Frequency domain variable

In

=

n × n identity matrix

I. Introduction
A time delay of τ ∈ R+ notes that dynamics of a system at time of t ∈ R+ depend on some state or input at a
previous time of t − τ . Such a delay can often arise due to hardware associated with sensing, actuation, and communication. It can be challenging to properly account for a time delay when designing a control system; consequently, the
time delay can introduce unexpected instabilities.
This paper considers a special case of time-delayed systems known as retarded time-delayed systems [1]. This case
* Graduate Student, Department of Mechanical and Aerospace Engineering, a1george@ufl.edu
† Associate Professor, Department of Mechanical and Aerospace Engineering, ricklind@ufl.edu

1
American Institute of Aeronautics and Astronautics

of system, shown in Eq. 1, relates a state vector of x ∈ Rn to state matrices of A0 , A1 ∈ Rn×n . The classification of
retarded notes the absence of any ẋ(t − τ ) term.

ẋ(t) = A0 x(t) + A1 x (t − τ )

(1)

The values of τ at which the system changes between stable and unstable relate to the poles of the dynmics, or the
roots of the determinant in Eq. 2, being purely imaginary. A finite number of purely-imaginary roots exist but each
root actually corresponds to an infinite number of time delays because of the transcendental term.

det [sIn − A0 − A1 eτ s ] = 0

(2)

All the time delays at which the system changes stability are actually computed using a single τ and pole value [2].
Essentially, the condition for any pole crossing the imaginary axis depends on the frequency of that pole at that crossing
instead of the time delay associated with that crossing. This relationship between all the time delays associated with
changing stability is given in Eq. 3 for a pole of s = ωi and the set of time delays of τ0 , τ1 , τ2 , . . . associated with
changes in stability. Note that the system may have several poles that cross the imaginary axis at different frequencies
so the stability analysis may want to compute all the sets of ωi and associated τ values.

τk = τ0 +

2π
k,
ωi

k = 0, 1, 2, ...

(3)

A symbolic method is formulated to find all sets of frequencies and time delays associated with change in stability [2, 3]. This method generates an algebraic expression for the characteristic polynomial in terms of s and e−τ s and
then utilizes the the Rekasius substitution to replace the transcendental exponential with a fraction of rational polynomials [4]. The Routh’s array is then used to find the conditions at which stability changes. The method generates an
accurate set of conditions; however, the array is exceptionally difficult to formulate and solve as the number of states
in the system increases beyond a few [5].
A numerical method is also formulated by introducing Kronecker multiplication to recast the n × n transcendental
determinant as a 2n2 × 2n2 eigenvalue problem [1, 7]. These eigenvalues indicate all the poles for which the stability
changes and are used to find the associated time delays. The computations use standard tools for eigenvalues; however,
the memory storage and processing time make the approach nearly infeasible for systems with more than 100 states.
Other approaches are formulated for decentralized systems [9, 10]. These approaches construct a Lyapunov function that generates an upper bound on the time delay that can be tolerated before the onset of instability. These
approaches are straightforward to apply to any system and the computations are not overly burdensome; however,
some function parameters are challenging to choose and the upper bound is often overly conservative.
This paper introduces an alternative method that combines elements of the Rekasius-based approach and the

2
American Institute of Aeronautics and Astronautics

Kronecker-based approach. In this alternative method, the Rekasius substitution is again used but it is introduced
earlier in the formulation than in the previous method. The new formulation is then converted from a determinant
computation to an eigenvalue computation involving a companion matrix in a similar process as used by the previous
method. The resulting formulation is solved numerically using a reasonable amount of memory and computations.

II. Existing Methods
A.

Rekasius Substitution

The characteristic equation of the system in Equation 1 is shown in Equation 4.

det [sIn − A0 − A1 eτ s ] = 0

(4)

The transcendental term makes the system infinite dimensional with infinitely many characteristic roots. The determinant can be carried out symbolically to yield the general form shown in Equation 5, where ak (s) are polynomials of
degree n − k.
det [sIn − A0 − A1 eτ s ] =

n
X

ak (s)e−kτ s

(5)

k=0

The foundation for the algorithm is the Rekasius substitution shown in Eq. 6 for a delay of τ ∈ R+ and an associated
T ∈ R and s ∈ C. The substitution allows the transcendental term of a complex exponential to be replaced with
a fraction of complex polynomials [4]. The critical issue for this substitution is its validity is restricted to purelyimaginary values of s such that s = ω for ω ∈ R; consequently, considering s as a pole of the system implies the
substitution is valid for crossings of a pole between the left-half plane and the right-half plane which are associated
with neutral stability.

e−τ s =

1 − Ts
1 + Ts

(6)

The mapping from T to τ is one-to-many for a given ω and is shown in Equation 7.

τ=


2  −1
tan (ωT ) ∓ ℓπ ,
ω

ℓ = 0, 1, 2, ...

(7)

Applying the Rekasius substitution into Equation 5 yields Equation 8.
n
X

k=0

ak (s)



1 − Ts
1 + Ts

k

=0

3
American Institute of Aeronautics and Astronautics

(8)

By multiplying both sides of Equation 8 by (1 + T s)n , the equation can be simplified to Equation 9.
n
X

ak (s)(1 + T s)n−k (1 − T s)k = 0

(9)

k=0

Equation 10 is obtained by grouping the terms in Equation 9 by powers of s, where bk (T ) are polynomials of T .
2n
X

bk (T )sk = 0

(10)

k=0

The polynomial in Equation 10 has the same roots as the polynomial of Equation 5 but does not contain an
exponential in s. Rather, each coefficient of s is parameterized by a polynomial in T . A Routh’s array can be
constructed for the polynomial of Equation 10, and the first column of the array can be used to determine the pairs,
(T, s) = (Tc , ωc ), that are roots of Equation 10, and Equation 5 by extension.
Finally, the time delay at which the dynamics in Eq. 1 become unstable are computed for those values of T at
which some eigenvalues are purely real. This time delay of τ results from the values of s = ω and T that satisfy
Eq. 10 using the equivalence in Eq. 11.

τ=


2  −1
tan (ωT ) ∓ ℓπ ,
ω

ℓ = 0, 1, 2, ...

(11)

Each value of T and corresponding value of ω are associated with an infinite number of time delays periodically
spaced by

2π
ω .

The transcendality in the original dynamics implies there are infinitely many roots in the characteristic

equation, but literature shows these roots are associated with only a finite number of frequencies. The Rekasius
substitution introduces a new variable, T , and shifts the transcendality from the dynamics of the system to the relation
between T and the associated time delays, τ , at frequencies, ω, where the substitution is valid. Literature also shows
that the direction of a pole crossing at a particular time delay is only a property of the associated frequency. In other
words, all the values of τ related to a particular pair of T and ω by Eq. 11 are either stabilizing or destabilizing [2].
Two properties of a Routh’s array can be taken advantage of to reduce the necessary computations in determining
the purely imaginary roots, s = ωc , and corresponding parameter values, Tc , of Equation 10 [3]. The first property
is that if there is a pair of imaginary roots, the only term on the row corresponding to s1 , defined here as R1 (T ), must
be zero. This property converts analyzing the entire first column into a single root finding problem in T . Also, only
purely real values of T are of interest, since the Rekasius substitution is only defined for T ∈ R. Define Tc as a real
root of R1 (T ).
The second helpful property involves the row of corresponding to s2 . This row has exactly two elements. Define

R21 (T ) and R22 (T ) as the first and second elements of this row. Then for Tc , R21 (Tc )s2 + R22 (Tc ) is a factor of

the characteristic equation in Equation 10. Therefore a subset of the roots of the original characteristic equation can
be found by solving Equation 12 for s.

4
American Institute of Aeronautics and Astronautics

R21 (Tc )s2 + R22 (Tc ) = 0

(12)

In order for the roots of Equation 12 to be purely imaginary, s = ωc , the product of R21 (Tc ) and R22 (Tc ) must
be positive. If this condition is met, the corresponding root crossing is given by Equation 13.

ωc =

s

R22 (Tc )
R21 (Tc )

(13)

An alternative method involving sum of squares optimization can be used to avoid constructing a Routh’s array
from Equation 10 in order to find the upper bound of the first stable time delay region [5]. This method searches
through all values of T and s = ω to find the largest value of T and smallest value of ω such that the magnitude of
Equation 10 is minimized. A sum of squares algorithm is used to efficiently find these values of T and ω given the
characteristic equation.

1. Advantages
The main advantage of this method is its exactness. The root crossings obtained from the Routh’s array are the
locations where the original system will cross the imaginary axis. This method also provides the direction that the
poles move across the imaginary axis at these root crossings and the time delays associated with each root. From all
of this information, a complete stability analysis can be performed which may reveal multiple regions of stability. Of
course, if one is only interested in the smallest time delay that causes instability, only the smallest time delay produced
by this method is needed, as this time delay is guaranteed to make the system unstable, assuming the delay-free system
is stable.

2. Disadvantages
The key disadvantage is that this method is symbolic. In the beginning of the formulation, the computation of a
characteristic equation in terms of two symbolic variables, s and e−τ s , is required. Then, the Rekasius substitution is
used to convert this symbolic equation into an equation in terms of s and a new variable, T . Finally, a Routh’s array
is constructed from this symbolic characteristic equation. While only the final terms of the array are necessary to find
the imaginary roots, the whole array must be constructed to obtain these terms. For a small system, the computational
efficiency of this method is comparable to a numerical method, but for a large system (i.e. a closed-loop distributed
system with multiple controllers and observers), computing the characteristic equation can be tedious or infeasible
depending on the strategy used.
The sum of squares approach does not require construction of a Routh’s array. However, it assumes that the
characteristic equation is symbolically computed for use in the sum of squares algorithm. As mentioned, computing
the characteristic equation of a large system analytically is infeasible.

5
American Institute of Aeronautics and Astronautics

B.

Kronecker Multiplication

The next method is a numerical approach to finding the imaginary root crossings that recasts the n × n transcendental
determinant as a 2n2 × 2n2 eigenvalue problem in terms of the system matrices that can be solved numerically [1, 7].
2

To begin, it is necessary to define two matrix operations. First, the operation ξ : C n×n → C n ×1 converts the square
T

n × n matrix, M = mT1 , mT2 , ..., mTn into a n2 × 1 column matrix as shown in Equation 14.


mT1





 T
m 
 2
:=  . 
 . 
 . 


mTn

(ξM )n2 ×1

(14)

The Kronecker product is also used in this formulation and is defined for any two matrices, An×m and Bp×q as shown
in Equation 15.


 a11 B
 .
.
A⊗B =
 .

an1 B


a1m B 
.. 
. 


· · · anm B

···
..
.

(15)

np×mq

The product of three n × n matrices, P1 , P2 , and P3 can be turned into a Kronecker product of dimension n2 × n2
using the ξ operation and property shown in Equation 16 [11].

ξ(P1 P2 P3 ) = (P3T ⊗ P1 )ξP2

(16)

This method begins with the original time-delayed ODE, shown again in Equation 17.

ẋ = A0 x + A1 x (t − τ )

(17)

The solution is assumed to have the form shown in Equation 18, where v is a vector and s is a scalar.

x(t) = est v

(18)

The derivative of the assumed solution is shown in Equation 19.

ẋ(t) = sest v

(19)

Substituting Equations 18 and 19 into Equation 17 yields Equation 20.

sest v = A0 est v + A1 es(t−τ ) v
6
American Institute of Aeronautics and Astronautics

(20)

Rearranging this equation and factoring est and v yields Equation 21.

est (sIn − A0 − A1 e−sτ )v = 0

(21)

If s = ω is a root, so is −s = −ω. The transpose of Equation 21 is taken to get Equation 22.

v ∗ (−sIn − AT0 ) = eτ s v ∗ AT1

(22)

Multiply Equations 21 and 22 to get Equation 23.




(sIn − A0 ) v (−v ∗ ) sIn + AT0 = e−τ s A1 v eτ s v ∗ AT1

(23)

Define V as the outer product of v with itself: V := vv ∗ . Substituting V into Equation 23 yields Equation 24.


(sIn − A0 ) V sIn + AT0 = −A1 V AT1

(24)

Taking the transpose of both sides of Equation 24 yields Equation 25.


(sIn + A0 ) V sIn − AT0 = −A1 V AT1

(25)

Using the ξ operator on both sides of the equation and applying the property in Equation 16 yields Equation 26.
h

sIn − AT0

T

i
i
h
T
⊗ (sIn + A0 ) ξV = − AT1
⊗ A1 ξV

(26)

Rearranging Equation 26 yields Equation 27.

[(sIn − A0 ) ⊗ (sIn + A) + A1 ⊗ A1 ] ξV = 0

(27)

Define λ(s) as the coefficient of ξV in Equation 27 and rearrange the terms in powers of s as shown in Equation 28.

λ(s) = s2 (In ⊗ In ) + s((In ⊗ A0 ) − (A0 ⊗ In )) + ((A1 ⊗ A1 ) − (A0 ⊗ A0 ))

(28)

Equation 27 simplifies to λ(s)ξV = 0. The nontrivial solution of this equation is obtained by setting det(λ(s)) = 0.
Note that λ(s) is a matrix polynomial, as shown in Equation 29.

λ (s) = G0 s2 + G1 s + G2

7
American Institute of Aeronautics and Astronautics

(29)

The coefficients G0 , G1 , and G2 are defined in Equation 30.
G0 := (In ⊗ In ) = In2
(30)

G1 := ((In ⊗ A0 ) − (A0 ⊗ In ))
G2 := ((A1 ⊗ A1 ) − (A0 ⊗ A0 ))

The goal of setting the determinant of λ(s) to zero can be achieved by linearlizing the polynomial [12]. The companion
matrix of Equation 29 is defined in Equation 31. The eigenvalues of this companion matrix are the values of s that
make the determinant of Equation 29 equal to zero.


0n2


C := 
−G−1
0 G2





  0n2
=
−G−1
G
−G2
1
0
In2


In2 

−G1

(31)

It can be easily shown that the determinant of λ(s) is equal to that of sI − C, as shown in Equations 32-35.

det(sI2n2 − C) =

=
=
=







In2 

 0n2
det sI2n2 − 

−G2 −G1


2
2
sI
−I
n

 n
det 

G2 sIn2 + G1

det In2 s2 + G1 s + G2
det(λ(s))

(32)

(33)
(34)
(35)

With this property, the numerically intensive problem of finding the imaginary root crossings of the system in
Equation 17 simplifies to determining the eigenvalues of C, for which there are many existing algorithms. Once all of
the imaginary roots are found, the corresponding initial time delay can be determined. For an imaginary root, s = ω,
the stability matrix, T = ωIn − A0 − A1 e−ωτ is only singular if |e−ωτ | = 0. The values of z := e−ωτ that cause T
to be singular are the generalized eigenvalues of the pair (A0 − ωIn , −A1 ). Let zd be the generalized eigenvalue that
satisfies |zd | = 1. Then the initial time delay corresponding to an imaginary axis crossing at s = ω is the smallest
positive delay, τ , that satisfies Equation 36, where k is an integer and imag(zd ) and real(zd ) represent the imaginary
and real components of zd , respectively.

τ=

1
ω





−imag(zd )
tan−1
+ 2πk
real(zd )

8
American Institute of Aeronautics and Astronautics

(36)

1. Advantages
This approach is very easy to formulate in code. Given the matrices, A0 and A1 , one only need compute the matrix
C given by Equation 31 using the block matrices obtained by Equation 30. In MATLAB, the kron command can be
used to compute the Kronecker products used by Equation 30. The eigenvalues of C can be computed to obtain all of
the imaginary axis crossings, and for each crossing, a smaller generalized eigenvalue problem is used to compute the
corresponding time delay. Both eigenvalue computations can be performed using the eig command in MATLAB.

2. Disadvantages
For large state-space problems (n ≥ 150), the companion matrix, C, from Equation 31 becomes very large. Each
block element of C has dimensions n2 × n2 as a result of the Kronecker product. Therefore, the total dimensions of
C are 2n2 × 2n2 . Storing such a matrix in memory becomes an issue for large n, and the eigenvalue operation, which
has a theoretical time complexity of ∼ O(n3 ), takes exponentially longer to compute. For an example of how much
system memory is required, consider a state-space size of 200. The corresponding companion matrix has dimensions
80, 000 × 80, 000, and therefore contains 6.4 billion elements. If each element requires 8 bytes of memory to store
(the standard size in MATLAB), 51.2 GB of system memory is required to store this companion matrix! The majority
of modern consumer grade computers have much less RAM (typically between 8 GB and 16 GB), and recall that the
starting state-space size is 200 states. Figure 1 shows required system memory to store the companion matrix as a
function of starting state-space size.

Figure 1: System memory needed to store companion matrix in Kronecker multiplication method

Similarly, computing the eigenvalues of the companion matrix becomes computationally infeasible as the matrix
size grows. In algorithm analysis, big-O notation, O(·), represents the worst case time complexity of an algorithm as

9
American Institute of Aeronautics and Astronautics

a function of problem size. For example, the time for an algorithm that is O(n2 ) to finish grows quadratically as a
function of the problem size, n, in the worst case. Literature has shown that the eigenvalue operation is O(n3 ) [13].
Therefore, the time to compute the eigenvalue of a matrix can be modeled as a cubic polynomial, with the parameters
estimated from measured times. In the following discussion, all computations were done in MATLAB using Windows
10 64-bit running on an Intel i7-4790K processor at 4.0 GHz with 16 GB of RAM. The time taken to compute the
eigenvalue of a random square matrix was measured as a function of matrix size for a 1 × 1 matrix to a 7200 × 7200
matrix. Then, knowing that the computation time is cubic with respect to matrix size, this result is extrapolated for
larger matrix sizes. Figure 2 shows the measured computation times along with the cubic curve fitted to this data.
Figure 3 shows the estimated computation times for much larger matrices.

Figure 2: Measured and curve-fitted CPU time necessary for eigenvalue calculations

Figure 4 shows the estimated computation time necessary to compute the eigenvalues of the companion matrix in
the Kronecker multiplication method based on the curve-fit in Figure 3.
While computation time becomes infeasible for large problem sizes, it is clear that the limiting factor is system
memory. Recall that a problem with 200 states requires 51.2 GB to store the companion matrix in system memory.
From Figure 4, computing the eigenvalues of said companion matrix is estimated to take one day. While this is not
ideal, it is still more feasible than the required system memory.
It may be possible to store a large companion matrix as a sparse matrix. Unfortunately, the sparsity guaranteed
by the companion matrix structure is not enough to guarantee that this storage method is more memory efficient
than storing the complete matrix. It is also possible to use iterative numerical methods to compute a subset of the
eigenvalues of the companion matrix, but none guarantee that the eigenvalues closest to the imaginary axis are found
first.

10
American Institute of Aeronautics and Astronautics

Figure 3: Estimated CPU time necessary for large matrix eigenvalue calculations

Figure 4: Computation time needed to compute eigenvalues of companion matrix in Kronecker multiplication method

11
American Institute of Aeronautics and Astronautics

C.

Lyapunov Methods

The Rekasius substitution method, Kronecker multiplication method, and other similar root finding methods were
primarily developed to analyze small systems (n < 20). For larger state-space systems, these methods become computationally infeasible. Research has been done in analyzing the stability of large time-delayed state space systems,
specifically decentralized systems. These methods typically involve the formulation of a Lyapunov function from
system matrices and tuning parameters. This Lyapunov function is then used to develop an upper-bound for tolerable
time delays. An area of interest in research involves formulating the time delay stability problem as a linear matrix
inequality problem.
The main advantage of these methods is that, once derived, it is easy to apply them to any system, and the computations involved are efficient for large systems. However, the upper-bound produced is typically very conservative
compared to the true value of the largest time-delay that the system can tolerate before it first becomes unstable, especially with poor choices for the tunable parameters. Even with ideal tuning, which can be time consuming, there is no
guarantee that the calculated upper-bound is the largest possible.

III. Methodology
A.

Formulation

The Rekasius substitution enables a formulation for stability to be derived that is computationally advantageous. This
derivation initially introduces the substitution from Eq. 6 to the condition for instability from Eq. 2, rather than after
symbolically expanding the determinant, to obtain Eq. 37. The Rekasius substitution only considers purely-complex
values of s so 1 + T s cannot equal 0 and is thus factored out to obtain Eq. 38 and Eq. 39. The formulation in Eq. 40
results by noting T 6= 0 because the time delay is non-zero. That formulation in Eq. 40 is equivalent to the formulation
in Eq. 41 by noting the determinant of a block matrix. Finally, the term of s is separated to derive the expression in
Eq. 42. Note that Eq. 42 is an eigenvalue problem for the companion matrix of the polynomial in Eq. 39 [11, 12].

12
American Institute of Aeronautics and Astronautics

0 =
=
=
=

=

=





1 − Ts
det sIn − A0 + A1
1 + Ts

n
1
det (s(1 + T s)In − (1 + T s)A0 − (1 − T s)A1 )
1 + Ts

det T In s2 + (In − T A0 + T A1 ) s − (A0 + A1 ) = 0




1
1
In − A0 + A1 s − (A0 + A1 )
det In s2 +
T
T


sI
−I
n
n


det 

1
1
− T (A0 + A1 ) sIn + T In − A0 + A1



0n
In



det sI2n − 


1
1
(A
+
A
)
−
I
−
A
+
A
0
1
n
0
1
T
T

(37)
(38)
(39)
(40)

(41)

(42)

The condition for stability in the presence of a time delay is thus related to an eigenvalue computation; however,
only the eigenvalues must be purely imaginary to indicate a change in stability. Some values of T result in a companion
matrix whose eigenvalues all have a non-zero real part so the system is never unstable for these values of T . Only
values of T resulting in a companion matrix having at least one complex-conjugate pair of eigenvalues without any
real part are associated with the system becoming unstable.
Finally, similar to the approach in Section A, the time delay at which the dynamics in Eq. 1 become unstable are
computed for those values of T at which some eigenvalues are purely real. This time delay of τ results from the values
of s = ω and T that satisfy Eq. 42 using the equivalence in Eq. 11.

B.

Algorithm

A search algorithm is utilized to compute the largest value of time delay for which the dynamics in Eq. 1 remain stable.
The approach searches over values of T to find instances of purely-imaginary values of s that satisfy Eq. 42 and then
substitute in Eq. 11 to find the associated values of time delay.
The concept of the search depends on the dynamics with τ = 0, and consequently the corresponding system in
Eq. 37 for T = 0, being stable. If the companion matrix for some larger magnitude of T has poles on the right-hand
side, then the companion matrix for some value between 0 and T must have poles that crossed the imaginary axis and
consequently the dynamics have become unstable for the associated value of τ .
A pair of searches are actually utilized to find the stabilizing range of time delay. The initial search evaluates
the companion matrix at a coarse grid of T to find when the poles of the companion matrix have changed between
the left-hand plane and the right-hand plane. The final search repeats that evaluation at a fine grid between any two
neighboring values of the coarse grid having the poles changing planes.
The set of T values can be limited to a finite set. This set should not include T = 0 because the system is assumed

13
American Institute of Aeronautics and Astronautics

stable at the value; additionally, the set should not include T = ∞ because the companion matrix converges to a
constant as T approaches infinity. The maximum magnitude of T simply needs to be sufficiently greater than either
kAo + A1 k or kA0 − A1 k−1 so the terms involving T in the companion matrix are negligibly small compared to the
other terms.
The largest value of time delay for which the dynamics remain stable is computed by considering those values of
T from the set, along with the associated ω for each T , for which the companion matrix has some purely-imaginary
eigenvalues. Essentially, a set of τ values are computed using Eq. 11 using all those T and ω. The smallest value of τ
from this set is the time delay at which the dynamics in Eq 1 become unstable.

C.

Computational Characteristics

This new algorithm is more memory efficient than the Kronecker multiplication method and other numerical techniques since it solves multiple small problems rather than one large problem. In each iteration of the algorithm, the
companion matrix of this method is only of size 2n compared to the 2n2 needed by other techniques, given an original state-space of size n. Figure 5 shows the amount of system memory required to store the companion matrix
during each iteration. Comparing this figure to Figure 1, it is clear that considerably less memory is required for this
algorithm.

Figure 5: Memory usage needed to find imaginary axis crossings using the new algorithm as a function of problem
size and number of T values searched through

The necessary CPU time is also reduced for larger problems, but it is now also dependent on the number of T
values used in the search. Obviously, reducing the search space can reduce the accuracy of the solution and lead to
missing root crossing values. Increasing the search space will increase the computation time linearly. Furthermore,
more accurate results require an even larger search space. A good strategy may be to perform a coarse resolution

14
American Institute of Aeronautics and Astronautics

search across a large range of T values to find regions of interest (i.e. |Re(λ)| ≤ 10−3 ), followed by a finer resolution
search within these regions of interest.
The parallelizability of this approach is also attractive. The eigenvalues of the companion matrix at any value of T
is computed independently of any other value of T . The initial overhead of sharing A0 and A1 data and final overhead
of sharing the set of τ values is negligible compared to the benefit of having computation time be a linear function of
number of processors.
Using parallel processing to reduce the total computation time of the problem comes at the cost of increasing the
necessary memory usage. Specifically, the memory usage increases by a factor of the number of processors used. In
most cases, this increase is negligible compared to the memory usage saved compared to most methods in literature and
can be compensated for by dividing the work across multiple machines rather than just multiple processors. Figure 6
shows the increased memory consumption per iteration if the algorithm is run in parallel on four processors of the
same computer.

Figure 6: Memory usage needed to find imaginary axis crossings using the new algorithm when the problem is
parallelized across four processors

IV. Examples
A.

Literature example

The new algorithm is tested with an example from literature to ensure proper operation of the algorithm and assess its
accuracy with respect to a problem whose solution is known. The algorithm is performed on the platform described in
Section II-B. The problem is characterized by the A0 and A1 matrices shown in Equation 43 [2].

15
American Institute of Aeronautics and Astronautics





−5.9
7.1
−70.3
−1
13.5
−1











A0 = −3 −1 −2 ; A1 =  2
−1
5 





2
0
6
−2 −1 −4

(43)

The (T, ω) pairs at which a pole pair crosses the imaginary axis, the corresponding time delays, and the effect on
stability are shown in Table 1. The smallest time delay at which the system becomes unstable is 0.1624 s. The system
becomes stable again in the range 0.1859 < τ < 0.222.
Table 1: Results of Literature Example
T (s)
-0.4269
-0.1332
0.0829
0.0953
0.6233

ω (rad/s)
15.5032
0.8407
3.0347
2.9123
2.1109

Stable/Unstable
Unstable
Stable
Unstable
Stable
Unstable

Time Delay (s)
0.2219, 0.6272, ...
7.208, 14.682, ...
0.1624, 2.233, ...
0.1859, 2.343, ...
0.8725, 3.849, ...

The algorithm reproduces the results shown in literature, with the smallest time delay causing instability as τ =
0.162 s. Based on existing results, an observation is made for applying this algorithm towards a full stability analysis
in the presence of time delays. As T increases near a purely imaginary pole, s = ω, the direction of the sign change
of the real part of s correlates to the effect on stability that the corresponding time delays, τ , have, even though the
Rekasius substitution is only valid at ω. A pole pair crossing from the left-half plane to the right-half plane correlates
to all corresponding time delays having a destabilizing effect, and vice versa.

B.

MAV example

The second example is a 428-state model of a closed-loop aircraft with a distributed control and estimation framework.
The details of this model are found in literature and in the addendum [14]. The model contains distributed sensors and
actuators along with distributed processors at nodes throughout the wing. Some estimation and control is computed at
each node in a decentralized fashion; however, all sensor data is also transferred through sequential-bus hardware to
the central node in the fuselage for estimation and control in a centralized fashion. The bus introduces an appreciable
delay in data transfer to this central node so the associated stability is of critical importance. The state dynamics of
the delay-free model are given in Eq. 44 (definitions of the submatrices are found in literature), with elements in bold
representing sources of time-delay [14].

16
American Institute of Aeronautics and Astronautics

 
ẋ(t)
  A

 

  LC
 x̂(t)
 
 ˙
 

 ˙1   1
 x̂ (t)   LS C1
 

 ..  
..
 .  
.
 

=

 ˙ nn   nn
x̂ (t) LS Cnn
 

 
 1
 ẋK (t)   0
 

 .  
..
 ..  
.
 

 

0
ẋnKn (t)


0

0

···

0

B1 CK1

···

A − LC

0

···

0

B1 C K 1

···

L1C
..
.

X1
..
.

···
..
.

0
..
.

BCK
..
.

···
..
.

LnCn

0

···

X nn

0

···

0
..
.

BK C
..
.

···
..
.

0
..
.

AK
..
.

···
..
.

0

0

0

···

· · · BK C



Bnn CKnn   x(t) 




Bnn CKnn 
  x̂(t) 


 1 
  x̂ (t) 
0


  .. 
..



.
 . 




BCK  x̂nn (t)



 1



0
  xK (t) 



..
  ... 
.




xnKn (t)
AK

(44)

A trio of existing methods is initially applied to this example. The number of states is too large for the symbolic
computation used by the existing Rekasius-substitution method since computing the determinant of a 428 × 428 symbolically to obtain the characteristic equation is not feasible. The Kronecker-multiplication method is also infeasible
for the available computer because the memory storage exceeds 1 TB of space. The only method to compute a valid
solution is the Lyapunov method; however, that method generated an upper bound on stability at a time delay of 10−13
seconds that is clearly conservative as shown by simulations.
The new algorithm is used to analyze the stability of the system by searching over T values in Equation 42 ranging
from -1000 s to 1000 s in increments of 0.001 s and excluding T = 0. The analysis, parallelized across 4 processing
cores, is completed in just over three days of CPU time. The algorithms computes 10.5 sec as the largest value of time
delay for which the system remains stable. Past this time delay, the lateral-directional dynamics of the system becomes
unstable at a frequency of 0.02 Hz. Simulations show that indeed the system is stable for a time delay of 10 sec and
unstable for a time delay of 11 sec. To illustrate this instability, Fig. 7 shows the yaw angle at the c.g. of the aircraft in
response to a step command to the wing-tip twist angle in the presence of no time delay, 10 sec time delay, and 11 sec
time delay.

V. Conclusion
A new algorithm is presented for the analysis of large state-space systems with time-delays. The algorithm is
adapted from the Rekasius substitution method and uses numerical techniques employed in the Kronecker product
method. Therefore, it maintains the exactness of these techniques while being computationally feasible for large
systems. An analysis is performed that compares the CPU cost and memory consumption of the Kronecker product
method with that of the new algorithm. The parallel nature of the algorithm is used to further reduce computation time
necessary for the stability analysis. The algorithm is used on a system with 428 states, and the results are compared
to a Lyapunov-based time-delay algorithm designed for large systems. It is shown that the largest stable time-delay

17
American Institute of Aeronautics and Astronautics

0

No delay
10 sec delay
11 sec delay

−0.05

Yaw angle at cg (deg)

−0.1

−0.15

−0.2

−0.25

−0.3

−0.35

−0.4
0

50

100

150

200

250

300

350

400

Time (s)
Figure 7: Yaw angle of aircraft in response to step command in wing tip twist angle with no time delay, 10 sec time
delay, and 11 sec time delay

computed by this new algorithm is much less conservative than the upper-bound obtained by the Lyapunov-based
algorithm.

VI.

Addendum

The purpose of this submission is to provide matrix information and code that can be used to recreate the results
and compare them to the results in the paper. The system is an aircraft, P , with distributed processors along segments
of the wing with local controllers, K i , and local Kalman filters, F i , and a central processor with a central Kalman
filter, F , receiving delayed inputs. Subscripts denote the subset of the matrix or vector corresponding to a node, and
superscripts denote values computed by that node’s processor.


ẋ(t) =
P

y(t) =
Ki



 ẋ(t)

=

Ax(t) +

P10

i=1

Bi uii (t)

(45)

Cx(t)
AK xK (t) + BK ŷ i (t)


ui (t) =

CK xK (t)

18
American Institute of Aeronautics and Astronautics

(46)

Fi



x̂˙ i (t)

y i (t)
F





= Ax̂i (t) + Bui (t) + LiC x̂(t) − x̂i (t) + LiS Ci x(t) − Ci x̂i (t)

˙
x̂(t)

(47)

C14,68 x̂i (t)

=

= Ax̂(t) +

P

10
i=1


Bi uii (t − τ ) + L (Cx(t − τ ) − C x̂(t))

(48)

The matrices A, B, C, AK , BK , CK , L, LS , and LC are stored in the file “matrices.mat”. The script “delayProof.m” constructs the matrices A0 and A1 corresponding to Equation 1 for this example and runs the new algorithm described in the paper. This new algorithm is in the function “delayStability.m” and accepts any two square
matrices, A and B of the same size as input. The search grid can be modified in the function. The output of the function
is a two-column matrix whose first column is a list of T values where the companion matrix contains an eigenvalue
(pair) near the imaginary axis and the second column contains the corresponding eigenvalue whose imaginary part is
positive. The output of delayProof.m is included for comparison purposes. The script “getDelay.m” performs the second half of the algorithm by processing the matrix output of delayStability.m to perform a finer search approximated
as interpolating the imaginary axis crossing points and computing the time-delays associated with each crossing.
To use the code for a general problem:
• Construct A0 and A1 for that problem
• Modify the grid defined in delayStability.m as desired
• Run the following line: “crossing = delayStability(A0, A1 )”
• Run getDelay with crossing
The delayProof script performs the first three steps for the example shown in the paper.

References
[1] Sipahi, R. and Olgac, N., “Stability Robustness of Retarded LTI Systems with Single Delay and Exhaustive
Determination of their Imaginary Spectra,” SIAM Journal on Control and Optimization, Vol. 45, No. 5, pp. 16801696, 2006, DOI:10.1137/050633238
[2] Olgac, N. and Sipahi, R.,“An Exact Method for the Stability Analysis of Time-Delayed Linear TimeInvariant (LTI) Systems,” IEEE Transactions on Automatic Control, Vol. 47, No. 5, pp. 793-797, 2002, DOI:
10.1109/TAC.2002.1000275
[3] Olgac, N. and Sipahi, R. “An Improved Procedure in Detecting the Stability Robustness of Systems with
Uncertain Delay,” IEEE Transactions on Automatic Control, Vol. 51, No. 7, pp. 1164-1165, 2006, DOI:
10.1109/TAC.2006.878774
19
American Institute of Aeronautics and Astronautics

[4] Rekasius, Z.V., “A Stability Test for Systems with Delays,” in Proc. Joint Automatic Control Conf., 1980, Paper
No. TP9-A, DOI:10.1109/JACC.1980.4232120
[5] Ebenbauer, C. and Allgöwer, F., “Stability Analysis for Time-Delay Systems using Rekasius’s Substitution and
Sum of Squares,” in Proceedings of the 45th IEEE Conference on Decision & Control, IEEE, 2006, pp. 53765381, DOI: 10.1109/CDC.2006.377471
[6] Pakzad, M.A., Pakzad, S., and Nekoui, M.A., “Stability Analysis of Time-Delayed Linear Fractional-Order
Systems,” International Journal of Control, Automation and Systems, Vol. 11, No. 3, pp. 519-525, 2013, DOI:
10.1007/s12555-012-0164-4
[7] Louisell, J., “A Matrix Method for Determining the Imaginary Axis Eigenvalues of a Delay System,” IEEE
Transactions on Automatic Control, Vol. 46, No. 12, pp. 2008-2012, 2001, DOI: 10.1109/9.975510
[8] Jarlebring, E., Meerbergen, K., and Michiels, W., “A Krylov Method for the Delay Eigenvalue Problem,” SIAM
Journal on Scientific Computing, Vol. 32, No. 6, pp. 3278-3300, 2010, DOI: 10.1137/10078270X
[9] Qu, J. and Gao, C., “Stability Analysis for the Large-Scale Systems with Time-Delay,” Journal of Systems
Science and Complexity, Vol. 19, No. 4, pp. 558-565, 2006, DOI: 10.1007/s11424-006-0558-6.
[10] Chen, J., “LMI-Based Stability Criteria for Large-Scale Time-Delay Systems,” JSME International Journal
Series C : Mechanical Systems, Machine Elements and Manufacturing, Vol. 49, No. 1, pp. 225-229, 2006,
DOI:10.1299/jsmec.49.225
[11] Bernstein, D., Matrix Mathematics, pp. 439-458, Princeton University Press, Princeton, 2009.
[12] Gohberg, I., Lancaster, P. and Rodman, L., Matrix Polynomials, pp. 11-22, Academic Press, London, 1982.
[13] Trefethen, L.N. and Bau, D., Numerical Linear Algebra, pp. 190-233, Society for Industrial and Applied Mathematics, Philadelphia, 1997.
[14] Armanious, G. and Lind, R. “Fly-by-Feel Control of an Aeroelastic Aircraft using Distributed Multi-Rate
Kalman Filtering,” Journal of Guidance, Control and Dynamics, Vol. 40, No. 9, pp. 2323-2329, 2017,
DOI:10.2514/1.G002799

20
American Institute of Aeronautics and Astronautics

