Stability of Fractional Order Systems

1 Departamento de Matemática Fundamental, University of La Laguna, Tenerife 38271 La Laguna, Spain 2Department of Economics, Belarusian State University, Nezavisimosti Avenue 4, 220030 Minsk, Belarus 3 Department of Electrical Engineering, Institute of Engineering, Polytechnic of Porto, 4200-072 Porto, Portugal 4Departamento de Análisis Matemático, University of La Laguna, Tenerife 38271 La Laguna, Spain


Classical Stability Analysis
The study of stability of polynomial and related questions for differential equations goes back to XIX century. Hurwitz (or Routh-Hurwitz) criterion [1] is a necessary and sufficient condition for all the roots of a polynomial with real coefficients and > 0 to have negative real parts. It consists of the following: all principal minors Δ , = 1, 2, . . . , , of the Hurwitz matrix are positive.
Here is a matrix of order whose th row is of the form where = 0 if > or < 0. Polynomial ( ) satisfying the Hurwitz condition is called a Hurwitz polynomial or, in applications of the Routh-Hurwitz criterion in the stability theory of oscillating systems, a stable polynomial. Exact and approximate methods of Hurwitz factorization were developed intensively (see, e.g., [2,3]).
Among other criteria concerning zeros distribution of polynomials, we have to mention Mikhailov stability criterion [4]. It states that all roots of a polynomial with real coefficients have strictly negative real part if and only if the complex-valued function = ( ), = √ −1, of a real variable ∈ [0,∞) describes a curve (the Mikhailov hodograph) in the complex -plane which starts on the positive real semiaxis and does not cross the origin and successively generates an anticlockwise motion through quadrants.
An equivalent condition is as follows: the radius vector ( ), as increases from 0 to +∞, never vanishes and monotonically rotates in a positive direction through an angle /2. The Mikhailov criterion gives a necessary and sufficient condition for the asymptotic stability of a linear differential equation of order with constant coefficients or of a linear system (where the prime symbol denotes first derivative with respect to time): with a constant matrix , the characteristic polynomial of which is ( ).
A very general result on the zero location problem for the polynomial is the Hermite-Biehler theorem [5], which states that the roots of ( ) + ℎ( ) are all on the same side of the real axis when ( ) and ℎ( ) are polynomials with real coefficients if and only if the zeros of ( ) and ℎ( ) are real and alternate. The Hermite-Biehler theorem provides necessary and sufficient conditions for Hurwitz stability of real polynomials in terms of an interlacing (alternating) property [6][7][8]. Notice that if a given real polynomial is not a Hurwitz one, then the Hermite-Biehler theorem does not provide information on its roots distribution.
During the last years, several surveys addressed this topic; see, for instance, [9][10][11][12]. This paper, without being exhaustive, is complementary to the contents of such other works. In our paper, we introduced formally a selected set of methods to characterize the stability of fractional systems. It is intended to form a comprehensive text so that readers can follow easily the concepts. Furthermore, the limitations of the known methods are also pointed out, giving readers the opportunity to consider the open problems.
Bearing these ideas in mind, this paper is organized as follows. Sections 2 and 3 introduce fundamental aspects, namely, the concepts of quasi-polynomials and fractional quasi-polynomials, respectively. Section 4 addresses the main core of the paper, the stability of fractional order systems, and is divided into eight subsections. The section starts by presenting general fractional order systems. In the next subsections issues associated with linear time-invariant systems are discussed in more details, such as controllers, positive systems, systems with delay, and distributed, discrete-time, and nonlinear systems. Stability of closed-loop linear control systems becomes a major motivation of the paper. Therefore, four examples in Section 5 deal with this topic. Finally, Section 6 outlines several techniques that are presently emerging and draws the main conclusions.

