Stability Analysis of Linear Time-Invariant Fractional Exponential ...
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: EXPRESS BRIEFS, VOL. 61, NO. 9, SEPTEMBER 2014
721
Stability Analysis of Linear Time-Invariant Fractional Exponential Delay Systems Mohammad Ali Pakzad and Mohammad Ali Nekoui
Abstract—This brief presents a new approach for the stability analysis of linear fractional exponential delay systems with commensurate orders and multiple commensurate delays that enable us to decide on some cases that were previously open problems. In the proposed approach, first an auxiliary polynomial is generated by mapping the principal sheet of the Riemann surface and a pseudodelay transformation. Next, this auxiliary polynomial is employed in determining all possible purely imaginary characteristic roots for any positive time delay. Then, the concept of stability is expressed as a function of delay. Two illustrative examples are provided to demonstrate the effectiveness of the proposed method and to gain a better understanding of the problem, and the root-locus curve of these systems has been plotted as a function of time delay. Index Terms—Fractional exponential systems, fractional-order systems, root locus, stability analysis, time-delay systems.
I. I NTRODUCTION
T
IME delay is an inherent part of many dynamic and physical systems. It is known that time delay could lead to inadequate performance and even instability of systems. Therefore, time-delay systems play significant roles in theoretical [1]–[7] and practical fields [8]–[10]. There are many linear infinite-dimensional systems with transfer functions that can be represented as fractional differential systems. When time delays and a fractional-order derivative appear in dynamic systems, we have fractional-delay systems. The characteristic equation of a fractional-delay system involves exponential-type transcen√ dental terms, i.e., e−τ s or e−τ s ; thus, a fractional-delay system has, in general, an infinite number of characteristic roots. In this brief, one of the most important and unresolved problems of fractional-delay systems is studied, i.e., the asymptotic stability of a general class of linear fractional exponential delay systems. Ozturk and Uraz [2], [4] are probably the pioneers in the stability investigation of fractional-order time-delay systems. They have developed the Ruth–Hurwitz criterion to evaluate the stability of a √ class of distributed systems whose transfer √ functions involve s and/or e−τ s . In [6], the necessary and sufficient conditions for the bounded-input–bounded-output Manuscript received January 17, 2014; revised April 30, 2014; accepted June 27, 2014. Date of publication July 8, 2014; date of current version September 1, 2014. This brief was recommended by Associate Editor J. Lu. M. A. Pakzad is with the Department of Electrical Engineering, Science and Research Branch, Islamic Azad University, Tehran 14778 93855, Iran (e-mail: [email protected]). M. A. Nekoui is with the Faculty of Electrical and Computer Engineering, Khajeh Nasir Toosi University of Technology, Tehran 19697 64499, Iran (e-mail: [email protected]). Color versions of one or more of the figures in this brief are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TCSII.2014.2335427
(BIBO) stability of fractional-order delay systems have been introduced. From the numerical analysis point of view, the effective numerical algorithms have been discussed in [11] for the evaluation of the stability of fractional-order delay systems. In [12], a heavy computation scheme based on Cauchy’s integral has been proposed to test the stability of such systems. Pakzad and Pakzad [13] have successfully extended the direct analytical method (presented by Walton and Marshal [14] for the stability analysis of integer-order time-delay systems) to fractionalorder delay systems. The stability of these systems can be also analytically analyzed by employing Hassard’s theorem proposed by Shi and Wang [15]. Based on the method they presented, the number of unstable poles of the characteristic equation for a given delay is determined. Recently, in [16], an analytical method for finding the stability regions of fractionaldelay systems with commensurate orders has been presented, which uses the bilinear Rekasius transformation to eliminate the exponential-type transcendental term in the characteristic equations of these systems. To gain more information about some of the recently developed methods and tools in this area, see [17] and [18]. A large number of circuits and systems have distributed parameters and/or delay elements. Thermal processes, the hole diffusion of transistors, electromagnetic devices, and transmission lines are typical examples. In some cases, the mathematical models in the Laplace-transform domain √ for the transfer function of a system may even consist of s in combination √ with fractional-delay term e−τ s . Many examples of fractional differential systems can √be found in literature. Simple examples, such as G(s) = e−τ s /s, with τ > 0 arising in the theory of transmission √ lines, √ are given in [19]. As another example, G(s) = tanh( s)/ s appears in a boundary-controlled and observed diffusion process in a bounded√domain In addi√ [20]. √ tion, transfer functions G(s) = cosh(x s)/ s sinh( s) and 0 √ √ G(s) = sinh(x0 s)/ sinh( s) appear in the boundary control of a 1-D heat equation with Neumann–Dirichlet √ √ boundary con√ ditions [20], or G(s) = 2e−τ s /b s(1−e−2τ s ) can be found in [21]. Other examples of this type can be found in [20]–[23]. The branch point of the transfer function of a system may be at a point rather than at the origin. As an example of this type, consider the asymmetric transmission line shown in Fig. 1. In this figure, R, C, L, and G show the resistance, capacitance, inductance, and conductance, respectively, in the unit of length, and the length of the line is assumed equal to L. According to Krishoff’s law, the equations governing this system are v(x, t) di(x, t) = Gv(x, t) + C (1) − dt dt i(x, t) dv(x, t) = Ri(x, t) + L (2) − dt dt
1549-7747 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
722
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: EXPRESS BRIEFS, VOL. 61, NO. 9, SEPTEMBER 2014
II. M ETHODOLOGY
Fig. 1. Model of an asymmetric transmission line [24]. Transfer function Vo (s)/Vs (s) has branch points at R/L and −G/C.
subject to the following boundary conditions: v(0, t) = vs (t)
i(L, t) = 0.
(3)
Taking the Laplace transform from both sides of (1)–(3) (assuming that the initial condition of the system is equal to zero) and after algebraic calculations, the transfer function of the system is obtained as follows [24]: G(s) =
Vo (s) 1 . = 1 1 Vs (s) cosh L(R + Ls) 2 (G + Cs) 2
(4)
In literature, the fractional-delay systems are those described by the following transfer function: 1 r q (s)e−βk s q0 (s) + nk=1 n2 k (5) G(s) = p0 (s) + k=1 pk (s)e−γk sr where r is a real number such that 0 < r < 1, pk are of the form lk mk αi δi i=1 qk ai s , and qk are of the form i=1 qk bi s . The fractional-delay characteristic equation is obtained by equating the denominator of a fractional-delay transfer function to zero. More formally, the fractional-delay characteristic equation under consideration in this brief can be expressed in the general form of CE(sα , τ ) = p0 (sα ) +
n
P (sα )e−τ
The main objective of this section is to present a new method for the evaluation of stability and the determination of the unstable roots of a fractional-order delay system. In this approach, we are trying to calculate the values of τ where there exists a crossing of poles through the imaginary axis. With these values in hand, we will be able to calculate the direction of crossing from the left half-plane to the right half-plane, which we will denote as a destabilizing crossing, or from the right to the left, which is a stabilizing crossing.
√
s
(6)
=1
where 0 < α < 1. We find out from the work in [22] that the transfer function of a system with a characteristic equation in the form of (6) will be BIBO stable if and only if it does not have any pole at (s) ≥ 0 (in particular, no poles of fractional order at s = 0). For fractional-order systems, if an auxiliary variable of v = sα is used, a practical test for the evaluation of stability can be obtained. In this new variable, the instability region of the original system is not given by the right half-plane but in fact by the region described as απ |∠v| ≤ (7) 2
A. Imaginary Axis Crossings In order to find the points on the imaginary axis where the crossing takes place, the auxiliary variable described in the previous section and a new substitution are used. As a result, the characteristic equation is converted into a solvable equation. To simplify this characteristic equation, we propose a new transformation that is defined as √ 1 − jT e−τ s , T ∈R (10) =r s=jω 1 + jT where r = exp(−2 tan−1 T ). It is important to note that this √ substitution is an exact expression of e−τ s for purely imaginary roots s = ±jω. Note that, for s = ±jω, the magnitudes of both sides in (10) always agree with each other. Furthermore, for the transformation to exactly hold, it is required that the phase condition developed from (10) holds. This condition can be found as √ 2 2 (11) τ = √ tan−1 (T ). ω which is the inverse mapping from the (T, ω) domain to τ . Note that there is a one-to-one mapping between T , ω, and τ . To summarize, the transformation in (10) is indeed exact at s = jω as long as (11) holds. By inserting (10) into (6), we have n 1 − jT CE(sα , T ) = p0 (sα ) + P (sα )r . (12) 1 + jT =1
By applying this auxiliary variable v = sα in (12), the following relation is obtained: n 1 − jT P (v)r . (13) CEv (v, T ) = p0 (v) + 1 + jT =1
with v ∈ C. Note that, under this transformation, the imaginary axis in the s-domain is mapped into απ (8) ∠v = ± . 2
This will transform the domain of the system from a multisheeted Riemann surface into a complex plane, in which the poles can be easily calculated. By multiplying (13) by (1 + jT )n , the new form of the characteristic equation is reached as
Let us assume that s = ±jω or, in other words, s = ωe±jπ/2 are the roots of characteristic (6) for τ ∈ R+ . Then, for the auxiliary variable, the roots are defined as απ απ ± j sin . (9) v = sα = ω α e±jαπ/2 = ω α cos 2 2
h(v, T ) = (1 + jT )n Cv (v, T ) = (1 + jT )n p0 (v) n P (v)r (1 − jT ) (1 + jT )(n−) . (14) +
Therefore, with auxiliary variable v = sα , there is a direct relation between the roots on the imaginary axis for the s-domain with the roots having argument ±απ/2 in the v-domain.
This expression is a polynomial in v of which the coefficients are parametric functions of T and r. When s = jω are the roots of (6) for some ω and τ ∈ R, then v = ω α ejαπ/2 will be also
=0
PAKZAD AND NEKOUI: ANALYSIS OF LINEAR TIME-INVARIANT FRACTIONAL EXPONENTIAL DELAY SYSTEMS
the roots of (14) for some ω and T ∈ R. To find the crossing frequencies from the imaginary axis in (6), v = ω α ejαπ/2 should be inserted into relation (14) and then the real and imaginary parts of the resulting equation should be separated as h(v, T ))|
jαπ v=ω α e 2
= h (ω α , T ) + jh (ω α , T ) = 0. (15)
In the given relations, h = (h) and h = (h). By replacing ω α with w in the given equation, we have h(w, T ) = h (w, T ) + jh (w, T ) = 0.
(16)
Making the real and imaginary parts of (16) equal to zero results in two algebraic equations; let us investigate those (ω, T ) solutions from h = h =
n i=0 n
ai (T )wi = 0
(17)
bi (T )wi = 0.
(18)
i=0
We utilize the resultant theory (see [25]) to eliminate w from the two multivariate polynomials h and h . Definition 1: Consider the two multivariate polynomials (17) and (18) in terms of w and T with real coefficients, where h and h have positive degrees in terms of w, and n > 0. The result of h and h with respect to w is defined by an (T ) an−1 (T ) . . . . . . . . . 0 .. .. .. .. .. .. . . . . . . .. ... . . . . a1 (T ) a0 (T ) 0 Rw (h , h ) = (19) 0 bn (T ) bn−1 (T ) . . . . . . . . . . .. .. .. .. .. . . . . . . . .. 0 ... . . . . b1 (T ) b0 (T ) which is the determinant of the well-known Sylvester matrix [25]. A 2n-order Sylvester matrix is constructed via (19) and its determinant Rw (h , h ) is a function of T as D(T ) = Rw (h , h ) = 0.
(20)
Hence, the set of all w and T that makes both equations zero can be obtained. Then, for every ω = w1/α and T , its corresponding τ can be determined using (11). The whole ω value, for which s = jω is a root of (6) for some nonnegative delays, is defined as the following crossing frequency set: 1
Ω = ω ∈ R+ C s α , τ = 0, for some τ ∈ R+ . (21) Corollary 1: If the system given as (12) is stable for τ = 0 (i.e., a system without delay) and Ω = φ, then the system will be stable for all positive values of τ . Proof 1: From the fact that there are no roots crossing the imaginary axis. This class of systems only exhibits a finite number of possible imaginary characteristic roots for any τ ∈ R+ at the given frequencies. Moreover, this method detected all of them. Let us call this set {ωc } = {ωc1 , ωc2 , ωc3 , . . . , ωcn }
(22)
723
where subscript c refers to “crossing” the imaginary axis. Furthermore, to each ωcm , m = 1, 2, . . . , n finitely correspond to the following τ values: {τc } = {τc1 , τc2 , τc3 , . . . , τcn }.
(23)
B. Direction of Crossing After the crossing points of characteristic (6) from the imaginary axis are obtained, the goal now is to determine whether each of these root crossings from the imaginary axis is a stabilizing cross or a destabilizing cross. Assume that (s, τ ) is a simple root of CE(sα , τ ) = 0. The root sensitivities associated with each purely imaginary characteristic root jω with respect to τ is defined as ∂CE ds s ∂τ Sτ |s=jωc = = − ∂CE . (24) dτ s=jωc ∂s s=jωc
The root tendency for each ωcm and τmk is defined as Root Tendency = RT|s=jωc = sign Sτs | s=jωcm . (25) τ =τmk
If it is positive, then it is a destabilizing crossing, whereas if it is negative, it is a stabilizing crossing. III. I LLUSTRATIVE E XAMPLES We present two cases, which display all the features discussed in this brief. Example 1: Consider the following linear time-invariant fractional-delay system: √
0.3e−τ s √ . G(s) = √ 6 √ 4 √ ( s) + 0.2( s) + ( s)2 + 0.3e−τ s The corresponding characteristic equation of G(s) is √ √ √ √ CE(s, τ ) = ( s)6 + 0.2( s)4 + ( s)2 + 0.3e−τ s . (26) This example has been taken from the work in [4]. The characteristic in (26) has poles in s = 0.0461 ± j1.0123 and s = −0.2921 for τ = 0. Therefore, this system is unstable for τ = 0. Our objective in this example is to find stability windows based on the method described in this brief for this system. By applying the criterion expressed in the previous section and by using auxiliary variable v = s0.5 , we can convert (26), according to (14), as h(v, τ ) = (v 6 + v 4 + v 2 )(1 − jT ) + 0.3r(1 − jT ). (27) √ √ √ By√ inserting expression v = ω ej(π/4) = ω(( 2/2) + j( 2/2)) into (27) and equating the real and imaginary parts of the obtained relation to zero, we get h = T ω 3 − 0.2ω 2 − T ω + 0.3r = 0 h = −ω 3 − 0.2T ω 2 + ω − 0.3T r = 0.
(28)
Using the homomorphism resultant algorithm [25] eliminates ω from h and h , and the cross points of (26) are calculated with the real solutions of the resulting equation for T . The real solutions for T and ω are T = 0.206109661 T = 0.987347284
ω = 0.957875475 ω = 0.063412910.
(29) (30)
724
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: EXPRESS BRIEFS, VOL. 61, NO. 9, SEPTEMBER 2014
Fig. 4.
Unity feedback control system.
Assume that a fractional phase-lead controller is used and that its transfer function is p1 (sp4 + p2 ) . (34) K(s) = (sp4 + p3 ) A schematic of the closed-loop system with a controller is given in Fig. 4. The following design solution is found by the modified backpropagation algorithm in [26]: p1 = 9.24
p2 = 7.513
p3 = 15.204
p4 = 1.1.
The closed-loop transfer function is given by √ 2p1 (sp4 +p2 )e−τ s √ √ . (35) H(s) = √ p s(s 4 +p3 )(1−e−2τ s )+2p1(sp4 +p2 )e−τ s The corresponding quasi-polynomial characteristic equation is obtained as
Fig. 2. Root-loci for CE(s, τ ) in (26) from τ = 0.1 to τ = 10.
√
−τ s CE(s, τ ) = (s1.6 +15.204s0.5 )+(18.48s1.1 +138.8402)e √ 1.6 0.5 −2τ s − (s +15.204s )e = 0. (36)
Furthermore, in view of relation (14), we get
Fig. 3. Step response of G(s) in Example 1.
The corresponding time delays of the cross points in (29) and (30) are obtained with regard to relation (11) as τ1 = 0.587420879 τ2 = 8.750064088.
(31) (32)
By applying the criterion expressed in the previous section, it is easy to find out that a stabilizing crossing of roots (RT = −1) has occurred at τ = 0.587420879 for s = ±j0.957875475 and that a destabilizing crossing (RT = +1) has taken place at τ = 8.750064088 for s = ±j0.06341291. Since the system is unstable for τ = 0, the only stability window for this system is 0.587420879 < τ < 8.750064088. In Fig. 2, the root-loci curve of this system for the changes in τ from τ = 0.1 (dark blue) to τ = 10 (dark red) has been presented for a better perception of the system. In Fig. 3, the step responses of G(s) has been plotted for τ = 0.5, τ = 3, τ = 5, τ = 8.75, and τ = 10. Example 2: Consider the plant whose transfer function G(s) is given by 1 √ . G(s) = √ s sinh(τ s)
(33)
The nonrational transfer function in (33) occurs when the plant is governed by a heat conduction (or diffusion) equation [26].
h(v, T ) = (v 16 + 15.204v 5 )(1 + jT )2 + (18.48v 11 + 138.8402)r(1 + T 2 ) − (v 16 + 15.204v 5 )r2 (1 − jT )2 = 0. (37) √ By inserting expression v = s1/10 = 10 ω ejπ/20 = w ejπ/20 into (37) and equating the real and imaginary parts of the obtained relation to zero, we get √ 1 + 5 16 15.204 5 h = − w + √ w (1 − T 2 )(1 − r2 ) 4 2 √ √ 5 5 − 5 16 √ w + 15.204 2w (1 + r2 )T − 2 √ √ √ 5− 5 2( 5 + 1) − w11 (1 + T 2 )r + 18.48 4 8 (38) + + T 2 )r = 0 138.8402(1 √ √ 1 + 5 16 h = − w + 15.204 2w5 (1 + r2 )T 2 √ 5 − 5 16 15.204 5 √ w + √ w (1 − T 2 )(1 − r2 ) + 2 2 2 √ √ √ 5− 5 2( 5+1) + w11 (1+T 2 )r = 0 + 18.48 4 8 (39) √ where w = 10 ω and r = exp(−2 tan−1 T ). Through the “resultant” command in the Maple software package, the resultant algorithm has been presented in [25] or in (19), w can be eliminated from (38) and (39), and the real and positive values of T and then ω = w10 can be calculated as ω = 0.67120191800
T = 2.60230883324.
(40)
The corresponding time delays of the cross points in (36) are obtained with regard to relation (11) as τ = 4.1563867684.
(41)
PAKZAD AND NEKOUI: ANALYSIS OF LINEAR TIME-INVARIANT FRACTIONAL EXPONENTIAL DELAY SYSTEMS
725
R EFERENCES
Fig. 5.
Root-loci curve for CE(s, τ ) in (36) from τ = 0.1 to τ = 5.
Fig. 6.
Step responses of H(s).
By applying the previously described method, it is realized that the crossing through the imaginary axis happens at τ = 4.1563867684 for ω = 0.671201918, i.e., destabilizing crosses (RT = +1). As the original system was delay-free stable, then it means that the only stability window for this system is 0 ≤ τ < 4.1563867683. To get a better understanding of the properties of this system, the root-loci curve for CE(s, τ ) in (36) has been plotted as a function of delay in Fig. 5. The color spectrum indicating the changes in τ from τ = 0.1 (dark blue) to τ = 5 (dark red) has been illustrated in the “color bar.” In addition, Fig. 6 shows the step responses of H(s) for τ = 1.15, τ = 2.15, τ = 3.15, τ = 4.15, and τ = 5. IV. C ONCLUSION A new method for the stability analysis of a large class of linear fractional exponential systems for both single- and commensurate-delay cases is proposed. In this method, a characteristic equation was converted into an auxiliary polynomial. 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. Finally, two illustrative examples are presented to highlight the proposed approach.
[1] M. A. Pakzad, S. Pakzad, and M. A. Nekoui, “Stability analysis of multiple time delayed fractional order systems,” in Proc. Amer. Control Conf., Washington, DC, USA, 2013, pp. 170–175. [2] N. Ozturk and A. Uraz, “An analysis stability test for a certain class of distributed parameter systems with delays,” IEEE Trans. Circuits Syst., vol. 34, no. 4, pp. 393–396, Apr. 1985. [3] M. A. Pakzad and M. A. Nekoui, “Direct method for stability analysis of fractional delay systems,” Int. J. Comput., Commun. Control, vol. 8, no. 6, pp. 863–868, Dec. 2013. [4] N. Ozturk, “Stability independent of distributed lag for a special class of distributed parameter systems,” in Proc. 34th IEEE Conf. Decis. Control, 1995, pp. 3245–3246. [5] S. Pakzad, M. A. Pakzad, and M. A. Nekoui, “Stability map of fractional delay systems in the parametric space of delays and coefficient,” in Proc. Amer. Control Conf., Washington, DC, USA, 2013, pp. 176–181. [6] C. Bonnet and J. R. Partington, “Analysis of fractional delay systems of retarded and neutral type,” Automatica, vol. 38, no. 7, pp. 1133–1138, Jul. 2002. [7] S. Pakzad and M. A. Pakzad, “Stability condition for discrete systems with multiple state delays,” WSEAS Trans. Syst. Control, vol. 6, no. 11, pp. 417–426, Nov. 2011. [8] C.-C. Hua, D. Liu, and X.-P. Guan, “Necessary and sufficient stability criteria for a class of fractional-order delayed systems,” IEEE Trans. Circuits Syst. II, Exp. Briefs, vol. 61, no. 1, pp. 59–63, Jan. 2014. [9] M. A. Pakzad and B. Moaveni, “Delay-Dependent state estimation for time delay systems,” WSEAS Trans. Syst. Control, vol. 8, no. 1, pp. 1–10, Jan. 2013. [10] M. A. Pakzad, “Kalman filter design for time delay systems,” WSEAS Trans. Syst., vol. 11, no. 10, pp. 551–560, Oct. 2012. [11] M. A. Pakzad and M. A. Nekoui, “Stability map of multiple time delayed fractional order systems,” Int. J. Control Autom. Syst., vol. 12, no. 1, pp. 37–43, Feb. 2014. [12] C. Hwang and Y. C. Cheng, “A numerical algorithm for stability testing of fractional delay systems,” Automatica, vol. 42, no. 5, pp. 825–831, May 2006. [13] M. A. Pakzad and S. Pakzad, “Stability map of fractional order time-delay systems,” WSEAS Trans. Syst., vol. 10, no. 11, pp. 541–550, Oct. 2012. [14] K. E. Walton and J. E. Marshal, “Direct method for TDS stability analysis,” Proc. Inst. Elect. Eng.-D, vol. 134, no. 2, pp. 101–107, Mar. 1987. [15] M. Shi and Z. H. Wang, “An effective analytical criterion for stability testing of fractional-delay systems,” Automatica, vol. 47, no. 9, pp. 2001– 2005, Sep. 2011. [16] M. A. Pakzad, M. A. Nekoui, and S. Pakzad, “Stability analysis of timedelayed linear fractional-order systems,” Int. J. Control Autom. Syst., vol. 11, no. 3, pp. 519–525, Jun. 2013. [17] M. A. Pakzad, S. Pakzad, and M. A. Nekoui, “On stability of linear timeinvariant fractional order systems with multiple delays,” in Proc. Amer. Control Conf., Portland, OR, USA, 2014, pp. 4883–4888. [18] M. A. Pakzad, S. Pakzad, and M. A. Nekoui, “On the stability of fractional order systems of neutral Type,” J. Comput. Nonlinear Dynam., 2014. doi:10.1115/1.4027593, to be published. [19] E. Weber, Linear Transient Analysis, vol. 2. New York, NY, USA: Wiley, 1956. [20] R. F. Curtain and H. J. Zwart, “An Introduction to Infinite-Dimensional Linear Systems Theory,” in Text in Applied Mathematics, vol. 21. Berlin, Germany: Springer-Verlag, 1995. [21] J.-J. Loiseau and H. Mounier, “Stabilisation de lequation de la chaleur commandee en flux,” in Proc. ESAIM, vol. 5, Systemes Differentiels Fractionnaires, Modeles, Methodes et Applications, 1998, pp. 131–144. [22] C. Bonnet and J. R. Partington, “Stabilization of fractional exponential systems including delays,” Kybernetika, vol. 37, no. 3, pp. 345–353, 2001. [23] R. Curtain and K. Morris, “Transfer functions of distributed parameter systems: A tutorial,” Automatica, vol. 45, no. 5, pp. 1101–1116, May 2009. [24] F. Merrikh-Bayat, “General formula for stability testing of linear systems with fractional-delay characteristic equation,” Central Eur. J. Phys., vol. 11, no. 6, Jun. 2013. [25] G. E. Collins, “The calculation of multivariate polynomial resultants,” J. Assoc. Comput. Mach., vol. 18, no. 4, pp. 515–532, 1971. [26] S. Arunsawatwong and V. Q. Nguyen, “Design of retarded fractional delay differential systems using the method of inequalities,” Int. J. Autom. Comput., vol. 6, no. 1, pp. 22–28, Feb. 2009.
721
Stability Analysis of Linear Time-Invariant Fractional Exponential Delay Systems Mohammad Ali Pakzad and Mohammad Ali Nekoui
Abstract—This brief presents a new approach for the stability analysis of linear fractional exponential delay systems with commensurate orders and multiple commensurate delays that enable us to decide on some cases that were previously open problems. In the proposed approach, first an auxiliary polynomial is generated by mapping the principal sheet of the Riemann surface and a pseudodelay transformation. Next, this auxiliary polynomial is employed in determining all possible purely imaginary characteristic roots for any positive time delay. Then, the concept of stability is expressed as a function of delay. Two illustrative examples are provided to demonstrate the effectiveness of the proposed method and to gain a better understanding of the problem, and the root-locus curve of these systems has been plotted as a function of time delay. Index Terms—Fractional exponential systems, fractional-order systems, root locus, stability analysis, time-delay systems.
I. I NTRODUCTION
T
IME delay is an inherent part of many dynamic and physical systems. It is known that time delay could lead to inadequate performance and even instability of systems. Therefore, time-delay systems play significant roles in theoretical [1]–[7] and practical fields [8]–[10]. There are many linear infinite-dimensional systems with transfer functions that can be represented as fractional differential systems. When time delays and a fractional-order derivative appear in dynamic systems, we have fractional-delay systems. The characteristic equation of a fractional-delay system involves exponential-type transcen√ dental terms, i.e., e−τ s or e−τ s ; thus, a fractional-delay system has, in general, an infinite number of characteristic roots. In this brief, one of the most important and unresolved problems of fractional-delay systems is studied, i.e., the asymptotic stability of a general class of linear fractional exponential delay systems. Ozturk and Uraz [2], [4] are probably the pioneers in the stability investigation of fractional-order time-delay systems. They have developed the Ruth–Hurwitz criterion to evaluate the stability of a √ class of distributed systems whose transfer √ functions involve s and/or e−τ s . In [6], the necessary and sufficient conditions for the bounded-input–bounded-output Manuscript received January 17, 2014; revised April 30, 2014; accepted June 27, 2014. Date of publication July 8, 2014; date of current version September 1, 2014. This brief was recommended by Associate Editor J. Lu. M. A. Pakzad is with the Department of Electrical Engineering, Science and Research Branch, Islamic Azad University, Tehran 14778 93855, Iran (e-mail: [email protected]). M. A. Nekoui is with the Faculty of Electrical and Computer Engineering, Khajeh Nasir Toosi University of Technology, Tehran 19697 64499, Iran (e-mail: [email protected]). Color versions of one or more of the figures in this brief are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TCSII.2014.2335427
(BIBO) stability of fractional-order delay systems have been introduced. From the numerical analysis point of view, the effective numerical algorithms have been discussed in [11] for the evaluation of the stability of fractional-order delay systems. In [12], a heavy computation scheme based on Cauchy’s integral has been proposed to test the stability of such systems. Pakzad and Pakzad [13] have successfully extended the direct analytical method (presented by Walton and Marshal [14] for the stability analysis of integer-order time-delay systems) to fractionalorder delay systems. The stability of these systems can be also analytically analyzed by employing Hassard’s theorem proposed by Shi and Wang [15]. Based on the method they presented, the number of unstable poles of the characteristic equation for a given delay is determined. Recently, in [16], an analytical method for finding the stability regions of fractionaldelay systems with commensurate orders has been presented, which uses the bilinear Rekasius transformation to eliminate the exponential-type transcendental term in the characteristic equations of these systems. To gain more information about some of the recently developed methods and tools in this area, see [17] and [18]. A large number of circuits and systems have distributed parameters and/or delay elements. Thermal processes, the hole diffusion of transistors, electromagnetic devices, and transmission lines are typical examples. In some cases, the mathematical models in the Laplace-transform domain √ for the transfer function of a system may even consist of s in combination √ with fractional-delay term e−τ s . Many examples of fractional differential systems can √be found in literature. Simple examples, such as G(s) = e−τ s /s, with τ > 0 arising in the theory of transmission √ lines, √ are given in [19]. As another example, G(s) = tanh( s)/ s appears in a boundary-controlled and observed diffusion process in a bounded√domain In addi√ [20]. √ tion, transfer functions G(s) = cosh(x s)/ s sinh( s) and 0 √ √ G(s) = sinh(x0 s)/ sinh( s) appear in the boundary control of a 1-D heat equation with Neumann–Dirichlet √ √ boundary con√ ditions [20], or G(s) = 2e−τ s /b s(1−e−2τ s ) can be found in [21]. Other examples of this type can be found in [20]–[23]. The branch point of the transfer function of a system may be at a point rather than at the origin. As an example of this type, consider the asymmetric transmission line shown in Fig. 1. In this figure, R, C, L, and G show the resistance, capacitance, inductance, and conductance, respectively, in the unit of length, and the length of the line is assumed equal to L. According to Krishoff’s law, the equations governing this system are v(x, t) di(x, t) = Gv(x, t) + C (1) − dt dt i(x, t) dv(x, t) = Ri(x, t) + L (2) − dt dt
1549-7747 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
722
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: EXPRESS BRIEFS, VOL. 61, NO. 9, SEPTEMBER 2014
II. M ETHODOLOGY
Fig. 1. Model of an asymmetric transmission line [24]. Transfer function Vo (s)/Vs (s) has branch points at R/L and −G/C.
subject to the following boundary conditions: v(0, t) = vs (t)
i(L, t) = 0.
(3)
Taking the Laplace transform from both sides of (1)–(3) (assuming that the initial condition of the system is equal to zero) and after algebraic calculations, the transfer function of the system is obtained as follows [24]: G(s) =
Vo (s) 1 . = 1 1 Vs (s) cosh L(R + Ls) 2 (G + Cs) 2
(4)
In literature, the fractional-delay systems are those described by the following transfer function: 1 r q (s)e−βk s q0 (s) + nk=1 n2 k (5) G(s) = p0 (s) + k=1 pk (s)e−γk sr where r is a real number such that 0 < r < 1, pk are of the form lk mk αi δi i=1 qk ai s , and qk are of the form i=1 qk bi s . The fractional-delay characteristic equation is obtained by equating the denominator of a fractional-delay transfer function to zero. More formally, the fractional-delay characteristic equation under consideration in this brief can be expressed in the general form of CE(sα , τ ) = p0 (sα ) +
n
P (sα )e−τ
The main objective of this section is to present a new method for the evaluation of stability and the determination of the unstable roots of a fractional-order delay system. In this approach, we are trying to calculate the values of τ where there exists a crossing of poles through the imaginary axis. With these values in hand, we will be able to calculate the direction of crossing from the left half-plane to the right half-plane, which we will denote as a destabilizing crossing, or from the right to the left, which is a stabilizing crossing.
√
s
(6)
=1
where 0 < α < 1. We find out from the work in [22] that the transfer function of a system with a characteristic equation in the form of (6) will be BIBO stable if and only if it does not have any pole at (s) ≥ 0 (in particular, no poles of fractional order at s = 0). For fractional-order systems, if an auxiliary variable of v = sα is used, a practical test for the evaluation of stability can be obtained. In this new variable, the instability region of the original system is not given by the right half-plane but in fact by the region described as απ |∠v| ≤ (7) 2
A. Imaginary Axis Crossings In order to find the points on the imaginary axis where the crossing takes place, the auxiliary variable described in the previous section and a new substitution are used. As a result, the characteristic equation is converted into a solvable equation. To simplify this characteristic equation, we propose a new transformation that is defined as √ 1 − jT e−τ s , T ∈R (10) =r s=jω 1 + jT where r = exp(−2 tan−1 T ). It is important to note that this √ substitution is an exact expression of e−τ s for purely imaginary roots s = ±jω. Note that, for s = ±jω, the magnitudes of both sides in (10) always agree with each other. Furthermore, for the transformation to exactly hold, it is required that the phase condition developed from (10) holds. This condition can be found as √ 2 2 (11) τ = √ tan−1 (T ). ω which is the inverse mapping from the (T, ω) domain to τ . Note that there is a one-to-one mapping between T , ω, and τ . To summarize, the transformation in (10) is indeed exact at s = jω as long as (11) holds. By inserting (10) into (6), we have n 1 − jT CE(sα , T ) = p0 (sα ) + P (sα )r . (12) 1 + jT =1
By applying this auxiliary variable v = sα in (12), the following relation is obtained: n 1 − jT P (v)r . (13) CEv (v, T ) = p0 (v) + 1 + jT =1
with v ∈ C. Note that, under this transformation, the imaginary axis in the s-domain is mapped into απ (8) ∠v = ± . 2
This will transform the domain of the system from a multisheeted Riemann surface into a complex plane, in which the poles can be easily calculated. By multiplying (13) by (1 + jT )n , the new form of the characteristic equation is reached as
Let us assume that s = ±jω or, in other words, s = ωe±jπ/2 are the roots of characteristic (6) for τ ∈ R+ . Then, for the auxiliary variable, the roots are defined as απ απ ± j sin . (9) v = sα = ω α e±jαπ/2 = ω α cos 2 2
h(v, T ) = (1 + jT )n Cv (v, T ) = (1 + jT )n p0 (v) n P (v)r (1 − jT ) (1 + jT )(n−) . (14) +
Therefore, with auxiliary variable v = sα , there is a direct relation between the roots on the imaginary axis for the s-domain with the roots having argument ±απ/2 in the v-domain.
This expression is a polynomial in v of which the coefficients are parametric functions of T and r. When s = jω are the roots of (6) for some ω and τ ∈ R, then v = ω α ejαπ/2 will be also
=0
PAKZAD AND NEKOUI: ANALYSIS OF LINEAR TIME-INVARIANT FRACTIONAL EXPONENTIAL DELAY SYSTEMS
the roots of (14) for some ω and T ∈ R. To find the crossing frequencies from the imaginary axis in (6), v = ω α ejαπ/2 should be inserted into relation (14) and then the real and imaginary parts of the resulting equation should be separated as h(v, T ))|
jαπ v=ω α e 2
= h (ω α , T ) + jh (ω α , T ) = 0. (15)
In the given relations, h = (h) and h = (h). By replacing ω α with w in the given equation, we have h(w, T ) = h (w, T ) + jh (w, T ) = 0.
(16)
Making the real and imaginary parts of (16) equal to zero results in two algebraic equations; let us investigate those (ω, T ) solutions from h = h =
n i=0 n
ai (T )wi = 0
(17)
bi (T )wi = 0.
(18)
i=0
We utilize the resultant theory (see [25]) to eliminate w from the two multivariate polynomials h and h . Definition 1: Consider the two multivariate polynomials (17) and (18) in terms of w and T with real coefficients, where h and h have positive degrees in terms of w, and n > 0. The result of h and h with respect to w is defined by an (T ) an−1 (T ) . . . . . . . . . 0 .. .. .. .. .. .. . . . . . . .. ... . . . . a1 (T ) a0 (T ) 0 Rw (h , h ) = (19) 0 bn (T ) bn−1 (T ) . . . . . . . . . . .. .. .. .. .. . . . . . . . .. 0 ... . . . . b1 (T ) b0 (T ) which is the determinant of the well-known Sylvester matrix [25]. A 2n-order Sylvester matrix is constructed via (19) and its determinant Rw (h , h ) is a function of T as D(T ) = Rw (h , h ) = 0.
(20)
Hence, the set of all w and T that makes both equations zero can be obtained. Then, for every ω = w1/α and T , its corresponding τ can be determined using (11). The whole ω value, for which s = jω is a root of (6) for some nonnegative delays, is defined as the following crossing frequency set: 1
Ω = ω ∈ R+ C s α , τ = 0, for some τ ∈ R+ . (21) Corollary 1: If the system given as (12) is stable for τ = 0 (i.e., a system without delay) and Ω = φ, then the system will be stable for all positive values of τ . Proof 1: From the fact that there are no roots crossing the imaginary axis. This class of systems only exhibits a finite number of possible imaginary characteristic roots for any τ ∈ R+ at the given frequencies. Moreover, this method detected all of them. Let us call this set {ωc } = {ωc1 , ωc2 , ωc3 , . . . , ωcn }
(22)
723
where subscript c refers to “crossing” the imaginary axis. Furthermore, to each ωcm , m = 1, 2, . . . , n finitely correspond to the following τ values: {τc } = {τc1 , τc2 , τc3 , . . . , τcn }.
(23)
B. Direction of Crossing After the crossing points of characteristic (6) from the imaginary axis are obtained, the goal now is to determine whether each of these root crossings from the imaginary axis is a stabilizing cross or a destabilizing cross. Assume that (s, τ ) is a simple root of CE(sα , τ ) = 0. The root sensitivities associated with each purely imaginary characteristic root jω with respect to τ is defined as ∂CE ds s ∂τ Sτ |s=jωc = = − ∂CE . (24) dτ s=jωc ∂s s=jωc
The root tendency for each ωcm and τmk is defined as Root Tendency = RT|s=jωc = sign Sτs | s=jωcm . (25) τ =τmk
If it is positive, then it is a destabilizing crossing, whereas if it is negative, it is a stabilizing crossing. III. I LLUSTRATIVE E XAMPLES We present two cases, which display all the features discussed in this brief. Example 1: Consider the following linear time-invariant fractional-delay system: √
0.3e−τ s √ . G(s) = √ 6 √ 4 √ ( s) + 0.2( s) + ( s)2 + 0.3e−τ s The corresponding characteristic equation of G(s) is √ √ √ √ CE(s, τ ) = ( s)6 + 0.2( s)4 + ( s)2 + 0.3e−τ s . (26) This example has been taken from the work in [4]. The characteristic in (26) has poles in s = 0.0461 ± j1.0123 and s = −0.2921 for τ = 0. Therefore, this system is unstable for τ = 0. Our objective in this example is to find stability windows based on the method described in this brief for this system. By applying the criterion expressed in the previous section and by using auxiliary variable v = s0.5 , we can convert (26), according to (14), as h(v, τ ) = (v 6 + v 4 + v 2 )(1 − jT ) + 0.3r(1 − jT ). (27) √ √ √ By√ inserting expression v = ω ej(π/4) = ω(( 2/2) + j( 2/2)) into (27) and equating the real and imaginary parts of the obtained relation to zero, we get h = T ω 3 − 0.2ω 2 − T ω + 0.3r = 0 h = −ω 3 − 0.2T ω 2 + ω − 0.3T r = 0.
(28)
Using the homomorphism resultant algorithm [25] eliminates ω from h and h , and the cross points of (26) are calculated with the real solutions of the resulting equation for T . The real solutions for T and ω are T = 0.206109661 T = 0.987347284
ω = 0.957875475 ω = 0.063412910.
(29) (30)
724
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: EXPRESS BRIEFS, VOL. 61, NO. 9, SEPTEMBER 2014
Fig. 4.
Unity feedback control system.
Assume that a fractional phase-lead controller is used and that its transfer function is p1 (sp4 + p2 ) . (34) K(s) = (sp4 + p3 ) A schematic of the closed-loop system with a controller is given in Fig. 4. The following design solution is found by the modified backpropagation algorithm in [26]: p1 = 9.24
p2 = 7.513
p3 = 15.204
p4 = 1.1.
The closed-loop transfer function is given by √ 2p1 (sp4 +p2 )e−τ s √ √ . (35) H(s) = √ p s(s 4 +p3 )(1−e−2τ s )+2p1(sp4 +p2 )e−τ s The corresponding quasi-polynomial characteristic equation is obtained as
Fig. 2. Root-loci for CE(s, τ ) in (26) from τ = 0.1 to τ = 10.
√
−τ s CE(s, τ ) = (s1.6 +15.204s0.5 )+(18.48s1.1 +138.8402)e √ 1.6 0.5 −2τ s − (s +15.204s )e = 0. (36)
Furthermore, in view of relation (14), we get
Fig. 3. Step response of G(s) in Example 1.
The corresponding time delays of the cross points in (29) and (30) are obtained with regard to relation (11) as τ1 = 0.587420879 τ2 = 8.750064088.
(31) (32)
By applying the criterion expressed in the previous section, it is easy to find out that a stabilizing crossing of roots (RT = −1) has occurred at τ = 0.587420879 for s = ±j0.957875475 and that a destabilizing crossing (RT = +1) has taken place at τ = 8.750064088 for s = ±j0.06341291. Since the system is unstable for τ = 0, the only stability window for this system is 0.587420879 < τ < 8.750064088. In Fig. 2, the root-loci curve of this system for the changes in τ from τ = 0.1 (dark blue) to τ = 10 (dark red) has been presented for a better perception of the system. In Fig. 3, the step responses of G(s) has been plotted for τ = 0.5, τ = 3, τ = 5, τ = 8.75, and τ = 10. Example 2: Consider the plant whose transfer function G(s) is given by 1 √ . G(s) = √ s sinh(τ s)
(33)
The nonrational transfer function in (33) occurs when the plant is governed by a heat conduction (or diffusion) equation [26].
h(v, T ) = (v 16 + 15.204v 5 )(1 + jT )2 + (18.48v 11 + 138.8402)r(1 + T 2 ) − (v 16 + 15.204v 5 )r2 (1 − jT )2 = 0. (37) √ By inserting expression v = s1/10 = 10 ω ejπ/20 = w ejπ/20 into (37) and equating the real and imaginary parts of the obtained relation to zero, we get √ 1 + 5 16 15.204 5 h = − w + √ w (1 − T 2 )(1 − r2 ) 4 2 √ √ 5 5 − 5 16 √ w + 15.204 2w (1 + r2 )T − 2 √ √ √ 5− 5 2( 5 + 1) − w11 (1 + T 2 )r + 18.48 4 8 (38) + + T 2 )r = 0 138.8402(1 √ √ 1 + 5 16 h = − w + 15.204 2w5 (1 + r2 )T 2 √ 5 − 5 16 15.204 5 √ w + √ w (1 − T 2 )(1 − r2 ) + 2 2 2 √ √ √ 5− 5 2( 5+1) + w11 (1+T 2 )r = 0 + 18.48 4 8 (39) √ where w = 10 ω and r = exp(−2 tan−1 T ). Through the “resultant” command in the Maple software package, the resultant algorithm has been presented in [25] or in (19), w can be eliminated from (38) and (39), and the real and positive values of T and then ω = w10 can be calculated as ω = 0.67120191800
T = 2.60230883324.
(40)
The corresponding time delays of the cross points in (36) are obtained with regard to relation (11) as τ = 4.1563867684.
(41)
PAKZAD AND NEKOUI: ANALYSIS OF LINEAR TIME-INVARIANT FRACTIONAL EXPONENTIAL DELAY SYSTEMS
725
R EFERENCES
Fig. 5.
Root-loci curve for CE(s, τ ) in (36) from τ = 0.1 to τ = 5.
Fig. 6.
Step responses of H(s).
By applying the previously described method, it is realized that the crossing through the imaginary axis happens at τ = 4.1563867684 for ω = 0.671201918, i.e., destabilizing crosses (RT = +1). As the original system was delay-free stable, then it means that the only stability window for this system is 0 ≤ τ < 4.1563867683. To get a better understanding of the properties of this system, the root-loci curve for CE(s, τ ) in (36) has been plotted as a function of delay in Fig. 5. The color spectrum indicating the changes in τ from τ = 0.1 (dark blue) to τ = 5 (dark red) has been illustrated in the “color bar.” In addition, Fig. 6 shows the step responses of H(s) for τ = 1.15, τ = 2.15, τ = 3.15, τ = 4.15, and τ = 5. IV. C ONCLUSION A new method for the stability analysis of a large class of linear fractional exponential systems for both single- and commensurate-delay cases is proposed. In this method, a characteristic equation was converted into an auxiliary polynomial. 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. Finally, two illustrative examples are presented to highlight the proposed approach.
[1] M. A. Pakzad, S. Pakzad, and M. A. Nekoui, “Stability analysis of multiple time delayed fractional order systems,” in Proc. Amer. Control Conf., Washington, DC, USA, 2013, pp. 170–175. [2] N. Ozturk and A. Uraz, “An analysis stability test for a certain class of distributed parameter systems with delays,” IEEE Trans. Circuits Syst., vol. 34, no. 4, pp. 393–396, Apr. 1985. [3] M. A. Pakzad and M. A. Nekoui, “Direct method for stability analysis of fractional delay systems,” Int. J. Comput., Commun. Control, vol. 8, no. 6, pp. 863–868, Dec. 2013. [4] N. Ozturk, “Stability independent of distributed lag for a special class of distributed parameter systems,” in Proc. 34th IEEE Conf. Decis. Control, 1995, pp. 3245–3246. [5] S. Pakzad, M. A. Pakzad, and M. A. Nekoui, “Stability map of fractional delay systems in the parametric space of delays and coefficient,” in Proc. Amer. Control Conf., Washington, DC, USA, 2013, pp. 176–181. [6] C. Bonnet and J. R. Partington, “Analysis of fractional delay systems of retarded and neutral type,” Automatica, vol. 38, no. 7, pp. 1133–1138, Jul. 2002. [7] S. Pakzad and M. A. Pakzad, “Stability condition for discrete systems with multiple state delays,” WSEAS Trans. Syst. Control, vol. 6, no. 11, pp. 417–426, Nov. 2011. [8] C.-C. Hua, D. Liu, and X.-P. Guan, “Necessary and sufficient stability criteria for a class of fractional-order delayed systems,” IEEE Trans. Circuits Syst. II, Exp. Briefs, vol. 61, no. 1, pp. 59–63, Jan. 2014. [9] M. A. Pakzad and B. Moaveni, “Delay-Dependent state estimation for time delay systems,” WSEAS Trans. Syst. Control, vol. 8, no. 1, pp. 1–10, Jan. 2013. [10] M. A. Pakzad, “Kalman filter design for time delay systems,” WSEAS Trans. Syst., vol. 11, no. 10, pp. 551–560, Oct. 2012. [11] M. A. Pakzad and M. A. Nekoui, “Stability map of multiple time delayed fractional order systems,” Int. J. Control Autom. Syst., vol. 12, no. 1, pp. 37–43, Feb. 2014. [12] C. Hwang and Y. C. Cheng, “A numerical algorithm for stability testing of fractional delay systems,” Automatica, vol. 42, no. 5, pp. 825–831, May 2006. [13] M. A. Pakzad and S. Pakzad, “Stability map of fractional order time-delay systems,” WSEAS Trans. Syst., vol. 10, no. 11, pp. 541–550, Oct. 2012. [14] K. E. Walton and J. E. Marshal, “Direct method for TDS stability analysis,” Proc. Inst. Elect. Eng.-D, vol. 134, no. 2, pp. 101–107, Mar. 1987. [15] M. Shi and Z. H. Wang, “An effective analytical criterion for stability testing of fractional-delay systems,” Automatica, vol. 47, no. 9, pp. 2001– 2005, Sep. 2011. [16] M. A. Pakzad, M. A. Nekoui, and S. Pakzad, “Stability analysis of timedelayed linear fractional-order systems,” Int. J. Control Autom. Syst., vol. 11, no. 3, pp. 519–525, Jun. 2013. [17] M. A. Pakzad, S. Pakzad, and M. A. Nekoui, “On stability of linear timeinvariant fractional order systems with multiple delays,” in Proc. Amer. Control Conf., Portland, OR, USA, 2014, pp. 4883–4888. [18] M. A. Pakzad, S. Pakzad, and M. A. Nekoui, “On the stability of fractional order systems of neutral Type,” J. Comput. Nonlinear Dynam., 2014. doi:10.1115/1.4027593, to be published. [19] E. Weber, Linear Transient Analysis, vol. 2. New York, NY, USA: Wiley, 1956. [20] R. F. Curtain and H. J. Zwart, “An Introduction to Infinite-Dimensional Linear Systems Theory,” in Text in Applied Mathematics, vol. 21. Berlin, Germany: Springer-Verlag, 1995. [21] J.-J. Loiseau and H. Mounier, “Stabilisation de lequation de la chaleur commandee en flux,” in Proc. ESAIM, vol. 5, Systemes Differentiels Fractionnaires, Modeles, Methodes et Applications, 1998, pp. 131–144. [22] C. Bonnet and J. R. Partington, “Stabilization of fractional exponential systems including delays,” Kybernetika, vol. 37, no. 3, pp. 345–353, 2001. [23] R. Curtain and K. Morris, “Transfer functions of distributed parameter systems: A tutorial,” Automatica, vol. 45, no. 5, pp. 1101–1116, May 2009. [24] F. Merrikh-Bayat, “General formula for stability testing of linear systems with fractional-delay characteristic equation,” Central Eur. J. Phys., vol. 11, no. 6, Jun. 2013. [25] G. E. Collins, “The calculation of multivariate polynomial resultants,” J. Assoc. Comput. Mach., vol. 18, no. 4, pp. 515–532, 1971. [26] S. Arunsawatwong and V. Q. Nguyen, “Design of retarded fractional delay differential systems using the method of inequalities,” Int. J. Autom. Comput., vol. 6, no. 1, pp. 22–28, Feb. 2009.
