Stability of Systems of Fractional-Order Differential Equations with Caputo Derivatives

Systems of fractional-order differential equations present stability properties which differ in a substantial way from those of systems of integer order. In this paper, a detailed analysis of the stability of linear systems of fractional differential equations with Caputo derivative is proposed. Starting from the well-known Matignon’s results on stability of single-order systems, for which a different proof is provided together with a clarification of a limit case, the investigation is moved towards multi-order systems as well. Due to the key role of the Mittag–Leffler function played in representing the solution of linear systems of FDEs, a detailed analysis of the asymptotic behavior of this function and of its derivatives is also proposed. Some numerical experiments are presented to illustrate the main results.


Introduction
The investigation of stability properties plays a prominent role in the qualitative theory of fractional-order systems, similarly as in the case of the classical theory of integer-order dynamical systems [1,2]. The classical Hartman-Grobman linearisation theorem, which states that the local behavior of a dynamical system in a neighborhood of a hyperbolic equilibrium is qualitatively equivalent to the behavior of its linearisation near the equilibrium, is extended to the case of fractional-order systems as well [3][4][5]. Consequently, linear stability analysis is of fundamental importance in the investigation of fractional-order systems, and, in particular, stability properties of linear autonomous systems of fractional-order differential equations play a key role in this context.
For single-order systems of fractional differential equations (FDEs), namely systems in which the FDEs have the same fractional order, the most important theoretical result, which may now be considered classical, is Matignon's stability theorem [6], recently generalized in [7] for the case when the fractional order belongs to the interval (0, 2).
Thus far, the investigation of stability properties of multi-order (incommensurate) fractional-order systems has unquestionably received less consideration. We refer to [8][9][10][11] for the stability analysis of incommensurate fractional-order systems with rational orders. Moreover, closely linked to this research topic, bounded input bounded output stability of systems with irrational transfer functions has been investigated in [12,13]. Very recently, the asymptotic properties of solutions of several classes of linear multi-order systems of fractional differential equations (such as systems with block triangular coefficient matrices) have been considered in [14].
The main difficulty in establishing necessary and sufficient conditions for the stability of multi-order linear systems of fractional differential equations (conceivably comparable to the classical Routh-Hurwitz conditions for integer-order systems) is due to the fact that a large number of parameters are involved: the system's coefficients, as well as multiple fractional orders. Undoubtedly, the complexity of the problem is positively correlated with the system's dimension.
The case of two-dimensional multi-order fractional-order systems has been fully investigated in [15][16][17]. On one hand, necessary and sufficient conditions for the asymptotic stability and instability of the fractional-order system have been obtained, in terms of the main diagonal elements and the determinant of the system's matrix, as well as the fractional orders of the Caputo derivatives. Moreover, necessary and sufficient fractionalorder independent conditions have also been presented, in terms of the main diagonal elements and the determinant of the system's matrix, which guarantee the asymptotic stability or instability of the considered two-dimensional system, regardless of the choice of the fractional-orders considered in the system. These latter results prove to be especially useful in practical applications where the exact fractional orders of the Caputo derivatives are not precisely known.
It is important to note that multi-term fractional-order differential equations [18] and their qualitative properties are sharply linked to multi-order systems of fractional differential equations. We refer to [11] for a thorough presentation of the relationship between these two concepts. The investigation of stability properties of multi-term FDEs is so far limited to two-term and three-term fractional-order differential and difference equations, which have been recently studied in [19][20][21][22]. However, due to the increasing complexity of the problem, equations with four or more fractional terms have not yet been investigated. This paper is organized as follows: Section 2 illustrates the statement of the problem and the main definitions. Due to the importance in the description of the solution of linear systems of FDEs, in Section 3, we provide a detailed description of the Mittag-Leffler function, of its derivatives and of the corresponding asymptotic behavior. Section 4 investigates the stability properties of single-order systems of FDEs, by presenting classical Matignon's theorem and some simulations illustrating the different stability behavior in dependence of the spectral properties of the matrix system. Stability analysis of multi-order systems is discussed in Section 5; since general results are far from being formulated in this case, we focus on some special cases and we separately investigate two-dimensional systems, higher dimensional systems with block-triangular structure, and higher dimensional systems with some special fractional orders. Some concluding remarks are hence provided in the concluding Section 6.