Quasi-Polynomials
Pontryagin [13] gave a generalization of the Hermite-Biehler Theorem, which appeared to be very relevant formal tool for the mathematical analysis of stability of quasi-polynomials, that is, of the functions of the following type: where ( ) are polynomials in with constant coefficients, and , = 0, . . . , , are real (or complex) numbers. By other words, ( ) is a sum, where the terms are the product of and exponential and polynomial function with constant coefficients. In control theory, such exponentials correspond to delays. If are commensurable real numbers, that is, = ⋅ , = 0, . . . , , and > 0, then the quasi-polynomial (6) can be written in the form with where ( , ) is a polynomial function of two variables, and then, if = , we get ( ). Thus, from this point of view, the determination of the zeros of a quasi-polynomial (7) by means of Pontryagin theorem can be considered to be a mathematical method for analysis of stabilization of a class of linear time invariant systems with time delay (see, e.g., [14]): where are constant coefficients.
Pontryagin Theorem (see [13]). Let ( ) = ( , ) be a quasi-polynomial of the type (7), where ( , ) is a polynomial function in two variables with real coefficients. Suppose that the "oldest" coefficient ̸ = 0. Let ( ) be the restriction of the quasi-polynomial ( ) to imaginary axis. One can express ( ) = ( ) + ( ), where the real functions (of a real variable) ( ) and ( ) are the real and imaginary parts of ( ), respectively. Let one denote by and , respectively, the zeros of the functions ( ) and ( ). If all the zeros of the quasi-polynomial ( ) lie to the left side of the imaginary axis, then the zeros of the functions ( ) and ( ) are real, alternating, and for each ∈ R.
Reciprocally, let one of the following conditions be satisfied: (1) all the zeros of the functions ( ) and ( ) are real and alternate, and the inequality (10) is satisfied for at least one value ; (2) all the zeros of the function ( ) are real, and for each zero of ( ) the inequality (10) is satisfied; that is, ( ) ( ) < 0; (3) all the zeros of the function ( ) are real, and for each zero of ( ) the inequality (10) is satisfied; that is, Then all the zeros of the quasi-polynomial ( ) lie to the left side of the imaginary axis.
In [15] quasi-polynomial of the type is studied where ( ) and ( ) are polynomials with constant coefficients given by Mathematical Problems in Engineering 3 A notion of the principal term, closely connected with the stability problem, is used in this study; namely, the principal term of quasi-polynomial (11) after premultiplying it by is the term in which the argument of the power and has the highest value for some = 0, 1, . . . , . From the Pontryagin criterion, follows that a quasi-polynomial with no principal term has infinitely many roots with arbitrary large, positive real parts. Hence, the presence of principal term in a quasi-polynomial is a necessary condition for its stability. In [15] the following formula, related zeros of quasi-polynomials and it coefficients are proved: Stability of systems of differential equations with delay (or, in other words, systems of differential-difference equations) with a constant matrix was also investigated by using properties of quasi-polynomials.
Thus, in [16] the following criterion was proved: let be a 2 × 2 matrix of real constants. Then all zeros of the quasipolynomial Δ( ) = 2 2 − tr( ) + | | have negative real parts if and only if where is the smallest positive root of the equation sin = −(1/2) tr( ). Distribution of zeros of quasi-polynomials, related to the coupled renewal-differential system is discussed in [17]). A numerical method for calculation of zeros of quasi-polynomials is proposed, for example, in [18]. Special attention was paid in the last years to (finite and infinite) Dirichlet series. In [19] it was proved that if is a convergent Taylor-Dirichlet series, where, as usual, C[ ] means the set of polynomial with constant complex coefficients. The symbol ↑ ∞ denotes 1 < 2 < ⋅ ⋅ ⋅ < < ⋅ ⋅ ⋅ → ∞, when → ∞, and satisfies an algebraic differential-difference equation ( , ( 1 ) ( + ℎ 1 ) , . . . , ( ) ( + ℎ )) = 0, ( , 1 , . . . , ) = ∑ 1 ,..., where 1 ,..., are constant coefficients, and then the set of its exponents { } ∞ =0 has a finite, linear, integral basis. Different questions related to the series of polynomial of exponents were discussed in [20] (see also references therein and [21]). In [22] the large asymptotic of zeros of sections of a generic exponential series are derived, where the th section of the mentioned series mean the the sum of the first terms of it. In [23] it is given an answer on the question when the reciprocal to the product of Gamma-functions coincide with a quasi-polynomial of type (6).

