INT J COMPUT COMMUN, ISSN 1841-9836
8(6):863-868, December, 2013.

Direct Method for Stability Analysis of Fractional Delay Systems

M.A. Pakzad, M.A. Nekoui

Mohammad Ali Pakzad*
Department of Electrical Engineering
Science and Research Branch, Islamic Azad University
Tehran, Iran
*Corresponding author: m.pakzad@srbiau.ac.ir

Mohammad Ali Nekoui
Faculty of Electrical and Computer Engineering
K.N.Toosi University of Technology, Seyed-Khandan
P.O. Box 16315-1355, Tehran, Iran
manekoui@eetd.kntu.ac.ir

Abstract: In this paper, a direct method is presented to analyze the stability of
fractional order systems with single and multiple commensurate time delays. Using
the approach presented in this study, first, without using any approximation, the tran-
scendental characteristic equation is converted to an algebraic one with some specific
crossing points. Then, an expression in terms of system parameters and imaginary
root of the characteristic equation is derived for computing the delay margin. Finally,
the concept of stability is expressed as a function of delay. An illustrative example is
presented to confirm the proposed method results.
Keywords: Control, fractional delay systems, stability windows, Root-Locus.

1 Introduction

Time delay is an inherent part of many dynamical and physical systems. The device delay
often appears in computer based control systems, the wireless and web-based control systems,
communication systems and so on. Also, the existence of various sensors and actuators in the
feedback loop is the cause of delay in many systems. Moreover, taking an automatic computer-
based control system or a process control with networked transmission for example, it is necessary
and more reasonable to simultaneously consider the possible transmission delay for the corre-
sponding control law and the device delay due to computing or processing. Recently, much
attention has been paid to the subjects of stability, stabilization and control of the time delay
systems [1] - [5]. Fractional order models would be more accurate than integer order models. In
fact, real world processes generally or most likely are fractional order systems.

In a generic sense, the fractional-order systems are identified by non-integer powers of the
Laplace variable s. When time delays and fractional-order derivative are involved in dynamical
systems, we have fractional-delay systems. The stability of fractional-delay systems can be
determined by the root location of the characteristic equations. The most natural way to find the
location of the roots of a linear fractional order system is to solve its corresponding characteristic
equation. But in our case this will be a transcendental one, being thus generally impossible to
solve it directly. For this reason, most of the existing approaches study stability of such systems
by finding the crossings of poles through the imaginary axis [8]. There has been a large effort
to deal with this problem, as can be seen by the large quantity of articles dealing with it for the
standard case (integer order systems); see [5], and others.

In [6], the necessary and sufficient conditions for the BIBO stability of fractional order delay
systems have been introduced. From the numerical analysis point of view, the effective numerical
algorithms have been discussed in [7] and [8] for the evaluation of BIBO stability of fractional

Copyright © 2006-2013 by CCC Publications



864 M.A. Pakzad, M.A. Nekoui

order delay systems. In [9] a heavy computation scheme based on the Cauchy’s integral has been
proposed to test the stability of such systems, and in [10], a technique based on the Lambert W
function was used for the same purpose. In this work, an analitical direct approach to determine
all possible stability regions in the parametric space of delay is proposed. The original idea in
this strategy is derived based on the method reported in [5] to achieve the stability criteria of
integer order delay systems. The proposed method is an analytically elegant procedure that
first converts the transcendental characteristic equation into an algebraic equation without the
transcendentality by eliminating the highest degree of commensuracy terms successively. The
resulting algebraic equation without the transcendentality also enables us to easily determine the
delay dependency of the system stability and the sensitivities of crossing roots (root tendency)
with respect to the time delay.

2 Preliminaries and Defintions

A general class of fractional order LTI (linear time invariant) systems with multiple commen-
surate delays of retarded type is taken into account:

D
1/α
t x(t) =

d1/αx(t)

dt1/α
= A0 +

n∑
ℓ=1

Aℓx(t − ℓτ)

The fractional delay characteristic equation of above system can be expressed in the general form
of

C(α
√
s,τ) = det

(
s1/αI − A0 −

n∑
ℓ=1

Aℓe
−ℓτs

)
=

n∑
ℓ=0

Pℓ(
α
√
s) e−ℓτs = 0 (1)

where parameter τ is non-negative, such that τ ∈ R+ and Pℓ(α
√
s) for ℓ ∈ NN is a real polynomial

in the complex variable α
√
s (where α ∈ N). Note that the zeros of characteristic equation (1)

