On the Stability of Fractional-Order Systems of Neutral Type

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.

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

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

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

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.

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.

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.

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

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, a_n(ω) = b_n(ω) = 0;

3. (3)

   all the coefficients in$hℜ$vanish, a_n(ω) = … = a₀(ω) = 0;

4. (4)

   all the coefficients in$hI$vanish, b_n(ω) = … = b₀(ω) = 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.

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

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.

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 + ε.

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ω_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.

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 jω_m regardless of τ_mk which creates it.

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

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

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.

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¹) to τ = 2.4 dark gray (dark red when viewed in color¹) 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²) to τ = 3.5 dark gray (dark red when viewed in color²).

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.

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