Fractional Quasi-Polynomials
Recently, an attention is paid to the study of linear fractional systems with delays described by the transfer function where is such a real number (0 < ≤ 1), and the fractional degree nontrivial polynomials ( ) and ( ) with real coefficients have the forms where , are real nonnegative numbers and 0 ̸ = 0, 0 ̸ = 0. The fractional degree characteristic quasi-polynomial of the system (19) has the form In the case of a system with delays of a fractional commensurate order (i.e., when = ⋅ ( = 0, 1, . . . , ); = ⋅ ( = 0, 1, . . . , )), one can consider the natural degree quasi-polynomial associated with the characteristic quasi-polynomial (21) of a fractional order.
In [24] new frequency domain methods for stability analysis of linear continuous-time fractional order systems with delays of the retarded type are proposed. The methods 4 Mathematical Problems in Engineering are obtained by generalization to the class of fractional order systems with delays of the Mikhailov stability criterion and the modified Mikhailov stability criterion known from the theory of natural order systems without and with delays.
The following results concerning stability of the considered system are proved in [24].
(1) The fractional quasi-polynomial (21) of commensurate degree satisfies the condition ( ) ̸ = 0, Re ≥ 0 if and only if all the zeros of the associated natural degree quasipolynomial (22) satisfy the condition where arg( ) means the principal branch of the multivalued function Arg( ), ∈ C; that is, arg( ) ∈ (− , ]. (2) The fractional quasi-polynomial (21) of commensurate degree is not stable for any > 1.
(3) The fractional characteristic quasi-polynomial (21) of commensurate degree is stable if and only if which means that the plot of ( ) with increasing from 0 to +∞ runs in the positive direction by quadrants of the complex plane, missing the origin of this plane.
(4) The fractional characteristic quasi-polynomial (21) (of commensurate or noncommensurate degree) is stable if and only if where and ( ) can be chosen, for example, as Remark 1. Stability of fractional polynomials is related to the stability of ordinary quasi-polynomials due to the following relation: where = log . Application of formula (28) needs to be very careful since is now the point on the Riemann surface of the logarithmic function (see, e.g., [22]).
An approach describing the stability of fractional quasipolynomials in terms of zeros distribution of these polynomials on certain Riemann surfaces is widely used now (see, e.g., survey paper [10] and references therein). Sometimes it is called "root-locus method. " The characteristic result obtained with the application of this method is the following [25].
(5) The fractional order system with characteristic polynomial ( ) is stable if and only if ( ) has no zeros in the closed right-half of the Riemann surface; that is, The fractional order polynomial ( ) is a multivalued function whose domain is a Riemann surface. In general, this surface has an infinite number of sheets, and thus the fractional polynomial has (in general) an infinite number of zeros. We are interested only in those zeros which are situated on the main sheet of the Riemann surface which can be fixed in the following way: − < arg( ) < .
Recently the notion of robust stability was introduced for systems with characteristic polynomials dependent on uncertainty parameter (see [26][27][28][29]). In [25] this notion is applied to the convex combination of two fractional degrees polynomials where is uncertainty parameter, and ( ), ( ) are fractional degree polynomials.
The family (30) of fractional degree polynomials is called robust stable if polynomial ( , ) is stable for all ∈ .
Generalization of the Mikhailov-type criterion (see [24]) to this case has the following form [25].
does not cross the nonpositive part (−∞, 0] of the real axis in the complex plane.

Stability of Fractional Order Systems
Several applications of the results on stability of the fractional polynomials to the systems describing different processes and phenomena are presented, for example, [30,31]. Some of these results are closely related to the recent achievements and the theory and applications of fractional calculus and fractional differential equations (see [32][33][34]). The results in this area need to be systematized as from the point of the ideas, technical point of view.