Preliminaries
Consider an n-dimensional fractional-order system with Caputo derivatives: with q = (q 1 , q 2 , ..., q n ) ∈ (0, 1] n , assuming that f : [0, ∞) × R n → R n is a continuous function on its domain of definition, Lipschitz-continuous with respect to the second variable and y : [0, ∞) → R n a vector-valued function. With C D q y(t), we denote the application of the Caputo derivative of order 0 < q i ≤ 1 to each component y i (t) of y(t), namely Existence and uniqueness of the solution of initial value problems associated with system (1) is ensured by Corollary 2.4 from [14].
Whenever q 1 = q 2 = . . . = q n , system (1) is said to be single-order; otherwise, the term multi-order will be used.

Definition 1.
Let α > 0 and denote by ϕ(t, y 0 ) the unique solution of (1) satisfying the initial condition y(0) = y 0 ∈ R n . Then: i. the trivial solution of (1) is called stable if for any ε > 0 there exists δ = δ(ε) > 0 such that, for every y 0 ∈ R n satisfying y 0 < δ, we have ϕ(t, y 0 ) ≤ ε for any t ≥ 0; ii. the trivial solution of (1) is called asymptotically stable if it is stable and there exists ρ > 0 such that lim t→∞ ϕ(t, y 0 ) = 0 for y 0 < ρ; iii. the trivial solution of (1) is called O(t −α )-asymptotically stable if it is stable and there exists ρ > 0 such that, for any y 0 < ρ, we have:

Remark 1.
In the particular case of linear systems of fractional-order differential equations with constant coefficients, we say that the system is stable/asymptotically stable/unstable if and only if its trivial solution is stable/asymptotically stable/unstable.
For the purposes of this paper (the reasons will be clearer later on), it is convenient to study and introduce a further generalization of the ML function.

The Prabhakar Function and Its Asymptotic Properties
For three real parameters α, β and γ, the three-parameter Mittag-Leffler (ML) function, also known as the Prabhakar function [24], is defined by its series representation This function is not only a generalization, to three parameters, of the two-parameter ML function E α,β (z) (indeed, when γ = 1, it is E 1 α,β (z) = E α,β (z)), but it also provides a simple and elegant way to represent derivatives of two-parameter ML functions since as one can easily check after a term-by-term differentiation of (2). In order to introduce a result about the Laplace transform (LT), it is necessary to introduce what is known as the Prabhakar kernel for which the following analytical representation of the LT is available: Having in mind the stability analysis of linear FDEs, whose solutions will be expressed in terms of Mittag-Leffler functions and their derivatives, it is of interest to recall some results about the asymptotic behavior of the Prabhakar function in the complex plane.
In particular, for large arguments and 0 < α ≤ 1, we first identify exponential and algebraic expansions, respectively given by and, thanks to the results obtained by Paris [25,26], we know that with the sign in e ∓πi which must taken negative for z in the upper complex half-plane and positive otherwise. Following the convention adopted in [27], in each sum, we have first indicated the dominant term, namely the exponential term F γ α,β (z) when | arg z| < απ 2 and the algebraic term A γ α,β (ze ∓πi ) when απ 2 < | arg z| < απ. The lines | arg z| = απ and | arg z| = απ 2 are, respectively, Stokes and anti-Stokes lines where asymptotic expansions change their behavior. The above result is graphically summarized in Figure 1.  We first recall that coefficients c j in asymptotic expansion F γ α,β (z) are obtained [25] from the inverse factorial expansion, for |s| → ∞ in | arg(s)| ≤ π − and any arbitrarily small > 0, of with (x) j = Γ(x + j)/Γ(x) the Pochhammer symbol and ψ = 1 − γ + β. They can be evaluated by means of a sophisticated algorithm introduced in in [25] and also explained in [28]. The first few entries of c k are available in [26]. Based on the asymptotic properties of the Prabhakar function, we obtain the asymptotic equivalence: as t → ∞, where the sign in the term e ±γπi is positive if λ is in the upper complex half-plane, and negative otherwise.

