## Abstract

The aim of this study is to offer a new analytical method for the stability testing of neutral type linear time-invariant (LTI) time-delayed fractional-order systems with commensurate orders and multiple commensurate delays. It is evident from the literature that the stability assessment of this class of dynamics remains unsolved yet and this is the first attempt to take up this challenging problem. The method starts with the determination of all possible purely imaginary characteristic roots for any positive time delay. To achieve this, the Rekasius transformation is used for the transcendental terms in the characteristic equation. An explicit analytical expression in terms of the system parameters which reveals the stability regions (pockets) in the domain of time delay has been presented. The number of unstable roots in each delay interval is calculated with the definition of root tendency (RT) on the boundary of each interval. Two example case studies are also provided, which are not possible to analyze using any other methodology known to the authors.

## Introduction

A lot of systems in nature have inherent delays in their functioning and this phenomenon also occurs in many industrial systems and processes, such as microwave oscillations, chemical reactors, dynamic models of population, thermal systems, long transmission lines, communication systems, etc. Moreover, the existence of various sensors and actuators in the feedback loop is the cause of delay in many systems, including the wireless and web-based control systems. Delay can appear in the input, output, or in the state variables. The existence of delay in input has no effect on system stability as long as the system has an open-loop form. However, if such systems are controlled via feedback, the delay will be transferred to the closed-loop system's characteristic equation and will cause it to have infinite-dimensional system and, then, it will have a considerable impact on the stability of the system. Therefore, time-delayed systems play significant roles in theoretical as well as practical fields; and this influence can be observed in numerous research articles written on various problems that involve this class of systems [1–8]. Time-delay systems can be divided into two types: retarded or neutral system. The retarded system contains delays only in its states, whereas the neutral system contains delays both in its states and in the derivatives of its states, so-called neutral delays.

The stability of every dynamic system is a basic question and a fundamental issue. In confronting the time-delay systems, we are curious to know what will happen if the amount of delay increases and how the stability feature will change. Regarding the systems with a single delay or with multiple commensurate delays, an interesting phenomenon called the “stability windows” [9] or “stability pockets” [10] may occur. And regarding the systems with multiple independent time delays, in view of the number of independent delays present, the stability map of the system can be expressed as stable and unstable regions in the 2D or 3D space [11,12]. In the LTI integer order systems, some useful time-domain methods [13,14] as well as frequency-domain methods [15,16] have been presented with the purpose of evaluating the stability of these dynamics. Regarding the fractional-order time-delay systems, due to the involvement of fractional mathematics in these problems, the stability analysis of these systems will be much more complicated than the integer order systems.

From the perspective of fractional-order systems of retarded type, the researchers of Refs. [17] and [18] may be regarded as vanguards and forerunners in the stability investigation of fractional-order time-delayed systems. They have developed the Ruth Hurwitz criterion for the evaluation of stability in certain delay systems in which the exponent of the involved fractional is $s$. From the numerical analysis point of view, the effective numerical algorithms have been discussed in Refs. [9] and [19] for the evaluation of bounded-input bounded-output (BIBO) stability of fractional-order delay systems. In Ref. [20], a heavy computation scheme based on the Cauchy's integral has been proposed to test the stability of such systems. The stability of these systems can also be analyzed analytically by employing the Hassard's theorem proposed by Shi and Wang [21]. Based on the method presented by them, the number of unstable poles of the characteristic equation for a given delay is determined. The problem with this approach is that the amount of delay must be known so that the stability of a system can be analyzed; thus, the method proposed by the mentioned authors cannot detect the stability windows of a system. The root-locus (RL) method is one of the most popular and powerful tools for both stability analysis and design of LTI systems. In Ref. [22], an efficient algorithm for obtaining the RL of fractional-order expressions of any type has been introduced. Pakzad and Nekoui in Ref. [23] have successfully extended the direct analytical method (presented by Walton and Marshal [16] for the stability analysis of integer order time-delay systems) to fractional-order delay systems. In addition, they have extended an analytical algorithm based on “advanced clustering with frequency sweeping” for testing the stability of fractional-order systems with multiple time delays [24] and with uncertain parameters in both time-delay space and coefficient space in Ref. [25]. Recently, Pakzad et al. in Ref. [26] and Hua et al. in Ref. [27] have presented an analytical method for finding the stability regions of fractional delay systems of retarded type with commensurate orders, which uses the bilinear Rekasius transformation to eliminate the exponential type transcendental term in the characteristic equations of these systems.