General Fractional Order Systems.
A general fractional order system can be described by a fractional differential equation of the form where = 0 denotes the Riemann-Liouville RL 0 or Caputo 0 fractional derivative [32,33]. Another form of general fractional order system is due to properties of the Laplace transform [32,34]. This form represents the system in terms of corresponding transfer function where is the Laplace variable. Here , . . . , 0 , , . . . , 0 are given real constants, and , . . . , 0 , , . . . , 0 are given real numbers (usually positive). Without loss of generality, these sets of parameters can be ordered as If both sets -and -constitute an arithmetical progression with the same difference, that is, = , = 0, . . . , , = , = 0, . . . , , then system (33) is called commensurate order system. Usually it is supposed that parameter satisfies the inequality 0 < < 1. In all other cases system (33) is called incommensurate order system. Anyway, if parameters and are rational numbers, then this case can be considered as commensurate one, with = 1/ being a least common multiple of denominators of fractions , . . . , 0 , , . . . , 0 [35]. For commensurate order system, its transfer function can be thought as certain branch of the following multivalued function: Since the right hand-side of this relation is a rational function of , then one can represent ( ) in the form of generalized simple fractions. The most descriptive representation of such a type is that for > , where − is a root of polynomial ( ) of multiplicity . In particular, if all roots are simple, then the representation (36) has the most simple form In this case an analytic solution to system (33) is given by the formula where the symbol " * " means the Laplace-type convolution, and ,] is the two-parametric Mittag-Leffler function [36] ,] ( ) = In the case of homogeneous fractional order system the analytical solution is given by the following formula (see, e.g., [10,37]): is the th derivative of two-parametric Mittag-Leffler function [32,38] ,] ( ) = The stability analysis of the fractional order system gives the following results (see [10,24,39,40]).

Theorem 2. A commensurate order system with transfer function (37) is stable if and only if
with − being the th root of the generalized polynomial ( ).
To formulate the result in incommensurate case, we use the concept of bounded input-bounded output (BIBO) or external stability (see [10,39]).

Theorem 3. Let the transfer function of an incommensurate order system be represented in the form
for some complex numbers , , positive , and positive integer .
Such system is BIBO stable if and only if parameters and arguments of numbers satisfy the following inequality: The result of Theorem 3 was obtained by using the stability results given in [40,41].