Asymptotic Behavior of Derivatives of the ML Function
Thanks to the relationship (3) between derivatives of the ML function and the Prabhakar function, the investigation of the asymptotic behavior of any m-th order derivative of E α,β (z) is hence possible by applying the corresponding results for the Prabhakar function and afterwards replacing β with αm + β and γ with m + 1.
To this purpose, we first observe that, after these replacements, the function F Hence, coefficients c j vanish for j = m + 1, m + 2, . . . and the exponential and algebraic expansions read Therefore, by taking into account just the dominant expansions in each sector of the complex plane delimited by Stokes and anti-Stokes lines, and just leading terms in each expansion, we can describe the asymptotic behavior of derivatives of the ML function as Behavior of Derivatives of the ML Function When | arg z| = απ 2 It remains to investigate the behavior along the anti-Stokes line | arg z| = απ 2 where both the exponential and the algebraic terms are present. We therefore consider Clearly, the second term asymptotically goes to zero when ρ → ∞. The first term, instead, in modulus asymptotically tends to zero only for suitable values of α and β such that 1 − αm − β + m − j ≤ 0 for any j ∈ {0, 1, . . . , m}, namely, when 1 − αm − β + m ≤ 0 or, equivalently, when When we consider the one-parameter ML function E α (z), namely β = 1, which is the instance of the ML involved in the stability analysis of linear FDEs, for arg z = ± απ 2 just E α (z) asymptotically converges to 1/α for |z| → ∞ but any m-th order derivative of This situation is illustrated in Figure 2, where we report the first derivatives of E α (z) for α = 0.6 and α = 0.8 evaluated for z along the anti-Stokes line arg z = απ 2 (results are similar when arg z = − απ 2 ). Modulus of E α (z) and its first and second derivatives with arg z = απ 2 and α = 0.6 (left plot) and α = 0.8 (right plot).

Remark 2.
The behavior on the anti-Stokes line | arg z| = απ 2 of m-th order derivatives of the ML function E α (z), which are unbounded as |z| → ∞ for m ≥ 1, is quite different from that of the exponential function e z (namely, the special instance of E α (z) for α = 1). Indeed, derivatives of the the exponential are never unbounded on the corresponding anti-Stokes lines | arg z| = π 2 since there it is d m /dz m e z = 1 for any m = 0, 1, . . .