are in fact the poles of the system under investigation. We find out from [6] that the transfer
function of a system with a characteristic equation in the form of (1) will be H∞ stable if, and
only if, it doesn’t have any pole at ℜ(s) ≥ 0. For fractional order systems, if a auxiliary variable
of v = α

√
s is used, a practical test for the evaluation of stability can be obtained. By applying

this auxiliary variable in characteristic equation (1), the following relation is obtained:

Cv(v,τ) =

n∑
ℓ=0

Pℓ(v) e
−ℓτvα (2)

This will transform the domain of the system from a multisheeted Riemann surface into the
complex plane, where the poles can be easier calculated. In this new variable, the instability
region of the original system is not given by the right half-plane, but in fact by the region
described as:

|∠v| ≤
π

2α
(3)

with v ∈ C . Let us assume that ; s = ±jω or in other words, s = ωe±jπ/2 are the roots of
characteristic equation (1) for a τ ∈ R+. Then for the auxiliary variable, the roots are defined
as follows:

v = α
√
s = α

√
ω e±jπ/2α (4)



Direct Method for Stability Analysis of Fractional Delay Systems 865

3 Crossing position

The main objective of this section is to present a new method for the evaluation of stability
and determination of the unstable roots of a fractional delay system. The proposed method elim-
inates the transcendental term of the characteristic equation without using any approximation
and converts it into a equation without the transcendentality such that its real roots coincide
with the imaginary roots of the characteristic equation exactly. If the characteristic equation (1)
has a solution of s = jωc then C(α

√
−s,τ) = 0 will have the same solution.

C(α
√
−s,τ) =

n∑
ℓ=0

Pℓ(
α
√
−s) eℓτs = 0 (5)

Characteristic equation (5) can be written in terms of the auxiliary parameter v as:

C(v,τ) =

n∑
ℓ=0

Pℓ(v) e
ℓτvα = 0 (6)

A recursive procedure should be developed to achieve that purpose. Therefore, let us define

C(1)(v,τ) =

n−1∑
ℓ=0

[P0(v)Pℓ(v) − Pn(v)Pn−ℓ(v)] e−ℓτv
α
=

n−1∑
ℓ=0

P
(1)
ℓ (v) e

−ℓτvα = 0 (7)

Then, we have

C(1)(v,τ) =

n−1∑
ℓ=0

[P0(v)Pℓ(v) − Pn(v)Pn−ℓ(v)] eℓτv
α
=

n−1∑
ℓ=0

P
(1)
ℓ (v) e

ℓτvα = 0 (8)

Where
P

(1)
ℓ (v) = P0(v)Pℓ(v) − Pn(v)Pn−ℓ(v) (9)

We can easily repeat this procedure to eliminate commensuracy terms successively by defining
a new polynomial

P
(r+1)
ℓ (v) = P

(r)
0 (v)P

(r)
0 (v) − P

(r)
n−r(v)P

(r)
n−r−ℓ(v) (10)

and an augmented characteristic equation

C(r)(v,τ) =

n−r∑
ℓ=0

P
(r)
ℓ (v) e

−ℓτvα = 0 (11)

By repeating this procedure n times, we eliminate the highest degree of commensuracy terms
and obtain the following augmented characteristic equation

C(n)(v) = P
(n)
0 (v) = 0 (12)

Where
P

(n)
0 (v) = P

(n−1)
0 (v)P

(n−1)
0 (v) − P

(n−1)
1 (v)P

(n−1)
1 (v) (13)

It should be emphasized that that if s = jωc is the solution of (1) for some τ, then it is also a
solution of (12). If we substitute v = (e−jπ/α)v and v = α

√
ωc e

jπ/2α in (13), we get the following
equation in ω

D(ω) =
(
P

(n−1)
0 (e

−jπ/α v) P
(n−1)
0 (v) − P

(n−1)
1 (v) P

(n−1)
1 (e

−jπ/α v)
)∣∣∣

v=α
√
ωc ejπ/2α

(14)



866 M.A. Pakzad, M.A. Nekoui

The corresponding value of time delay is then computed by

τ∗ =
1

ωc
tan−1

( ℑ[(P(n−1)0 (v))/(P(n−1)1 (v))]
ℜ[−(P(n−1)0 (v))/(P

(n−1)
1 (v))]

)∣∣∣
v=α

√
ωc e

jπ
2α

+
2kπ

ωc
;k ∈ Z+ (15)

Theorem 1. A system with characteristic equation (1) has finite crossing points for any τ ∈ R+.

Proof: Let assume s = jωc be a pair of roots for C(α
√
s,τ) then v = α

√
s = α