Fractional Order Linear Time-Invariant Systems. Besides the conception of stability, for fractional order linear timeinvariant systems
the conceptions of controllability and observability (known as linear and nonlinear differential systems [42,43]) are introduced too. In (47) x ∈ R is an unknown state vector, and u ∈ R , y ∈ R are the control vector and output vector, respectively. Given (constant) matrices A, B, and C are of the following size A ∈ R × , B ∈ R × , and C ∈ R × , . Positive vector q = [ 1 , . . . , ] denotes the (fractional) order of system (47). If 1 = ⋅ ⋅ ⋅ = = , then system (47) is called a commensurate order system.
As in case of ordinary linear differential time-invariant systems, controllability and observability conditions [44] are represented in terms of controllability We have to mention also the stability criterion for the system (47) (see [40,[45][46][47][48]).

Fractional Order Controllers.
The fractional-order controller (FOC) PI D (also known as PI D controller) was proposed in [33] as a generalization of the PID controller with integrator of real order and differentiator of real order . The transfer function of such controller in the Laplace domain has this form where is the proportional constant, is the integration constant, and is the differentiation constant.
In [69] a classification of different modifications of the fractional PI D controllers (see also [49,[70][71][72]): (i) CRONE controller (1st generation), characterized by the bandlimited lead effect: There are a number of real-life applications of three generations of the CRONE controller [73].
(iii) Noninteger integral and its application to control as a reference function [74,75]; Bode suggested an ideal shape of the loop transfer function in his work on design of feedback amplifiers in 1945. Ideal loop transfer function has the form where is desired crossover frequency and is the slope of the ideal cut-off characteristic. The Nyquist curve for ideal Bode transfer function is simply a straight line through the origin with arg( ( )) = /2.
(iv) TID compensator [76], which has structure similar to a PID controller but the proportional component is replaced with a tilted component having a transfer function to the power of (−1/ ). The resulting transfer function of the TID controller has the form where , , and are the controller constants, and is a non-zero real number, preferably between 2 and 3. The transfer function of TID compensator more closely approximates an optimal transfer function, and an overall response is achieved, which is closer to the theoretical optimal response determined by Bode [74].
Different methods for determination of PI D controller parameters satisfying the given requirements are proposed (see, e.g., [69] and references therein).

Positive Fractional Order Systems.
A new concept (notion) of the practical stability of the positive fractional 2D linear systems is proposed in [77]. Necessary and sufficient conditions for the practical stability of the positive fractional 2D systems are established. It is shown that the positive fractional 2D systems are practically unstable (1) if a corresponding positive 2D system is asymptotically unstable and (2) if some matrices of the 2D system are nonnegative.
Simple necessary and sufficient conditions for practical stability independent of the length of practical implementation are established in [78]. It is shown that practical stability of the system is equivalent to asymptotic stability of the corresponding standard positive discrete-time systems of the same order.

Fractional Order Systems with Delay.
Fractional order systems with delay meet a number of important applications (see, e.g., [31,49]). There are several important works about stability of closed-loop fractional order systems/controllers with time delays. Some relevant examples can be found in [79][80][81].
To describe simplest fractional order systems with delay, let us introduce some notations (see [82]).
Let C ([ , ], R ) be the set of continuous functions mapping the interval [ , ] to R ). One may wish to identify a maximum time delay of a system. In this case, we are interested in the set of continuous function mapping [− , 0] to R ), for which we simplify the notation to C = C([− , 0], R ). For any > 0 and any continuous function of time x ∈ C([ 0 − , 0 + ], R ), 0 ≤ 0 + , let x ( ) ∈ C be a segment of function defined as x ( ) = x( + ), − ≤ ≤ 0.
Let the fractional nonlinear time-delay system be the system of the following type: where x ∈ C([ 0 − , 0 + ], R ) for any > 0, q = ( , . . . , ), 0 < < 1, and : R × C → R . To determine the future evolution of the state, it is necessary to specify the initial state variables x( ) in a time interval of length , say from 0 − to 0 ; that is, where ∈ C is given. In other words x( 0 )( ) = ( ), − ≤ ≤ 0. Several stability results for fractional order systems with delay were obtained in [82][83][84]. In particular, in [84], the linear and time-invariant differential-functional Caputo fractional differential systems of order are considered: − 1 < ≤ , − 1 ∈ Z 0+ , 0 < ℎ 0 < ℎ 1 < ⋅ ⋅ ⋅ < ℎ = ℎ < ∞. A 0 , A ∈ R × are matrices of dynamics for each delay ℎ , and B ∈ R × is the control matrix. Under standard initial conditions, the solution to this problem is represented via Mittag-Leffler function, and the dependence of the different delay parameters is studied.

Distributed Order Fractional
Systems. An example of distributed order fractional systems is the system of the following type (see, e.g., [85]) containing the so-called distributed order fractional derivative: where 0 < < ≤ 1, and (⋅) = (⋅)/ is a standard (single order) fractional derivative.
Application of such systems to the description of the ultraslow diffusion is given in series of articles by Kochubei (see, e.g., [86]).