Stability of Linear Systems of Single-Order FDEs
We first consider the following linear system of Caputo-type fractional-order differential equations of the same fractional order: where q ∈ (0, 1] and A ∈ R n×n , coupled with the initial condition y(0) = y 0 ∈ R n . In is important to emphasize that system (6) is equivalent to the following system of weakly singular Volterra integral equations of convolution type (see, for example, [29,30]): For the most important advances regarding the general theory of linear Volterra integral equations, including the case when the convolution kernel is completely monotonic, we refer to [31][32][33][34].
The characteristic equation associated with system (6) is where, according to [35], the principal value (first branch) of the complex power function is considered. Therefore, it is easy to see that s is a root of the characteristic Equation (8) if and only if there exists an eigenvalue λ of the matrix A such that Hence, this leads to the following characterization of the stability properties of system (6), in terms of the roots of its characteristic equation: With the aim of investigating the stability properties of system (6) by characterizing the stability region S q , and presenting a concise proof of Matignon's theorem [6], it is convenient to use the Jordan normal form of the matrix A. Indeed, let us consider a nonsingular matrix P ∈ C n×n such that and λ k are eigenvalues of the matrix A. The size of the largest Jordan block J k of A associated with the eigenvalue λ k is called the index of λ k [36]. On the other hand, the total number of Jordan blocks associated with a given eigenvalue λ k in the Jordan normal form of the matrix A is the geometric multiplicity of the eigenvalue λ k . Moreover, the sum of the sizes of all Jordan blocks corresponding to λ k is the algebraic multiplicity of λ k . Therefore, the index of an eigenvalue λ k is equal to 1 if and only if its algebraic and geometric multiplicities are equal. With these observations, we next give a slightly modified version of the classical result of Matignon, to fix a small imprecision in the second statement, related to the use of the geometric multiplicity instead of the index of an eigenvalue: Theorem 1 (Matignon, 1996 [6]). The linear system (6) ii. stable if and only if σ(A) ⊂ S q and the eigenvalues of A which satisfy | arg(λ)| = qπ 2 have index 1.
Proof. With the notations introduced previously, denoting z(t) = Py(t), it is easy to verify that system (6) is equivalent to Applying the LT to the linear system (10) leads to the following formula for the LT of the vector function z(t): Since the Jordan normal form J is a block diagonal matrix, the matrix (s q I − J) −1 is also block diagonal, and its blocks are upper triangular matrices of the form: where d k represents the dimension of the k-th Jordan block J k . Correspondingly, the Laplace transform Z(s) is made up of "blocks" (of size d k ) of the form Applying the inverse LT, and taking into account that we obtain: Based on (5), we obtain the following asymptotic equivalence: Therefore, the following conclusions can be drawn: • e m q,(m−1)q+1 (t; λ) converges to 0 as t → ∞, if and only if | arg(λ)| > qπ 2 ; moreover, in this case, e m q, the function e m q,(m−1)q+1 (t; λ) is bounded if and only if m = 1. With the above observations, the conclusions of Matignon's theorem readily follow. We emphasize that, for the case of statement ii., if there exists an eigenvalue of A which satisfies | arg(λ)| = qπ 2 , the solutions of (6) are bounded if and only if the size of the largest Jordan block associated with this critical eigenvalue is equal to 1, i.e., the index of the eigenvalue is 1.

Remark 3.
The above proof slightly differs from the one in [6]. Matignon's proof, indeed, makes use of derivatives of the ML function instead of the Prabhakar kernel e γ α,β (t; λ) as in the proof of Theorem 1. A link between the two proofs can be, however, easily established in view of the relationship (3) between derivatives of the ML function and the Prabhakar function.

Remark 4.
Matignon's theorem implies that, if 0 < q 1 < q 2 ≤ 1 and system (6) is asymptotically stable for q = q 2 , then it will be asymptotically stable for q = q 1 as well. In particular, if the classical integer-order systemẏ = Ay is asymptotically stable (i.e., all eigenvalues of A have negative real part), it follows that the fractional-order system (6) is asymptotically stable, for any fractional-order q ∈ (0, 1).

Example 1.
To present numerical evidences of the above results, we consider here the linear systems of FDEs (6), with fractional order q = 2/3 and the coefficient matrix A chosen from one of the following four matrices: The solution y(t) = E q (t q A)y 0 evaluates by direct computation the matrix ML function thanks to the algorithm described in [37], after using the initial condition y 0 = 1, −4, −2, 4 T .
The value = 0.1 is used in A 3 and A 4 . The asymptotic behavior of the solution y(t) depends on the spectral properties of the matrix. In particular, we observe that: A 1 has two eigenvalues λ 1/2 = e ±q π 2 i laying on the border of the stability sector S q and both having index 2; according to Theorem 1, the system produces unbounded solutions as clearly shown in the left plot of Figure 3; A 2 has the same two eigenvalues λ 1/2 = e ±q π 2 i of A 1 , laying on the border of the stability sector S q , but their index is now 1; the expected bounded solutions are shown in the right plot of Figure 3; A 3 has two eigenvalues λ 1/2 with index 2, as A 1 , but now they lay inside the stability sector S q ; the asymptotically stable solutions are illustrated in the left plot of Figure 4; A 4 has two eigenvalues λ 1/2 with index 1, as A 2 , but lying outside the stability sector S q ; the resulting unbounded solutions are illustrated in the right plot of Figure 4.