The main contribution of this document is that we demonstrate for the first time that the stability maps of a neutral type fractional-order system with time delays can be obtained efficiently. We present one of the important methods of the frequency domain for the stability analysis of fractional-order delay systems. Using the approach presented in this study, first, the transcendental terms have been eliminated from the characteristic equation and then, all the locations where roots cross the imaginary axis are determined. The number of unstable roots in each interval is calculated with the definition of RT on the boundary of each interval. Finally, the concept of stability is expressed as a function of delay. Prior to the conclusion, two examples are presented to confirm the effectiveness of the proposed method results.

## Preliminaries and Definitions

Operator $aDtα$ denotes fractional integral and fractional derivative operators where α < 0 and α > 0, respectively
$aDtα={dαdtαα>01α=0∫at(dτ)-αα<0$
(1)
and a and t are the start and end values of the integral, respectively. The αth-order fractional integral of function x(t) is defined as follows [28]:
$0Dt-αx(t)=1Γ(α)∫0t(t-τ)α-1x(τ)dτ, 0<α∈R$
(2)
where Γ is the Gamma function
$Γ(α)=∫0ttα-1e-tdt$
(3)

For fractional derivative, three definitions have been presented in Ref. [28]. These definitions are as follows:

Grunwald–Letnikov definition
$0Dtαx(t)=limh→0nh=th-α∑j=0n(-1)j(αj)x(t-jh)$
(4)
Riemann–Liouville definition
$0Dtαx(t)=1Γ(n-α)dndtn∫0t(t-τ)n-α-1x(τ)dτ$
(5)
where n is the first integer which is not less than α. Laplace transform of the Riemann–Liouville derivative is given as follows:
$L{dαx(t)dtα}=sαL{x(t)}-∑k=0m-1skdα-k-1x(0)dtα-k-1, m-1≤α≤m$
(6)
Unfortunately, the Riemann–Liouville fractional derivative appears unsuitable to be treated by the Laplace transform technique in that it requires knowledge of the noninteger order derivatives of the function at t = 0. The mentioned problem does not exist in the Caputo definition of the fractional derivative [29]. The Caputo definition of the fractional derivative, which sometimes is called smooth fractional derivative, is described as follows:
$0Dtαx(t)={1Γ(m-α)∫0tx(m)(τ)(t-τ)α+1-mdτ m-1<α
(7)
where m is the first integer that is larger than α. The Laplace transform of the Caputo fractional derivative is
$L{dαx(t)dtα}=sαL{x(t)}-∑k=0m-1sα-k-1x(k)(0), m-1≤α≤m$
(8)
Contrary to the Laplace transform of the Riemann–Liouville fractional derivative, only integer order derivatives of function x appear in the Laplace transform of the Caputo fractional derivative. For zero initial conditions, the Laplace transform of the fractional derivative (8) reduces to
$L{dαx(t)dtα}=sαL{x(t)}$
(9)

In the rest of this paper, the notation $Dtα$ represents the Caputo fractional derivative of order α.

A general class of neutral type fractional-order LTI (linear time invariant) systems with multiple commensurate delays is taken into account
$Dtαx(t)=Ax(t)+Bx(t-τ)+CDtαx(t-τ)$
(10)
where x(n × 1), $A,B,C∈Rn×n$ and $τ∈R+$. Taking A, B, and C as constant matrices, the single parameter that influences the stability of Eq. (10) is the time delay, τ. Notice that the highest order dynamics, $Dtαx$, is time delayed here. This is the attribute of the neutral time-delay systems. For the retarded time-delay systems, this property does not hold, i.e., C = 0. Clearly, the delay injects exponential transcendentality to the characteristic equation. The characteristic equation of the system in Eq. (10) is
$C(s,τ)=det(sαI-A-Be-τs-Csαe-τs)=0$
(11)
which yields a scalar equation in the general form of
$C(s,τ)=Pn(sα)e-nτs+Pn-1(sα)e-(n-1)τs+…+P0(s)=∑ℓ=0nPℓ(sα)e-ℓτs$
(12)

In general case, the stability of Eq. (10) is studied over its characteristic equation given by Eq. (11), where parameter τ is non-negative, such that $τ∈R+$ and $Pℓ(sα)$ for $ℓ∈NN$ are real polynomials in the complex variable sα. We find out from Ref. [30] that a characteristic equation in the form of Eq. (12) will be H stable if, and only if, it does not have any pole at $ℜ(s)≥0$ (in particular, no poles of fractional order at s = 0).

Let us assume that s = ± or in other words, $s=ωe±jπ/2$ is the root of characteristic equation (11) for a $τ∈R+$. Then for sα, the roots are defined as follows:
$sα=ωαe±jαπ/2=ωα(cos(απ2)±jsin(απ2))$
(13)

There is a direct relation between the roots on the imaginary axis for the s-domain with the ones having argument ± απ/2 in the sα - domain.

## Methodology

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-order delay system of neutral type. One of the most common methods for the stability analysis of a linear delay system is to examine the locations of the roots of the system's characteristic equation. In systems whose characteristic equation has a delay term, it would be impossible to solve the equation directly. In dealing with time-delay systems, one of the well-known methods of the frequency domain for the stability analysis of these systems is to find the points at which the poles cross the imaginary axis and then to examine the root tendencies at these point to see whether the roots have a tendency to move to the right side of the imaginary axis (which is called a destabilizing crossing) or to the left side of the imaginary axis (which is called a stabilizing crossing). Note that the expressions of “destabilizing crossing” and “stabilizing crossing” indicate that a pair of poles has crossed the imaginary axis in a defined direction, and not that the system is turning unstable or stable, respectively. For that, it is necessary to know the number of unstable poles before the crossings.

### Crossing Position.

In order to find the points on the imaginary axis, where the crossing takes place, the bilinear Rekasius substitution is used and as a result, the transcendental characteristic equation is converted into an algebraic equation. System (10) will be asymptotically stable if, and only if, all its eigenvalues are on the open left-half complex plane. Characteristic equation (12) is transcendental. To simplify this characteristic equation, one can convert the infinite-dimensional characteristic equation (12) to a finite-dimensional characteristic equation with continuous coefficients in $T∈R$ as proposed in Ref. [26]. This conversion is established by the Rekasius transformation [31], which has been defined as follows:
$e-τs=1-Ts1+Ts, T∈R$
(14)
It is important to note that this substitution is an exact expression of eτs for purely imaginary roots s = ±. Moreover, transformation (14) is different from the first-order Pad'e approximation of eτs which is
$e-τs≈1-0.5τs1+0.5τs$
(15)
For multiple commensurate delays, the transformation (14) can be modified as follows:
$e-ℓτs=(1-Ts1+Ts)ℓ, T∈R$
(16)
Note that for s = ±, the magnitudes of both sides in Eq. (14) always agree with each other. Furthermore, for the transformation to exactly hold, it is required that the phase condition developed from Eq. (14) holds. This condition can be found as
$τ=2ω(tan-1(Tω)+kπ), k=0,1,…$
(17)

Equation (17) is the inverse mapping from (T, ω) domain to τ. This equation describes an asymmetric mapping in which one T is mapped into countably infinite τk which are distributed with a periodicity of 2π/ωc. Note that there is a one-to-one mapping between T and ω. To summarize, the transformation in Eq. (14) is indeed exact at s =  so long as Eq. (17) holds.

By inserting Eq. (16) into Eq. (12), we have
$C(s,T)=∑ℓ=0nPℓ(sα)(1-Ts1+Ts)ℓ$
(18)
By multiplying Eq. (18) by (1 + Ts)n, the new form of the characteristic equation is reached
$h(s,T)=(1+Ts)nC(s,T)=∑ℓ=0nPℓ(sα)(1-Ts)ℓ(1+Ts)(n-ℓ)$
(19)
This expression is a equation in s of which the coefficients are parametric functions of T. As is observed, characteristic equation (12), which had transcendental terms, has been converted into algebraic equation (19). To find the crossing frequencies from the imaginary axis in Eq. (12), s =  should be inserted into relation (19) and then the real and imaginary parts of the resulting equation should be separated as follows:
$h(s,T))|s=jω=hℜ(ω,T)+jhI(ω,T)=0$
(20)
In the above relations, $hℜ=ℜ(h),hI=I(h)$. Making the real and imaginary parts of Eq. (20) equal to zero results in two algebraic equations; let us investigate those (ω, T) solutions from
$hℜ=∑i=0nai(ω)Ti=0$
(21)
and
$hI=∑i=0nbi(ω)Ti=0$
(22)

We utilize the resultant theory to eliminate T from the two multivariate polynomials $hℜ$ and $hI$.

Definition 1.Consider the two multivariate polynomials (21) and (22) in terms of ω, T with real coefficients, where$hℜ$and$hI$have positive degrees in terms of T, and n>0. The resultant of$hℜ$and$hI$with respect to T is defined by
$RT(hℜ,hI)=|an(ω)an-1(ω)………00an(ω)an-1(ω)……0⋮⋮⋮⋮⋮⋮0……⋮a1(ω)a0(ω)bn(ω)bn-1(ω)………00bn(ω)bn-1(ω)……0⋮⋮⋮⋮⋮⋮0……⋮b1(ω)b0(ω)|$
(23)

which is the determinant of the well-known Sylvester matrix [32 ].

Theorem 1 [32]. If (ω, T) is a common root of Eqs.(21)and(22), then$RT(hℜ,hI)=0$for ω. Conversely, if$RT(hℜ,hI)=0$for some ω, then at least one of the following four conditions holds:

1. (1)

there exists (ω, T) that is a common root of both Eqs.(21)and(22);

2. (2)

leading coefficients of both$hℜ$and$hI$vanish, an(ω) = bn(ω) = 0;

3. (3)

all the coefficients in$hℜ$vanish, an(ω) = … = a0(ω) = 0;

4. (4)

all the coefficients in$hI$vanish, bn(ω) = … = b0(ω) = 0.

Detection of the common roots of Eqs. (21) and (22) corresponds to case (1) in Theorem 1, and the remaining cases (2) and (3) can be identified for a given (ω, T) double.

A 2n-order Sylvester matrix is constructed via (23), and its determinant $RT(hℜ,hI)$ is a function of ω.
$D(ω)=RT(hℜ,hI)=0$
(24)

hence the set of all ω and T which makes both equations zero can be obtained. And then for every ω and T, its corresponding τ set can be determined using Eq. (17).

Theorem 2.A given time-delayed fractional system; characteristic equation(12)can exhibit only a finite number of purely imaginary characteristic roots ±jωcfor all possible$τ∈R+$.

Proof. Let us assume s = c be a pair of roots for C(s, τ) then s = c would be a pair of roots for D(ω). since D(ω) is a finite degree polynomial with maximum degree of max(deg(h(s, T))) then the number of crossing points of Eq. (12) that are the real roots of D(ω) is finite.

The whole ω values, for which s =  is a root of Eq. (12) for some non-negative delays, are defined as the crossing frequency set
$Ω={ω∈R+|C(s1/α,τ)=0, for some τ∈R+}$
(25)

Corollary 1.If the system given as Eq.(10)is stable for τ=0 (i.e., system without delay) and$Ω=φ$, then the system will be stable for all positive values of τ.

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

This class of systems exhibits only a finite number of possible imaginary characteristic roots for all $τ∈R+$ at given frequencies. And this method detected all of them. Let us call this set
${ωc}={ωc1,ωc2,ωc3,…,ωcn}$
(26)
where subscript c refers to “crossing” the imaginary axis. Furthermore to each ωcm, m = 1, 2,…,n correspond infinitely many, periodically spaced τ values. All this set
${τm}={τm1,τm2,τm3,…,τm∞}, m=1,2,…,n$
(27)

where τm,k+1 − τm,k = 2π/ωm is the apparent period of repetition. After finding cross points of characteristic equation (12), the direction of cross frequency transition ωcm in τmk should be determined as τmk increases from τmkε to τmk + ε.

### Direction of Crossing.

After the crossing points of characteristic equation (12) 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. As it was shown in Ref. [26], this is constant with respect to subsequential crossings (in Eq. (27)), and therefore, it is denoted as RT. Assume that (s, τ) is a simple root of C(s, τ) = 0. The root sensitivities associated with each purely imaginary characteristic root with respect to τ are defined as
$Sτs|s=jωc=dsdτ|s=jωc=-∂C/∂τ∂C/∂s|s=jωc$
(28)
The RT for each ωcm and τmk is defined as
$Root tendency=RT|s=jωc=sign(ℜ(Sτs|s=jωcmτ=τmk))$
(29)

If it is positive, then it is a destabilizing crossing, whereas if it is negative, this means a stabilizing crossing. In case the result is 0, a higher order analysis is needed, since this might be the case where the root just touches the imaginary axis and returns to its original half-plane. Notice that RT represents the direction of transition of the roots at cm as τ increases from τmkε to τmkε, 0 ε  ≪ 1.

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

Proof. The simple roots of Eq. (12) are continuously differentiable with respect to $τ∈R+$. Then one can find ds/ for simple roots of Eq. (12) as follows:
$dsdτ=-∂C/∂τ∂C/∂s=∑ℓ=0nPℓ(sα)ℓse-ℓτs∑ℓ=0nPℓ(sα)dse-ℓτs-∑ℓ=0nPℓ(sα)ℓτse-ℓτs$
(30)
Based on Eq. (30) and definition given in Eq. (29), the root tendency of each time-delay τmk is written as follows:
$RT|s=jωcτ=sign(ℜ(Sτs|s=jωc))=sign(ℜ(-∂C/∂τ∂C/∂s|s=jωc))=sign(ℜ(∑ℓ=0nPℓ(sα)ℓse-ℓτs∑ℓ=0nPℓ(sα)dse-ℓτs-Pℓ(sα)ℓτe-ℓτs))|s=jωcτ=τm1+2kπωc=sign(ℜ(∑ℓ=0nPℓ(sα)dse-ℓτs∑ℓ=0nPℓ(sα)ℓse-ℓτs-τs)-1)|s=jωcτ=τm1+2kπωc=sign(ℜ(∑ℓ=0nPℓ(sα)dse-ℓτs∑ℓ=0nPℓ(sα)ℓse-ℓτs))|s=jωcτ=τm1+2kπωc=sign(I(∑ℓ=0nPℓ(sα)dse-ℓτs∑ℓ=0nPℓ(sα)ℓe-ℓτs))|s=jωcτ=τm1+2kπωc$
(31)

The proof can easily be completed by observing the fact that expression in Eq. (31) is invariant with respect to the τmk values, which are obtained from the same τm1. Therefore, the RT remains the same at a given crossing point m regardless of τmk which creates it.

## Illustrative Examples

We present two example cases, which display all the features discussed in the text.

Example 1. A LTI fractional delay system of neutral type is given as follows:
$Dt0.9x(t)=-x(t)-1.5x(t-τ)-0.5Dt0.9x(t-τ)$
(32)
The corresponding quasi-polynomial characteristic equation is given by
$C1(s,τ)=(s0.9+1)+(0.5s0.9+1.5)e-τs$
(33)

This system has no unstable pole for τ = 0 (in fact, it has no pole in the physical Riemann sheet). Our objective in this example is to find all the stability windows based on the method described in this article for this system.

By applying the criterion expressed in Sec. 3, we can eliminate exponential term from Eq. (33) as follows:
$h(s,T)=(s0.9+1)(1+Ts)+(s0.9+1)(1-Ts)=0.5Ts1.9-0.5Ts+1.5s0.9+2.5=0$
(34)
And, since $s=jω=ωejπ2$ we will have
$s1.9=ω1.9ej1.9π2=ω1.9cos(1.9π2)+jω1.9sin(1.9π2)s0.9=ω0.9ej0.9π2=ω0.9cos(0.9π2)+jω0.9sin(0.9π2)$
(35)
By inserting expression s =  and Eq. (35) into Eq. (34) and equating the real and imaginary parts of the obtained relation to zero, we get
$hℜ=0.5ω1.9cos(1.9π2)T+1.5ω0.9cos(0.9π2)+2.5=0hI=(0.5ω1.9sin(1.9π2)-0.5ω)T+1.5ω0.9sin(0.9π2)=0$
(36)
Using homomorphism resultant algorithm [32] eliminates T from $hℜ$ and $hI$ and cross points of Eq. (36) are calculated with the real solutions of resulting equation for ω
$D(ω)=ω(0.75ω1.8+0.5ω0.9cos(0.9π2)-1.25)=0$
(37)
The real solutions of Eq. (37) for ω are
$ω=1.2698834065, T=3.5893025094$
(38)
The corresponding infinite countable time delays of the cross points in Eq. (38) are obtained with regards to relation (17)
$τk=2.1337761636+4.9478442469k, k∈Z+$
(39)

By applying the criterion expressed in Sec. 3, it is easy to find out that a destabilizing crossing of roots (RT = +1) has occurred at τ = 2.1337761636 + 4.9478442469k for s = ±j1.2698834065 for all values of $k∈Z+$. Therefore, since the system is stable for τ = 0, the only stability window for this system is 0 ≤ τ < 2. 1337761636. In Table 1, the stability map of system is given. The number of unstable roots in each interval of unstable region has been determined as well. In Fig. 1, the root-loci curve of this system for the changes of τ from τ = 0.4 light gray (dark blue when viewed in color) to τ = 2.4 dark gray (dark red when viewed in color) has been presented for a better perception of the system. Also for some values of the delay, the natural responses of system (32) with initial condition x(t) = 1, −τ < t < 0 are plotted in Fig. 2. It is clear that the stable and unstable regions in Table 1 are in agreement with Fig. 2.

Fig. 1
Fig. 1
Close modal
Fig. 2
Fig. 2
Close modal
Table 1

Stability region (shaded) for example 1

τω (rad/s)Root tendencyNumber of unstable roots
0
0
2.13377616361.2698834065+
2
7.08162041051.2698834065+
4
12.02946465751.2698834065+
6
16.97730890441.2698834065+
τω (rad/s)Root tendencyNumber of unstable roots
0
0
2.13377616361.2698834065+
2
7.08162041051.2698834065+
4
12.02946465751.2698834065+
6
16.97730890441.2698834065+
Example 2. Consider the following fractional-order system of neutral type with order 0.5 and commensurate delays:
$Dt0.5x(t)=[-0.45523.4656-1.49051.7552]x(t) +[-2.79932.75460.399-1.1548]x(t-τ) +[1.6684-2.04281.3108-1.7684]Dt0.5x(t-τ)$
(40)
The corresponding quasi-polynomial is obtained as
$C2(s,τ)=s-1.3s0.5+4.3664+(0.1s+0.1s0.5-1.6647)e-τs +(-0.2728s+0.2279s0.5+2.1335)e-2τs=0$
(41)
The characteristic equation (41) has a pair of poles (s = −5.1547 ±j2.7560) in the left half-plane of the complex plane for τ = 0. Therefore, the system is stable without delay. And in view of relation (19), we can eliminate exponential term from Eq. (41) as follows:
$h(s,T)=(s-1.3s0.5+4.3664)(1+Ts)2 +(0.1s+0.1s0.5-1.6647)(1-Ts)(1+Ts) +(-0.2728s+0.2279s0.5+2.1335)(1-Ts)2$
(42)
By inserting expression s =  and $s0.5=ω0.5(22+j22)$ into Eq. (42) and equating the real and imaginary parts of the obtained relation to zero, we get
$hℜ=(0.586052ω2.5-8.1646ω2)T2 +(-2.5456ω2+1.52792ω1.5)T +(-0.486052ω0.5+4.8352)=0$
(43)
$hI=(-0.6272ω3+0.586052ω2.5)T2 +(-1.52792ω1.5+4.4658ω)T +(0.8272ω-0.486052ω0.5)=0$
(44)
Through the “resultant” command in the Maple software package, T can be eliminated from Eqs. (43) and (44) and the real and positive values of ω can be calculated (the resultant algorithm has been presented in Ref. [21] and in Eq. (23)). And then for the calculated value of ω, the amount of T can be determined from Eqs. (43) and (44). Thus, the time delay set associated with each (ω, T) is obtained with regards to relation (17)
$ω=3.0210352253; T=-0.35895753$
(45)
$ω=4.9580635306; T=-0.316490734$
(46)
The corresponding infinite countable time delays of the cross points in Eqs. (45) and (46) are
$τ=-0.54675168+2.07981199k$
(47)
$τ=-0.404762327+1.267265993k$
(48)

By applying the previously described method, it is easy to find out that a stabilizing crossing of roots (RT = −1) has occurred at τ = −0.54675168 + 2.07981199k for s = ±j3.0210352253 and a destabilizing crossing (RT = +1) has taken place at τ = −0.404762327 + 1.267265993 for s = ±j4.9580635306 for all values of $k∈Z+$. Therefore, we will have two stability windows as follows: 0 ≤ τ < 0.862503666 and 1.53306031 < τ < 2.129769659. In Table 2, we could determine the picture of stability mapping of the system in more detail. Note that at τ = 2.129769659, 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 τ = 3.61287230, another unstable pair of poles goes toward the right half-plane at τ = 3.397035652; and thus, the system cannot recover the stability.

Table 2

Stability regions (shaded) for example 2

τω (rad/s)Root tendencyNumber of unstable roots
0
0
0.8625036664.9580635306+
2
1.533060313.0210352253
0
2.1297696594.9580635306+
2
3.3970356524.9580635306+
4
3.612872303.0210352253
2
4.6643016454.9580635306+
4
τω (rad/s)Root tendencyNumber of unstable roots
0
0
0.8625036664.9580635306+
2
1.533060313.0210352253
0
2.1297696594.9580635306+
2
3.3970356524.9580635306+
4
3.612872303.0210352253
2
4.6643016454.9580635306+
4

To get a better understanding of the properties of this system, its root-loci curve has been plotted as a function of delay in Fig. 3. The color spectrum indicating the changes of τ from τ = 0.1 light gray (dark blue when viewed in color1) to τ = 2.4 dark gray (dark red when viewed in color1) has been illustrated in the “color bar.” One can see that although the depth inside the right half-plane for each destabilizing cross becomes less expressive when τ increases, the stabilizing crossing happens closer to the following destabilizing one, this can be better seen in Fig. 4, where it is shown a zoom around the crossing points through the imaginary axis for the changes of τ from τ = 0.1 light gray (dark blue when viewed in color2) to τ = 3.5 dark gray (dark red when viewed in color2).

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal

## 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 of neutral type is proposed. In this method, the transcendental term was eliminated from the characteristic equation of the original system, and this equation was converted into an algebraic equation. Next, the crossing points through the imaginary axis, and their direction of crossing, were determined. Then, system stability was expressed as a function of delay, based on the information obtained from the system. Finally, two illustrative examples are presented to highlight the proposed approach.

1,2

For interpretation of the references to color in figures 3 and 4, the reader is referred to the web version of this article

### Nomenclature

Nomenclature
$C$ =

set of complex numbers

$j=-1$ =

the imaginary unit

$I(z)$ =

imaginary part of a complex number z

$N,Z$ =

sets of natural and integer numbers

$ℜ(z)$ =

real part of a complex number z

$R(R+,R-)$ =

set of real numbers (larger or equal to zero, smaller or equal to zero)

$z¯$ =

complex conjugate of a complex number z

$|z|,∠z$ =

magnitude and argument of a complex number z

## References

1.
,
M. A.
, and
,
S.
,
2012
, “
Stability Map of Fractional Order Time-Delay Systems
,”
WSEAS Trans. Syst.
,
10
(
11
), pp.
541
550
.
2.
,
M. M.
,
2013
, “
The Use of Generalized Laguerre Polynomials in Spectral Methods for Solving Fractional Delay Differential Equations
,”
ASME J. Comput. Nonlinear Dyn.
,
8
(
4
), p.
041018
.10.1115/1.4024852
3.
Suh
,
Y. S.
,
Ro
,
Y. S.
,
Kang
,
H. J.
, and
Lee
,
H. H.
,
2004
, “
Necessary and Sufficient Stability Condition of Discrete State Delay Systems
,”
Int. J. Control Autom. Syst.
,
2
(
4
), pp.
501
508
.
4.
,
M. A.
,
2012
, “
Kalman Filter Design for Time Delay Systems
,”
WSEAS Trans. Syst.
,
11
(
10
), pp.
551
560
.
5.
Leung
,
A. Y. T.
,
Yang
,
H. X.
, and
Zhu
,
P.
,
2014
, “
Periodic Bifurcation of Duffing-van der Pol Oscillators Having Fractional Derivatives and Time Delay
,”
Commun. Nonlinear Sci. Numer. Simul.
,
19
(
4
), pp.
1142
1155
.10.1016/j.cnsns.2013.08.020
6.
Gu
,
K.
, and
Zheng
,
X.
,
2014
, “
Stability Crossing Set for Systems With Three Scalar Delay Channels
,”
Int. J. Dyn. Control
(submitted)10.1007/s40435-014-0088-3
7.
,
M. A.
, and
Moaveni
,
B.
,
2013
, “
Delay-Dependent State Estimation for Time Delay Systems
,”
WSEAS Trans. Syst.
,
8
(
1
), pp.
1
10
.
8.
Gu
,
K.
,
2013
, “
Complete Quadratic Lyapunov-Krasovskii Functional: Limitations, Computational Efficiency, and Convergence
,”
Advances in Analysis and Control of Time-Delayed Dynamical Systems
, New Jersey World Scientific, pp.
1
19
.10.1142/8878
9.
Fioravanti
,
A. R.
,
Bonnet
,
C.
,
Özbay
,
H.
, and
Niculescu
,
S. I.
,
2012
, “
A Numerical Method for Stability Windows and Unstable Root-Locus Calculation for Linear Fractional Time-Delay Systems
,”
Automatica
,
48
(
11
), pp.
2824
2830
.10.1016/j.automatica.2012.04.009
10.
Olgac
,
N.
, and
Sipahi
,
R.
, “
An Exact Method for the Stability Analysis of Time Delayed Lineartime-Invariant (LTI) Systems
,”
IEEE Trans. Autom. Control
,
47
(
5
), pp.
793
797
.10.1109/TAC.2002.1000275
11.
,
M. A.
,
,
S.
, and
Nekoui
,
M. A.
,
2013
, “
Stability Analysis of Multiple Time Delayed Fractional Order Systems
,”
American Control Conference
,
Washington, DC
, pp.
170
175
.
12.
Sipahi
,
R.
, and
Delice
,
I. I.
, “
Extraction of 3D Stability Switching Hypersurfaces of a Time Delay System With Multiple Fixed Delays
,”
Automatica
,
45
(
6
), pp.
1449
1454
.10.1016/j.automatica.2009.01.017
13.
Suh
,
Y. S.
, and
Lee
,
M. H.
,
1999
, “
Stability of State Delay Systems Based on a Lyapunov Functional
,”
Proceedings of IEEE International Symposium on Industrial Electronics
,
Bled, Slovenia
, Vol.
3
, pp.
1093
1098
.
14.
,
S.
, and
,
M. A.
,
2011
, “
Stability Condition for Discrete Systems With Multiple State Delays
,”
WSEAS Trans. Syst. Control
,
6
(
11
), pp.
417
426
.
15.
Thowsen
,
A.
,
1981
, “
The Routh-Hurwitz Method for Stability Determination of Linear Differential-Difference Systems
,”
Int. J. Control
,
33
(
5
), pp.
991
995
.10.1080/00207178108922971
16.
Walton
,
K. E.
,
,
Marshal
,
J. E.
,
1987
Direct Method for TDS Stability Analysis
,”
IEE Proc., Part D
,
134
, pp.
101
107
.10.1049/ip-d.1987.0018
17.
Öztürk
,
N.
, and
Uraz
,
A.
,
1985
, “
An Analysis Stability Test for a Certain Class of Distributed Parameter Systems With Delays
,”
IEEE Trans. Circuits Syst.
,
34
(
4
), pp.
393
396
.10.1109/TCS.1985.1085704
18.
Jury
,
E. I.
, and
Zeheb
,
E.
,
1986
, “
On a Stability Test for a Class of Distributed System With Delays
,”
IEEE Trans. Circuits Syst.
,
37
(
10
), pp.
1027
1028
.10.1109/TCS.1986.1085839
19.
Buslowicz
,
M.
,
2008
, “
Stability of Linear Continuous Time Fractional Order Systems With Delays of the Retarded Type
,”
Bull. Pol. Acad. Sci.: Tech. Sci.
,
56
(
4
), pp.
319
324
.
20.
Hwang
,
C.
, and
Cheng
,
Y. C.
,
2006
, “
A Numerical Algorithm for Stability Testing of Fractional Delay Systems
,”
Automatica
,
42
(
5
), pp.
825
831
.10.1016/j.automatica.2006.01.008
21.
Shi
,
M.
, and
Wang
,
Z. H.
,
2011
, “
An Effective Analytical Criterion for Stability Testing of Fractional-Delay Systems
,”
Automatica
,
47
(
9
), pp.
2001
2005
.10.1016/j.automatica.2011.05.018
22.
,
J. A.
,
2011
, “
Root Locus of Fractional Linear Systems
,”
Commun. Nonlinear Sci. Numer. Simul.
,
16
(
10
), pp.
3855
3862
.10.1016/j.cnsns.2011.01.020
23.
,
M. A.
, and
Nekoui
,
M. A.
,
2013
, “
Direct Method for Stability Analysis of Fractional Delay Systems
,”
Int. J. Comput. Commun. Control
,
7
(
6
), pp.
863
868
.
24.
,
M. A.
, and
Nekoui
,
M. A.
,
2014
, “
Stability Map of Multiple Time Delayed Fractional Order Systems
,”
Int. J. Control Autom. Syst.
,
12
(
1
), pp.
37
43
.10.1007/s12555-012-0481-7
25.
,
S.
,
,
M. A.
, and
Nekoui
,
M. A.
, “
Stability Map of Fractional Delay Systems in the Parametric Space of Delays and Coefficient
,”
American Control Conference
,
Washington, DC
, pp.
176
181
.
26.
,
M. A.
,
Nekoui
,
M. A.
, and
,
S.
, “
Stability Analysis of Time-Delayed Linear Fractional-Order Systems
,”
Int. J. Control Autom. Syst.
,
11
(
3
), pp.
519
525
.10.1007/s12555-012-0164-4
27.
Hua
,
Ch.
,
Liu
,
D.
, and
Guan
,
X. P.
,
2014
, “
Necessary and Sufficient Stability Criteria for a Class of Fractional-Order Delayed Systems
,”
IEEE Trans. Circuits Syst.
,
61
(
1
), pp.
59
63
.10.1109/TCSII.2013.2291137
28.
Podlubny
,
I.
,
1999
,
Fractional Differential Equations
,
,
New York
.
29.
Tavazoei
,
M. S.
, and
Haeri
,
M.
,
2008
, “
Chaotic Attractors in Incommensurate Fractional Order Systems
,”
Phys. D
,
237
(
20
), pp.
2628
2637
.10.1016/j.physd.2008.03.037
30.
Bonnet
,
C.
, and
Partington
,
J. R.
,
2002
, “
Analysis of Fractional Delay Systems of Retarded and Neutral Type
,”
Automatica
,
38
(
7
), pp.
1133
1138
.10.1016/S0005-1098(01)00306-5
31.
Rekasius
,
Z. V.
,
1980
, “
A Stability Test for Systems With Delays
,”
Proceedings of the Joint Automation Control Conference
,
San Francisco, CA
.
32.
Collins
,
G. E.
,
1971
, “
The Calculation of Multivariate Polynomial Resultants
,”
J. Assoc. Comput. Mach.
,
18
(
4
), pp.
515
532
.10.1145/321662.321666