Discrete-Time Fractional Systems.
There are different definitions of the fractional derivative (see, e.g., [34,87]). The Grünwald-Letnikov definition, which is the discrete approximation of the fractional order derivative, is used here. The Grünwald-Letnikov fractional order derivative of a given function ( ) is given by where the real number denotes the order of the derivative, is the initial time, and ℎ is a sampling time. The difference operator Δ is given by where ( ) is the Pochhammer symbol and [⋅] denotes an integer part of a number. Traditional discrete-time state-space model of integer order, that is, when is equal to unity has the form, where u( ) ∈ R and y( ) ∈ R are, respectively, the input and the output vectors, and x( ) = [ 1 ( ), . . . , ( )] ∈ R is the state vector. Its initial value is denoted by 0 = (0) and can be set equal to zero without loss of generality. A, B, C, and D are the conventional state space matrices with appropriate dimensions.
The generalization of the integer-order difference to a noninteger order (or fractional-order) difference has been addressed in [88] where the discrete fractional-order difference operator with the initial time taken equal to zero is defined as follows: In the sequel, the sampling time ℎ is taken equal to 1. These results conducted to conceive the linear discrete-time fractional-order state-space model, using The discrete-time fractional order system is represented by the following state space model: Mathematical Problems in Engineering where A 0 = A − 1 I and A = − +1 for ≥ 2 with = −(−1) ( ). This description can be extended to the case of noncommensurate fractional-order systems modeled in [88] by introducing the following vector difference operator: ] .
Stability analysis of such system is performed, for example, [89] (see also [90,91]). The paper [92] is devoted to controllability analysis of discrete-time fractional systems.

Fractional Nonlinear Systems.
The simplest fractional nonlinear system in the incommensurate case is the system of the type 0 q x ( ) = f (x ( ) , ) , It has been mentioned in [39] that exponential stability cannot be used to characterize the asymptotic stability of fractional order systems. A new definition of power law stability was introduced [93].
Power law − asymptotic stability is a special case of the Mittag-Leffler stability [94], which has the following form.
Among the many methods that have been proposed for the study of "different kinds of stability definitions" of nonlinear fractional order systems (66), we mention perturbation analysis.
Thus, in [95], it is investigates the qualitative behaviour of a perturbed fractional order differential equations with Caputo derivatives that differs in initial position and initial time with respect to the unperturbed fractional order differential equation with Caputo derivatives. In [96], the stability of -dimensional linear fractional differential systems with commensurate order and the corresponding perturbed systems is investigated. By using the Laplace transform, the asymptotic expansion of the Mittag-Leffler function, and the Grönwall's inequality, some conditions on stability and asymptotic stability are given. In [97], the stability of nonlinear fractional differential systems with Caputo derivatives by utilizing a Lyapunov-type function is studied. Taking into account the relation between asymptotic stability and generalized Mittag-Leffler stability, the condition on Lyapunov-type function is weakened.

Some Examples
In this section, we analyze the stability of the four systems by means of the root locus and the polar diagram. For calculating the root-locus, the algorithm proposed in [98] is adopted. The closed-loop system is constituted by a controller and a plant with transfer functions ( ) and ( ), respectively, and unit feedback.
For studying the stability, the following classical criteria are applied

Future Directions of Research and Conclusions
In the previous discussion, one can outline the main methods applied at the study of stability of ordinary and fractional order systems. These directions are as follows: (i) complex analytic methods related to the properties of single-and multivalued analytic functions; (ii) methods of geometric functions theory describing the behaviour of polynomials or systems in geometrical terms; (iii) methods of linear algebra (mainly matrix analysis); (iv) methods of stochastic analysis; (v) methods of fuzzy data analysis; (vi) perturbation analysis of nonlinear systems; (vii) methods of differential equations of fractional order.      We can also mention the possibility to apply the study of fractional order stability methods of Padé (or Hermite-Padé) approximation (see, e.g., [100][101][102]).
In conclusion, this paper reviewed the main contributions that were proposed during the last years for analysing the stability of fractional order systems. Different problems were addressed such as control, systems including a delay or with a distributed nature, as well as discrete-time and nonlinear systems. The paper presented in a comprehensive and concise way many details that are scattered in the literature and provide researchers a reference text for work in this area.