Stability of Linear Multi-Order Systems of FDEs
Extending Matignon's theorem to the case of systems of FDEs with multiple fractional orders raises several technical difficulties, and, consequently, with the current state of the art, we are unable to present an exhaustive theory regarding this matter.
One of the technical difficulties that should be mentioned in this context is the fact that, for a general multi-order system of the form C D q y(t) = Ay(t), where A ∈ R n×n , q = (q 1 , q 2 , . . . , q n ) ∈ (0, 1] n (such that not all q i are equal), considering the Jordan normal form J of the matrix A and a nonsingular matrix P ∈ C n×n such that A = PJP −1 (similarly as in the previous section), the transformation z(t) = Py(t) does not lead to an equivalent system of the form C D q z(t) = Jz(t).
Therefore, different theoretical approaches should be used to tackle linear multi-order systems of FDEs.
Using another approach, namely the Laplace transform method, we first obtain the following system: where Y i (s) is the Laplace transform of the i-th component y i (t) of the solution y(t). System (13) is equivalent to the following system: where b i (s) = s q i −1 y i (0), for any i = 1, n and ∆(s) = diag(s q 1 , s q 2 , . . . , s q n ) − A.
Using standard properties of the Laplace transform [8,14,35], the following result holds: Theorem 2. The multi-order system (12) is asymptotically stable if all the roots of the characteristic equation det ∆(s) = 0 have negative real parts.
It is important to point out that, for large-scale systems with many different fractional orders for the Caputo derivatives, the analysis of the roots of the characteristic Equation (14) is a very difficult and complex task.
Nevertheless, the case of two-dimensional linear multi-order systems has been fully analyzed in [17], and a summary of the main results will be presented in the next section.

Stability of Two-Dimensional Systems of FDEs with Different Fractional Orders
In the general case of a two-dimensional linear system of fractional-order differential equations: C D q 1 y 1 (t) = a 11 y 1 (t) + a 12 y 2 (t) C D q 2 y 2 (t) = a 21 y 1 (t) + a 22 y 2 (t) (15) where A = (a ij ) ∈ R 2×2 and q 1 , q 2 ∈ (0, 1], applying the LT leads to the following characteristic equation: which can be written as where s q 1 and s q 2 represent the principal values (first branches) of the corresponding complex power functions [35]. Employing asymptotic properties and the Final Value Theorem of the LT [12,35], the following result [16] holds: Proposition 2.