√
ωc e

±jπ/2α would
be a pair of roots for D(ω). since D(ω) is a finite degree polynomial with maximum degree of
(P

(n)
0 (v)) then the number of crossing points of (1) that are the real roots of D(ω) is finite. 2

The whole ω values, for which s = jω is a root of equation (1) for some non-negative delays,
is defined as the crossing frequency set.

Ω =
{
ω ∈ R+

∣∣∣ C(α√s,τ) = 0, for some τ ∈ R+} (16)
Corollary 2. If the system given as (1) is stable for τ = 0 (i.e. system without delay) and
Ω = ϕ, then the system will be stable for all positive values of τ ∈ R+.

Proof: From the fact that there are no roots crossing the imaginary axis. 2

4 Direction of crossing

After the crossing points of characteristic equation (1) from the imaginary axis are obtained,
the goal now is to determine whether each of these root crossings from the imaginary axis is a
stabilizing cross or a destabilizing cross. Assume that (s,τ) is a simple root of C(α

√
s,τ) = 0.

The root tendency for each ωcm and τmk is defined as:

Root Tendency = RT |s=jωc = sgn
(
ℜ
(
Ssτ|s=jωcm

τ=τmk

))
= sgn

(
ℜ
(
−

∂C/∂τ
∂C/∂s

∣∣∣∣
s=jωc

))
(17)

If it is positive, then it is a destabilizing crossing, whereas if it is negative, this means a stabilizing
crossing. Notice that root tendency represents the direction of transition of the roots at jωcm as
τ increases from τmk−ϵ to τmk−ϵ , 0 < ϵ ≪ 1.

Theorem 3. The root tendency at a crossing, jωc is invariant with respect to time delay τmk .

Proof: One can find ds/dτ for simple roots of (1) as follows:

ds

dτ
= −

∂C/∂τ
∂C/∂s

=

∑n
ℓ=0 Pℓ(

α
√
s) ℓ s e−ℓτs∑n

ℓ=0
Pℓ(

α
√
s)

ds
e−ℓτs −

∑n
ℓ=0 Pℓ(

α
√
s)ℓτse−ℓτs

(18)

Based on (18) and definition given in (17), the root tendency of each time delay τmk is written
as follows:

RT|τs=jωc = sgn
(
ℜ
(
Ssτ|s=jωc

))
= sgn

(
ℜ
(
−

∂C/∂τ
∂C/∂s

∣∣∣∣
s=jωc

))
= sgn

(
ℜ
( ∑n

ℓ=0 Pℓ(
α
√
s)ℓse−ℓτs∑n

ℓ=0
Pℓ(

α
√
s)

ds
e−ℓτs − Pℓ(α

√
s)ℓe−ℓτs

))
= sgn

(
ℜ
( ∑n

ℓ=0
Pℓ(

α
√
s)

ds
e−ℓτs∑n

ℓ=0 Pℓ(
α
√
s) ℓ s e−ℓτs

−
τ

s

)−1)
= sgn

(
ℜ
( ∑n

ℓ=0
Pℓ(

α
√
s)

ds
e−ℓτs∑n

ℓ=0 Pℓ(
α
√
s)ℓse−ℓτs

))
= sgn

(
ℑ
( ∑n

ℓ=0
Pℓ(

α
√
s)

ds
e−ℓτs∑n

ℓ=0 Pℓ(
α
√
s)ℓe−ℓτs

) )∣∣
s=jωc

τ=τm1+
2kπ
ωc

(19)

The root tendency in each time delay τmk is independent from the time delay itself and constant
for each crossing frequency, because e−ℓτs and Pℓ(α

√
s) do not depend on τmk. 2



Direct Method for Stability Analysis of Fractional Delay Systems 867

5 Illustrative Example

We present an example case, which display all the features discussed in the text.

Example 4. This example has been taken from [3] and [8]. Consider the following linear time-
invariant fractional order system with one delay:

C(
√
s,τ) = (

√
s)3 − 1.5(

√
s)2 + 4(

√
s) + 8 − 1.5(

√
s)2e−τs (20)

This system has a pair of poles (s = ±8j) on the imaginary axis for τ = 0. A very involved
calculation scheme based on Cauchy’s integral has been used in [9] to show that this system is
unstable for τ = 0.99 and stable for τ = 1. Our objective in this example is to find all the stability
windows based on the method described in this article for this system. Using auxiliary variable
v =

√
s into (20) the characteristic equation is obtained as:

Cv(v,τ) = v
3 − 1.5v2 + 4v + 8 − 1.5v2 e−τv

