Stability Analysis of Fractional Differential Equations with Unknown Parameters

In this paper, the stability of fractional differential equations (FDEs) with unknown parameters is studied. FDEs bring many advantages to model the physical systems in the nature or man-made systems in the industry. Because this representation has a property between linear differential equations and nonlinear differential equations. Therefore, the designer may use the FDEs to model complex systems instead of nonlinear differential equations which have hard mathematical background. Using the graphical based D-decomposition method, we investigate the parametric stability analysis of FDEs without complicated mathematical analysis. To achieve this, stability boundaries are obtained firstly, and then the stability region set depending on the unknown parameters is found. The applicability of the presented method is shown considering some benchmark equations which are often used to verify the results of a new method. Simulation examples shown that the method is simple and give reliable stability results.


INTRODUCTION
Fractional differential equations [1,[2][3][4][5] are a generalization of classical integer order differential equations through the application of fractional calculus [6,[7][8][9][10] which has been developed by pure mathematicians firstly since after half of the 19th century though engineers and physicist found applications of fractional calculus for various concepts 100 years later [11]. As a field of mathematical analysis, fractional calculus studies the possibility of taking real or complex number powers of differential operators. It may be considered an old branch of mathematical analysis, but it is a novel topic yet [12]. Especially, fractional calculus has gained a great deal of popularity in modelling some physical and engineering systems as well as fractal phenomena in the last few decades [13,[14][15][16][17][18][19]. In fact, many systems in the real world are now better characterized by FDEs and analysed by numerical techniques developed for solving differential equations involving non-integer order derivatives [20]. FDEs are also known as extraordinary differential equations.
Continuing technological developments have required new methods in basic sciences, especially in mathematics for analysis and design of physical systems and their control tools. These methods which are easily implemented with the advancement of high speed computers aimed to better and better characterization, design tools and control performance of modern technological products of engineering systems of developing civilization. These development which had covered only static system models involving geometry and algebra until 1965, had started using dynamical models involving differential and integral calculus since 1965; and now have been accelerating since the 1960's with the fractional order modelling involving FDEs which have gained force and dare with high speed computers [21,22].
Hence FDEs have become a powerful tool in studying, designing and control of physical systems and engineering products of the present world and it still constitutes a research area in need.
Stability is one of the most important objects in the analysis and design of dynamical systems. If the differential equation of a system has not a stable property, the system may burn out, disintegrate or saturate when a signal is applied [23]. Therefore an unstable system is useless in practice and needs a stabilization process via an additional control element mostly [24,25]. If a system has unknown parameters, the stability analysis is called as parametric stability analysis. It is more difficult than the classical stability analysis which has simple analysis methods such as Routh-Hurwitz method, Nyquist stability theorem, etc. In the literature, there are some methods on the stability of FDEs with uncertain parameters [26,27,28]. The uncertainty in these studies is considered by a certain interval of the unknown parameters. However, to consider the unknown parameters' values in a whole parametric interval from zero to infinity is more useful than a certain little interval, which is the subject of this paper.
Motivated by the need of stability analysis for FDEs with unknown parameters, we suggest in this paper an efficient graphical based stability analysis using the D-decomposition method [29,30]. The Ddecomposition method provides a powerful and simple stability work environment to the analyst. This method is based on a conformal mapping from frequency domain to parametric domain of unknown parameters. With this mapping, the imaginary axis which is the stability bound of complex s-plane converts to three types of stability boundaries, which are named as real, complex and infinite root boundaries, in the parametric space. These boundaries give us the stability regions which are important tools including useful stability knowledge. The algorithm presented in this paper has a reliable result which is illustrated by several examples, and is practically useful in the computerized analysis of FDEs having unknown parameters.
The paper is organized as follows: The basic principles of FDEs are revisited in the next section. The stability concept for FDEs is explained in Section 3. A derivation of the stability boundary formulae for the stability regions of fractional differential equations with unknown parameters is given in Section 4. The next section illustrates the effectiveness of the stability analysis proposed with three simulation examples. Finally, Section 6 gives some concluding remarks.

PRINCIPLES OF FRACTIONAL DIFFERENTIAL EQUA-TIONS
In the most general case, FDEs are expressed by the following form [31]: where F and G are fractional differential functions, α i (i = 1 ∼ n) and β k (k = 1 ∼ m) are positive real numbers such that 0 < α 1 < α 2 < · · · < α n , 0 < β 1 < β 2 < · · · < β m and m < n. α D γ t is fractional order derivative and integral operator and it is defined as follows [32]: Here γ is fractional order, (γ) is the real part of fractional order and a is a constant coefficient related with initial conditions. Commonly, t is an independent variable representing time, u(t) is input exiting function and y(t) is the output response function of a dynamical system. There are various definitions for fractional derivative. Riemann-Liouville, Grünwald-Letnikov, Caputo and Mittag-Leffler are the well-known and common definitions among them. (For a more detail of these definitions, the reader can see [3,32]). One of the most common types of FDEs is linear time-invariant fractional differential equations n i=0 where a i and b k are real numbers. Eq. (3) is also called as non-commensurate order FDE. As a special case, fractional orders a i and b k may be multiple of same real number ε like a i = iε and b k = kε. In this case, Eq. (3) is named by commensurate order FDE [3,33,34]. The Laplace transform method is commonly used in engineering systems and their analysis. According to Grünwald-Letnikov definition, Laplace transform of differential operator a D γ t is given by where s is the Laplace operator. For the differential equation in Eq. (3), the transfer function giving the input-output expression of a system with zero initial conditions is given by In this equation, U (s) is the Laplace transform of the exciting function u(t); similarly Y (s) is that of the response y(t). N (s) and D(s) are the numerator and denominator polynomials of the transfer function, respectively. Being the stability as a first, the function given in Eq. (5) contains many important system information and concepts.

STABILITY ANALYSIS OF FRACTIONAL DIFFEREN-TIAL EQUATIONS
There are many ways of testing the stability of a linear time-invariant differential equation. Stability of the differential equation can be examined by applying Routh-Hurwitz test on its denominator polynomial, checking the locations of the poles of its transfer function whether they are on the left half s-plane or not and exploring its output whether the output remains bounded with an impulse or step input excitation [35]. If the differential equation involves time-delay or fractional order terms, in this case, Routh-Hurwitz' criteria cannot be applied.
The denominator seen in Eq. (5) of a FDE is in the form of a quasi-polynomial and it is expressed by D(s) = a n s an + a n−1 s an−1 + · · · + a 1 s a1 + a 0 .
For the stability analysis, the quasi-polynomial in Eq. (6) is transformed to the following commensurate order quasi-polynomial where α is the least common multiple of α 1 , α 2 , · · · , α n . Stability condition of this fractional order polynomial was given by Matignon [36] in 1996 as follows where −λ i (i = 1, 2, · · · , n) are the roots of the pseudo-polynomial P (s α ) [37,38]. Matignon's stability analysis is applicable for only linear FDEs whose coefficients are known and invariant. For the FDEs changing their coefficients/parameters in an interval, Chen et al. [27] proposed a very effective method for the stability analysis. However, to investigate the stability may be more useful in the case that the parameters of a FDE changes between minus and plus infinity. In this paper, a method is presented for this type of stability analysis. The method is based on obtaining stability boundaries and it contains a graphical presentation. The most important property of this method is to construct a conformal mapping from s plane to parameter space composed by unknown parameters of the FDE. Therefore, the method is called as parametric stability analysis and based on the D-decomposition method (see [39] for more detail).
In the D-decomposition method, there are three stability boundaries of a polynomial [39,40]. The first boundary belongs to a real pole which changes its stability property when passing through origin and crossing the opposite half of s−plane with the parameter changes. Therefore, this boundary is called real root boundary. It is obtained by putting zero instead of s in Eq. (6) with Second boundary is named by infinite root boundary because of the fact that the boundary belongs to a pole which changes stability property at infinity with the parameter changes. Infinite root boundary is determined by equalizing the coefficient of the greatest order term to zero; a n = 0.
In order to obtain the boundaries from Eqs. (9) and (10), the coefficients a 0 and a n should contain unknown parameters. Otherwise, these boundaries do not exist for the considered FDE.
The last boundary results a couple of complex poles passing to one half plane from the other half plane over the imaginary axis of the s−plane with the parameter changes. This is the main boundary which determines the stability region of the FDE and is named by complex root boundary. To obtain this boundary, s in Eq. (6) is replaced by jω as follows; D(jw) = a n (jω) an + a n−1 (jω) an−1 + · · · + a 1 (jω) a1 + a 0 .
Using the expansion j x = cos(0.5πx) + j sin(0.5πx) in Eq. (11) for the fractional order powers of the complex number j, we get By decomposing D(jw) into real and imaginary parts, we obtain By equalizing D t (ω) and D i (ω) to zero separately, two variable equations system whose variables are the unknown parameters of the FDE are obtained. The complex root boundary is found by solving of this system with respect to ω for 0 < ω < ∞. The stability region of the FDE is enclosed by these three boundaries in the parameter space. The main property of this region is that all parameters in this area make the FDE stable.

PARAMETRIC STABILITY ANALYSIS OF FRACTIONAL DIFFERENTIAL EQUATIONS
In this section, we consider the FDEs which are often encountered in the engineering systems have the following form: where u(t) and y(t) are forcing and response functions, respectively; a, b, c and k are real coefficients, α 1 and α 2 are the fractional order powers such that 0 < α 1 < α 2 < 2.
As it is pointed out in the previous section, the transfer function of a FDE encloses important stability information. By taking the Laplace transform of the FDE in Eq. (14), the transfer function giving input-output relation is obtained as If the coefficients a, b, c and the fractional orders α 1 and α 2 are constant or they change in a specific interval, the stability of Eq. (15) can be determined by Matignon's method [36] and Chen et al's method [27] easily as mentioned in Section 3. Here, we investigate the parametric stability for the full variation range of these coefficients using the D-decomposition method. For obtaining the stability boundaries, we consider the denominator of Eq. (15) D(s) = as α2 + bs α1 + c.
The real root boundary for Eq. (16) is determined as Since α 2 > α 1 , the infinite root boundary is found by applying Eq. (10) to Eq. (16) as the follows For the complex root boundary, the real and imaginary parts of Eq. (13) are equalized to zero and the following system of equations is obtained: aω α2 sin(0.5πα 2 ) + bω α1 sin(0.5πα 1 ) = 0, By solving this system of equations with respect to b, α 1 , α 2 for the coefficients a and c, we get the following formulas: By changing ω from zero to infinity for different values of b, α 1 and α 2 , the complex root boundaries are obtained. Definition 4.1: When three boundary conditions defined in Eqs. (17), (18), (21) and (22) are drawn in (a, c)−parameter plane together, the plane separates to many regions. The most basic property of these regions is that all points in any region have the same number of stable or unstable roots [39]. With this reference, the region whose all poles are stable are called as stability region. Every point in this region makes the FDE in Eq. (14) stable. The stability of each region can be determined by selecting a test point in the region and checking the stability of Eq. (16) in every region [36]. Plot of stability boundaries and stability regions on (a, c)−parameter plane will yield the parametric investigation of the stability of the given fractional order system. This provides for example to design a more robust system.

SIMULATION EXAMPLES
In this section, stability analysis of some fractional differential equations commonly used in the literature are given for illustration of the validity of the presented method. In the first example, fractional Basset equation defining the dynamical motion of an object which submerged into a fluid is considered. In the next example, as a more general case, stability analysis is investigated for a commensurate order FDE. In the last example, the stability ranges of the parameters a, b and c for an industrial heating furnace is investigated.

EXAMPLE 1
The motion dynamics of an object which is submerged into an incompressible fluid is one of the frequently studied topics in engineering literature. Basset [41] proposed an equation and its solution for a sphere moving in a viscous liquid when the sphere is moving in a straight line under the action of a constant force, such as gravity, and also when the sphere is surrounded by viscous liquid and is set in rotation about a fixed diameter and then left to itself [42]. This equation is named as Basset equation and it is expressed by the following fractional order differential equation where α ∈ (0, 1), and a = 0, b, c are arbitrary real coefficients [43]. While the equation in (23) is called by classical Basset equation for α = 0.5, it is named as generalized Basset equation for 0 < α < 1 [44].
For the stability analysis of Basset equation, Govindaraj and Balachandran proposed an analytical method in [45]. Even though they found solutions for the changes in the coefficients a,b and c separately, they have not presented a general solution with respect to changing in the coefficients. The goal in this example is to obtain a simple graphical result giving the stability or instability of the equation according to the changing of the parameters for the system given in Eq. (23).
Real and infinite root boundaries are the same with in Eqs. (17) and (18), respectively. For the complex root boundary, putting α 2 = 1 and 0 < α 1 = α < 1 in Eqs. (21) and (22), the expressions giving complex root boundary are found by For the stability analysis of Eq. For changing ω from 0 to ∞, these boundaries decompose the (a, c)−plane into many number of regions for various values of b. For example, the stability boundaries constitute five regions in the (a, c)−plane for the value of b = −2 as shown in Fig. 1 where complex root boundary is ac = b 2 2 = 2, a, c > 0 from Eq. (28). Since the regions are unlimited throughout the axes of a and c, the figure is limited for good visibility in the interval of [−10, 10] for these axes. The most important characteristic of these regions is that all points in every region have the same number of stable and unstable roots. Because of this reason, to determine which regions are stable or not among these five areas, it is sufficient to select only one testing point from every region and checking the stability of Eq. (23) according to these points. As shown in Fig. 2, it is found that the second and fourth regions are the stability regions. For verification of these regions, it can be seen that the results are suitable with the following results  which are given by Govindaraj and Balachandran [45] for the classical Basset equation. By varying b and repeating the above procedure, different stability regions are obtained for each b value as shown in Fig. 3. It is seen from this figure that smaller values of |b| provide bigger stability regions. To illustrate the graphical results more clearly, the overall stability region can then be visualized in a 3-D plot as shown in Fig. 4.

EXAMPLE 2
Commensurate order FDEs are commonly used for modelling of physical systems and industrial processes [46]. In this example, the following equation containing two fractional terms is considered: This equation represents a FDEs family for different values of α and it also contains classical Basset equation for α = 0.5. Members of this family are named by multi-term differential equations if the power of the greatest derivative term is greater than 1 (or α > 0.5) and single-term differential equations if the power of the greatest derivative term is less than 1 (or α < 0.5) [47].
For the various values of α, the stability regions can be easily obtained according to the values of b. For example, the stability regions for α = 0.2 and α = 0.8 are plotted for b ∈ [−5, −1] as shown in Fig.  8 where the complex root boundary calculated from (30) and (31) is ; a, c > 0.
Stability regions for the values of b changing in the interval [1,5] are replacements of the same  regions with respect to origin symmetrically. In Fig. 9, the stability regions for five different values of α are seen when b = −2 and b = 3. For b = −2, when the values of α are changed more often for instance, 100 times in the range of (0, 1), the change of stability region is appeared from Fig. 10. It is seen from this figure that the stability region fills right upper part of the (a, c)−plane if the value of α approaches to 0 and the stability region is getting smaller if the value of α

EXAMPLE 3
In this example, the following incommensurate order FDE is considered for modelling an industrial heating furnace. In this equation, nominal values of a, b and c parameters are given as a = 14994, b = 6009.5 and c = 1.69 in [48]. Sondhi and Hote [49] has shown the stability of the equation for these nominal values. The goal in this example is to verify this stability results and to determine the stability intervals by assuming these parameters are varying. Stability boundaries for Eq. (32) are as follows: Real root boundary: Infinite root boundary: a = 0, Complex root boundary: a = −1.1303bω −0.34 and c = −0.576bω 0.97 .
The stability region of the system for b = 6009.5 is shown in Fig. 11 where from (33c) the lower left stability region is found bounded by the curve

CONCLUSIONS
In this paper, a graphical based stability analysis method is presented for the FDEs. The analysis concept is based on the derivation of the stability boundaries and then the determination of the stability region including the parameter set which makes the FDE stable. One of the most important advantages of the method is that the stability analysis is done in a visual environment without considering complex analytical solutions. This method can be used not only for the investigation of the stability of a differential equation whose parameters are not changed but also for observation of the parametric robust stability of the equation whose parameters are varying in a large interval. From this aspect, having a large usage perspective of the method in comparison with the other methods is one of the other advantages of the method. This provides opportunity to engineers in their analyses about discussion more detail. Simulation examples have been selected from the benchmark problems encountered in engineering systems. As evidenced by the results given in these examples, it can be concluded that the proposed graphical based method is reliable method not only for stability analysis but also parametric robust stability analysis according to the parameter changes. The presented method can be generalized for stability analyses of fractional differential equations having time delay which is a very popular subject in the last decade. Moreover, the differential equations having more fractional terms can be also studied. Here, when the number of unknown parameters increases, the three dimensional graphs will be insufficient. In this case, more than one graphs or four dimensional graphs with the fourth dimension assuming by colour can be used.