1.
System (15)  In general, computing the roots of the characteristic Equation (16) is not a straightforward task. Thus, departing from Proposition 2, we seek to obtain necessary and sufficient conditions involving the coefficients a 11 and a 22 of the main diagonal of the matrix A as well as the determinant det(A), which guarantee the stability or instability of system (15).
We first concentrate our attention on fractional-order-dependent stability and instability conditions, as described below. The proof of the following results is rather elaborate, involving the root locus method, and has been presented in detail in [17]. Note that only the case det(A) > 0 is discussed here, as det(A) < 0 implies that system (15) is unstable, for any fractional orders (q 1 , q 2 ) ∈ (0, 1] 2 (in fact, it is trivial to show that, if det(A) < 0, the characteristic Equation (16) has at least one positive real root).
Theorem 3 provides a relatively simple algebraic criterion (in the form of inequalities comprising the elements of the main diagonal of the system's matrix A as well as its determinant and the fractional orders) that enables us to immediately decide the question of asymptotic stability or instability for a given two-dimensional multi-order system of fractional differential equations. In fact, Theorem 3 may be seen as a generalization of the Routh-Hurwitz stability criterion.

Remark 5.
If q 1 = q 2 := q, the curve Γ(δ, q 1 , q 2 ) reduces to the straight line: Therefore, Theorem 3 provides that, for equal fractional orders, system (15) is asymptotically stable if and only if The eigenvalues of the system's matrix A are and, hence, inequality (17) is equivalent to | arg λ 1,2 | > qπ 2 . Consequently, for two-dimensional systems, the conclusion of Matignon's theorem is recovered as a particular case of Theorem 3.
The next step is to seek necessary and sufficient conditions which ensure the asymptotic stability or instability of system (15) for any choice of the fractional orders. A complete investigation of the family of curves Γ(δ, q 1 , q 2 ) leads to the following fractional-order independent stability and instability results [17]: Theorem 4 (Fractional-order independent instability results). i.
If det(A) < 0, system (15) is unstable, regardless of the fractional orders q 1 and q 2 . ii. If det(A) > 0, system (15) is unstable regardless of the fractional orders q 1 and q 2 if and only if one of the following conditions holds: a 11 + a 22 ≥ det(A) + 1 or a 11 > 0, a 22 > 0, a 11 a 22 ≥ det(A).
The previous theorems provide easily verifiable necessary and sufficient conditions which ensure the asymptotic stability or instability of the two-dimensional system (15), for any choice of the fractional orders q 1 , q 2 ∈ (0, 1]. These conditions are expressed as simple inequalities involving the main diagonal elements a 11 and a 22 as well as the determinant det(A) of the system's matrix. On one hand, if det(A) < 0, Theorem 4 provides that system (15) is unstable, for any choice of the fractional orders q 1 , q 2 ∈ (0, 1]. Hence, we will focus our attention on the case δ = det(A) > 0. Let us denote by R s and by R u the fractional-order independent stability and instability regions given by Theorems 4 and 5: R u ={(a 11 , a 22 , δ) ∈ R 2 × (0, ∞) : a 11 + a 22 ≥ δ + 1 or a 11 > 0, a 22 > 0, a 11 a 22 ≥ δ} R s ={(a 11 , a 22 , δ) ∈ R 2 × (0, ∞) : a 11 + a 22 < 0 and max{a 11 , a 22 } < min{1, δ}} The regions R u and R s are plotted in Figure 8. The intersections of these regions with the δ = det(A) = 6 plane are shown in Figure 9. Moreover, it can be verified [17] that the union of all the curves Γ(δ, q 1 , q 2 ) (for δ > 0 and q 1 , q 2 ∈ (0, 1]) represents the complementary of R s ∪ R u (see Figure 9).  . Curves Γ(δ, q 1 , q 2 ) given by Lemma 1, for det(A) = δ = 6 and q i ∈ k 40 , k = 1, 40 , i = 1, 2 (1600 curves), color-coded from red to violet according to increasing values of q 1 q 2 . The red/blue shaded regions represent the intersections of the fractional-order independent stability and instability regions (see Figure 8) with the det(A) = 6 plane.

Remark 7.
The classical Routh-Hurwitz stability test for two-dimensional systems of the forṁ y = Ay provide that the system is asymptotically stable if and only if Tr(A) < 0 and det(A) > 0. However, based on Theorem 5, the additional inequality max{a 11 , a 22 } < min{1, det(A)} has to be verified in order to ensure the asymptotic stability of the fractional-order system (15), regardless of the choice of fractional orders q 1 and q 2 .
In conclusion, based on the previously described results, the following steps should be undertaken for the stability analysis of a two-dimensional system of FDEs: if det(A) < 0, then the system is unstable, for any choice of the fractional orders q 1 , q 2 ∈ (0, 1], based on Theorem 4; 2.
if (a 11 , a 22 , det(A)) / ∈ R s ∪ R u , then the stability properties of the system depend on choice of the fractional orders q 1 , q 2 ∈ (0, 1] and Theorem 3 should be applied. The results described in this section, particularly Theorems 3-5, give a comprehensive method to assess stability properties of two-dimensional fractional-order systems. However, the generalization of these results to higher dimensional fractional-order systems still remains an open question.

Stability of Higher Dimensional Systems of FDEs with Specific Structures
Consider that the matrix A of the linear system (12) has a block-triangular structure: We also assume that q = (q 1 , q 1 , . . . , q 1 d 1 times , . . . , q m , q m , . . . , q m d m times , q 1 m+1 , q 2 m+1 , . . . q 1 p , q 2 p ) ∈ (0, 1] n . In this case, the characteristic equation associated with system (12) is Therefore, combining Matignon's theorem (Theorem 1) and Theorem 3, the following statements are obtained: • system (12) is asymptotically stable if and only if for any i = 1, m and - , for any i = m + 1, p, where A 11 ii and A 22 ii are the main diagonal elements of matrix A ii , δ i = det(A ii ) and φ is defined in Lemma 1.
• system (12) is unstable if at least one of the following holds: there exists i ∈ {1, 2, . . . m} such that the matrix A ii has at least one eigenvalue λ such that | arg(λ)| < q i π 2 or -there exists i ∈ {m + 1, . . . , p} such that ii are the main diagonal elements of matrix A ii , δ i = det(A ii ) and φ is defined in Lemma 1.

Stability of Higher Dimensional Systems of FDEs with Special Fractional Orders
Let us consider the following n-dimensional linear multi-order system of fractional differential equations: where q i ∈ (0, 1], a ij ∈ R and n ≥ 3.
If the coefficient matrix of the system is not of a triangular or block triangular form as considered in the previous section, one can not provide a comprehensive stability theory. Still, an approach that works under certain restrictions on the fractional orders of the Caputo derivatives has been developed in [14]. We will next recall the general results obtained by the mentioned authors.
Suppose q j ∈ (0, 1], for any j = 1, n and that there exists q * ∈ (0, 1] and ρ j ∈ Q such that q j = ρ j q * . It follows that there exists r j , s j ∈ N for j = 1, n such that gcd(r j , s j ) = 1 and ρ j = r j s j . Let s be the least common multiple of the denominators s j . Then, for any j, there exists α j ∈ N such that We can rewrite the j-th equation of system (20) as an equivalent system of α j differential equations having the order q * s . It follows that system (20) can be expressed as a system of n * = n ∑ j=1 α j equations of order q * s : where A * has the following block structure Even though the dimension n * of the system may be significantly higher than the dimension n of the original system (20), resulting in higher computational costs, all the equations of the new system (21) now have the same fractional order, giving an advantage in studying the stability properties of the solutions of the system.
We expose the main result of this section, based on [14], which gives us stability criteria involving the components of the matrix A * . Theorem 6. Suppose that q j ∈ (0, 1] for any j and there exists q * ∈ (0, 1] and ρ j ∈ Q such that q j = ρ j q * , for all j. Then, all the solutions of system (20) converge to zero at infinity if the eigenvalues λ * j of the associated system's coefficient matrix A * satisfy |arg λ * j | > πq * 2s , ∀j = 1, n, with s being the least common multiple of the denominators of ρ j .
Example 4. Again, we reconsider system (18) with q 1 = 0.2 and q 2 = 1. In this case, the matrix A * given by the above procedure is  (18) is equivalent to a system of six fractional-order differential equations with the same order q = q 1 = 0.2. The matrix A * has a pair of complex conjugated eigenvalues (λ, λ), λ = 0.543842 + i 0.133131 such that | arg(λ)| = 0.240076 < 0.1π = qπ 2 . Hence, based on Matignon's theorem, system (18) is unstable for q 1 = 0.2 and q 2 = 1. Therefore, this is in good agreement with the results obtained in Example 2, based on Theorem 3.
However, it is important to note that cases q 1 = 1 π and q 2 = 1 cannot be investigated using the technique provided by Theorem 6.

Conclusions
An extensive analysis of stability properties of linear systems of FDEs has been provided. This analysis is of importance to describe the asymptotic behavior of physical systems when modeled by means of FDEs. Both single-order and multi-order systems have been studied, reviewing the most important theoretical results that have been obtained so far in the literature. The role of the Mittag-Leffler function, and of its derivatives, has been highlighted and a presentation of their asymptotic behavior has been proposed. We have seen that, unlike systems of integer order, coefficients of the systems are not sufficient to describe stability properties of solutions, due to the tight dependence on the order of the fractional derivatives. This dependence becomes more and more difficult to investigate in systems incorporating derivatives of different order, as we have observed from the analysis of two-dimensional systems. Stability analysis of multi-order higher dimensional systems is still an open problem which deserves to be investigated with more attention; with this work, a first contribution has been provided by examining systems with some specific structures, and we hope these results will stimulate the analysis of more general systems.