2
(21)

By applying the criterion expressed in the previous section, we can eliminate exponential term
from (21) as follows:

jv6 + 1.5(1 − j)v5 + 14(1 + j)v3 − j16v2 + 32(1 − j)v + 64 = 0 (22)

By inserting expression v =
√
ω ej

π
4 = a(1 + j) in the above equation , we get:

8a6 − 12a5 − 56a3 + 32a2 + 64a + 64 = 0 (23)

Where a =
√
ω/2. The real solutions of (23) for a is

a = 1.82 , v2 = 1.82(1 + j) , ω = 6.6248 (24)

a = 2 , v1 = 2(1 + j) , ω = 8 (25)

The corresponding infinite countable time delays of the cross points in (24) and (25) are obtained
with regards to relation (15) as τ1k = 0.0499 + 0.3019kπ and τ2k = 0.25kπ, respectively.

By applying the criterion expressed in the previous section, it is easy to find out that a
destabilizing crossing of roots (RT = +1) has occurred at τ = 0.25kπ for s = ±j6.6248 and a
stabilizing crossing (RT = −1) has taken place at τ = 0.0499 + 0.3019kπ for s = ±j8 for all
values of k ∈ Z+. Therefore, we will have 5 stability windows as follows: 0.0499 < τ < 0.7854
, 0.9983 < τ < 1.5708 , 1.9486 < τ < 2.3562 , 2.8953 < τ < 3.1416 and 3.8437 < τ < 3.9270
which agree with the results presented by [3] and [8]. Note that at τ = 3.9269 , an unstable pair
of poles crosses toward the right half-plane, and before this unstable pole pair can turn to the
left half-plane at τ = 4.7922 , another unstable pair of poles goes toward the right half-plane at
τ = 4.7123 ; and thus, the system can not recover the stability. To get a better understanding
of the properties of this system, its Root-Locus curve has been plotted as a function of delay in
Fig.1.

6 Conclusions

An efficient method to analyze the BIBO stability of a large class of time-delayed fractional
order systems for both single and commensurate-delay cases is proposed. The method introduces
an augmented equation whose real roots give the finite values of crossing frequencies at which
stability feature of the system change. According to the infinitely countable time delays corre-
sponding to each crossing point, the parametric space of τ is discretized to investigate stability
in each interval. Finally, an illustrative example is presented to highlight the proposed approach.



868 M.A. Pakzad, M.A. Nekoui

Figure 1: Root-loci for C(
√
s,τ) until τ=3.9

Bibliography

[1] Stojanovic, S.B.; Debeljkovic,D.LJ.; Dimitrijevic, N. (2012); Stability of Discrete-Time Sys-
tems with Time-Varying Delay: Delay Decomposition Approach, Int J Comput Commun,
ISSN 1841-9836, 7(4): 775-783.

[2] Liu, C.L.; Liu,F. (2010); Consensus Problem of Second-order Dynamic Agents with Heteroge-
neous Input and Communication Delays, Int J Comput Commun, ISSN 1841-9836, V(3):325-
335.

[3] Pakzad, M. A.; Pakzad, S.; Nekoui, M.A. (2013); Stability analysis of time-delayed linear
fractional-order systems, International Journal of Control, Automation, and Systems, 11(3):
519-525.

[4] Pakzad, S.; Pakzad, M. A. (2011), Stability condition for discrete systems with multiple state
delays, WSEAS Trans. on Systems and Control, 6(11): 417-426.

[5] Walton, J KE.; Marshal, JE. (1987), Direct method for TDS stability analysis, IEE Pro-
ceeding Part D. 134, 101-107.

[6] Bonnet, C.; Partington, J.R. (2002), Analysis of fractional delay systems of retarded and
neutral type, Automatica, 38(7):1133-1138.

[7] Pakzad, M. A.; Pakzad, S.; Nekoui, M.A. (2013); Stability analysis of multiple time delayed
fractional order systems, American Control Conference, Washington, DC, 170-175.

[8] A. R. Fioravanti, C. Bonnet, H. Ozbay and S. I. Niculescu, (2012), A numerical method for
stability windows and unstable root-locus calculation for linear fractional time-delay systems,
Automatica, 48(11):2824-2830.

[9] Hwang, C.; Cheng, Y.C. (2006), A numerical algorithm for stability testing of fractional delay
systems, Automatica, 42(5): 825-831.

[10] Hwang, C.; Cheng, Y.C. (2005), A note on the use of the lambert w function in the stability
analysis of time-delay systems, Automatica, 41(11):1979-1985.