## 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

*α*< 0 and

*α*> 0, respectively

*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]:

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

*n*is the first integer which is not less than

*α*. Laplace transform of the Riemann–Liouville derivative is given as follows:

*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:

*m*is the first integer that is larger than

*α*. The Laplace transform of the Caputo fractional derivative is

*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

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

*x*(

*n*× 1), $A,B,C\u2208Rn\xd7n$ and $\tau \u2208R+$. 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\alpha 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

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 $\tau \u2208R+$ and $P\u2113(s\alpha )$ for $\u2113\u2208NN$ 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 $\u211c(s)\u22650$ (in particular, no poles of fractional order at

*s*= 0).

*s*= ±

*jω*or in other words, $s=\omega e\xb1j\pi /2$ is the root of characteristic equation (11) for a $\tau \u2208R+$. Then for

*s*, the roots are defined as follows:

^{α}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.

^{–}

*for purely imaginary roots*

^{τs}*s*= ±

*jω*. Moreover, transformation (14) is different from the first-order Pad'e approximation of e

^{–}

*which is*

^{τs}*s*= ±

*jω*, 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

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*=

*jω*so long as Eq. (17) holds.

*Ts*)

*, the new form of the characteristic equation is reached*

^{n}*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*=

*jω*should be inserted into relation (19) and then the real and imaginary parts of the resulting equation should be separated as follows:

*ω*,

*T*) solutions from

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

*Consider the two multivariate polynomials (21) and (22) in terms of ω, T with real coefficients, where*$h\u211c$

*and*$hI$

*have positive degrees in terms of T, and n*

*>*

*0. The resultant of*$h\u211c$

*and*$hI$

*with respect to T is defined by*

*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\u211c,hI)=0$*for ω. Conversely, if*$RT(h\u211c,hI)=0$*for some ω, then at least one of the following four conditions
holds:*

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.

*n*-order Sylvester matrix is constructed via (23), and its determinant $RT(h\u211c,hI)$ is a function of

*ω*.

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ω*_{c}*for all possible*$\tau \u2208R+$.

*Proof.* Let us assume *s* = *jω*_{c} be a pair of roots for *C*(*s*, *τ*) then *s* = *jω*_{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.

*ω*values, for which

*s*=

*jω*is a root of Eq. (12) for some non-negative delays, are defined as the crossing frequency set

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

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

*ω*

_{c}

*,*

_{m}*m*= 1, 2,…,

*n*correspond infinitely many, periodically spaced

*τ*values. All this set

where *τ _{m}*

_{,}

_{k}_{+1}−

*τ*

_{m}_{,}

*= 2*

_{k}*π*/

*ω*is the apparent period of repetition. After finding cross points of characteristic equation (12), the direction of cross frequency transition

_{m}*ω*

_{c}

*in*

_{m}*τ*should be determined as

_{mk}*τ*increases from

_{mk}*τ*−

_{mk}*ε*to

*τ*+

_{mk}*ε*.

### Direction of Crossing.

*s*,

*τ*) is a simple root of

*C*(

*s*,

*τ*) = 0. The root sensitivities associated with each purely imaginary characteristic root

*jω*with respect to

*τ*are defined as

*ω*

_{c}

*and*

_{m}*τ*is defined as

_{mk}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 *jω*_{c}* _{m}* 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.*The simple roots of Eq. (12) are continuously differentiable with respect to $\tau \u2208R+$. Then one can find

*ds*/

*dτ*for simple roots of Eq. (12) as follows:

*τ*is written as follows:

_{mk}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

*τ*

_{m}_{1}. Therefore, the RT remains the same at a given crossing point

*jω*regardless of

_{m}*τ*which creates it.

_{mk}## 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:

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,
it is easy to find out that a destabilizing crossing of roots (RT = +1) has occurred
at *τ* = 2.1337761636 + 4.9478442469*k* for *s* = ±*j*1.2698834065 for all values of $k\u2208Z+$.
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.

τ | ω (rad/s) | Root tendency | Number of unstable roots |
---|---|---|---|

0 | |||

0 | |||

2.1337761636 | 1.2698834065 | + | |

2 | |||

7.0816204105 | 1.2698834065 | + | |

4 | |||

12.0294646575 | 1.2698834065 | + | |

6 | |||

16.9773089044 | 1.2698834065 | + | |

⋮ | ⋮ | ⋮ | ⋮ |

τ | ω (rad/s) | Root tendency | Number of unstable roots |
---|---|---|---|

0 | |||

0 | |||

2.1337761636 | 1.2698834065 | + | |

2 | |||

7.0816204105 | 1.2698834065 | + | |

4 | |||

12.0294646575 | 1.2698834065 | + | |

6 | |||

16.9773089044 | 1.2698834065 | + | |

⋮ | ⋮ | ⋮ | ⋮ |

*Example*

*2*

*.*Consider the following fractional-order system of neutral type with order 0.5 and commensurate delays:

*τ*= 0. Therefore, the system is stable without delay. And in view of relation (19), we can eliminate exponential term from Eq. (41) as follows:

*s*=

*jω*and $s0.5=\omega 0.5(22+j22)$ into Eq. (42) and equating the real and imaginary parts of the obtained relation to zero, we get

*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)

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.07981199*k* for *s* = ±*j*3.0210352253 and a destabilizing
crossing (RT = +1) has taken place at *τ* = −0.404762327 + 1.267265993 for *s* = ±*j*4.9580635306 for all values of $k\u2208Z+$.
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.

τ | ω (rad/s) | Root tendency | Number of unstable roots |
---|---|---|---|

0 | |||

0 | |||

0.862503666 | 4.9580635306 | + | |

2 | |||

1.53306031 | 3.0210352253 | − | |

0 | |||

2.129769659 | 4.9580635306 | + | |

2 | |||

3.397035652 | 4.9580635306 | + | |

4 | |||

3.61287230 | 3.0210352253 | − | |

2 | |||

4.664301645 | 4.9580635306 | + | |

4 | |||

⋮ | ⋮ | ⋮ | ⋮ |

τ | ω (rad/s) | Root tendency | Number of unstable roots |
---|---|---|---|

0 | |||

0 | |||

0.862503666 | 4.9580635306 | + | |

2 | |||

1.53306031 | 3.0210352253 | − | |

0 | |||

2.129769659 | 4.9580635306 | + | |

2 | |||

3.397035652 | 4.9580635306 | + | |

4 | |||

3.61287230 | 3.0210352253 | − | |

2 | |||

4.664301645 | 4.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 color^{1}) to *τ* = 2.4 dark gray (dark
red when viewed in color^{1}) 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 color^{2}) to *τ* = 3.5 dark gray (dark red when viewed in color^{2}).

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

- $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

- $\u211c(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\xaf$ =
complex conjugate of a complex number

*z*- $|z|,\u2220z$ =
magnitude and argument of a complex number

*z*