Stability Analysis of Multiple Time Delayed Fractional Order Systems
2013 American Control Conference (ACC) Washington, DC, USA, June 17-19, 2013
Stability Analysis of Multiple Time Delayed Fractional Order Systems Mohammad Ali Pakzad1 , Sara Pakzad2 , and Mohammad Ali Nekoui3
Obviously, the accuracy and completeness of the analysis strongly depends on finding all the existing PSSC without any approximations [1]. There has been a large effort to deal with this problem, for the standard case (integer order systems); see [1]-[6], and others. The researchers of [7] and [8] are probably the pioneer to consider stability analysis of the fractional delay system with single-delay. They have extended the Ruth-Hurwitz criteria for analyzing the stability of some √ special delay systems to those involve fractional power s. For single delay case; In [9], necessary and sufficient conditions for BIBO stability of the retarded fractional-order delay systems and sufficient conditions for some neutral types have been introduced; and also from the analitical analysis point of view, the effective analitical approaches to check BIBO stability of fractionalorder delay systems have been discussed in [17] and [18]. The original idea in this strategy is derived based on the method (advanced clustering with frequency sweeping) reported in Reference [1] to achieve the stability map of integer order systems with multiple time delays. In this paper, we extend the approach of [1] to fractionalorder systems with multiple delays. Necessary and sufficient conditions which yield the exact lower and upper bounds of the crossing frequency set (CFS) can be computed via an automated sequential formula. These bounds are crucial as they determine the sweeping range of the only parameter, the frequency that ACFS sweeps.
Abstract— A new methodology, based on advanced clustering with frequency sweeping (ACFS), is presented for the stability analysis of fractional-order systems with multiple time delays against delay uncertainties. The problem is known to be notoriously complex, primarily because the systems are infinite dimensional due to delays. Multiplicity of the delays in this study complicates the analysis even further. And “fractionalorder” feature of the systems makes the problem much more challenging compared to integer order systems. ACFS does not impose any restrictions in the number of delays, and it can directly extract the 2-D cross sections of the stability views in any two delay domain. We show that this procedure analytically reveals all possible stability regions exclusively in the space of the delays. As an added strength, it does not require the delay-free system under consideration to be stable. The main contribution of this document is that we demonstrate for the first time that the stability maps of a fractional-order system with multiple time delays can be obtained efficiently. Two illustrative examples are presented to confirm the proposed method results.
I. INTRODUCTION In a generic sense, the fractional-order systems are identified by non-integer powers of the Laplace variable s. When time delays and fractional-order derivative are involved in dynamical systems, we have fractional-delay systems. The characteristic equation of a fractional-delay system involves exponential type transcendental terms, so a fractional delay system has infinitely many roots. This makes the stability analysis of fractional-delay systems a challenging task. In this paper, one of the most important and unresolved problems of fractional delay systems is studied: the asymptotic stability of a general class of fractional-order systems with multiple time delays, against delay uncertainties. We make two qualitative comments first. (a) There is no methodology presently to the knowledge of the authors, which can reveal the aimed “stability map” of this system. (b) As a natural consequence of (a) there is no methodology to present this stability map exhaustively. Exhaustiveness means including all the stable regions in a given segment of interest in (τ1 , τ2 ) space attached to or detached from the origin (0, 0). The main objective in 2-D stability analysis is to construct all the potential stability switching curves (PSSC) which partition the delay space into stable and unstable regions.
II. PRELIMINARIES AND DEFINITIONS Standard notation has been used throughout the article: N , Z : sets of natural and integer numbers. R(R+ , R− ) : set of real (positive real, negative real) numbers. C:√ set of complex numbers. j = −1 : the imaginary unit. z : complex conjugate of a complex number z ℜ(z) : real part of a complex number z ℑ(z) : imaginary part of a complex number z |z|, ∠z:magnitude and argument of a complex number z RTℓ (p1 , p2 ): denotes the resultant of bivariate polynomials p1 (T1 , T2 ) and p2 (T1 , T2 ) with eliminating Tℓ , where ℓ = 1, 2 Consider a fractional-order system with the following characteristic equation:
1 M.A. Pakzad is with Department of Electrical Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
n
C(s, τ1 , τ2 ) =
[email protected] 2 S.
South
ℓ=0
Pakzad is with the Departments of Electrical Engineering, Tehran Branch, Islamic Azad University, Tehran, Iran
(1)
i=1
where parameters τ1 and τ2 are non-negative, such that (τ1 , τ2 ) ∈√ R2+ and pℓ and qi are real polynomials in complex variable α s with arbitrary order (where α ∈ N). Note that the zeros of characteristic equation (1) are in fact the poles
[email protected] 3 M.A. Nekoui is with the Faculty of Electrical and Computer Engineering, K.N.Toosi University of Technology, Tehran, Iran
[email protected] 978-1-4799-0176-0/$31.00 ©2013 AACC
m √ √ −ℓτ1 s α s) e + p ( ∑ qi ( α s) e−iτ2 s ∑ ℓ
170
of the system under investigation. We find out from [10] that the transfer function of a system with a characteristic equation in the form of (1) will be H∞ stable if, and only if, it doesn’t 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. By applying this auxiliary variable in characteristic equation (1), the following relation is obtained: n
Cv (v, τ1 , τ2 ) =
α
∑ pℓ (v) e−ℓτ1 v
ℓ=0
m
α
+ ∑ qi (v) e−iτ2 v
A. ACFS Methodology If there exists an imaginary root of equation (1) at s = ± jωc (c for crossing) for a given set of time delays {τ } = (τ1 , τ2 ) the same imaginary root will also exist at all the countably infinite grid points of 2π 2π l, τ2 + k) ωc ωc l = 0, 1, 2, ... , k = 0, 1, 2, ... (6) 2π 2π ≤ 0, τ20 − ≤0 τ10 − ωc ωc This signifies that τi0 is the smallest positive τil , τi0 = min(τil ), i = 1, 2, l = 0, 1, 2, ..., (τil > 0). Notice that for a fixed ωc the distinct points of (6) generate a grid in {τ } ∈ R2+ space with equal grid size, 2π /ωc in both dimensions. Definition 1 (“kernel curves”): Assume that the set of (τ10 , τ20 )|ωc is determined exhaustively in {τ } = (τ1 , τ2 ) space for all possible ωc values satisfying (1) and (6).These curves as a group are called the “kernel curves” of system described by the characteristic equation (1). We denote these curves by ℘0 (τ1 , τ2 ). Definition 2 (“offspring curves”): The trajectories of (τ1 , τ2 ) grid points in (6) for l = 1, 2, ..., k = 1, 2, ... corresponding to the kernel are called the ”offspring curves” or “offspring” in short. They are represented by ℘lk (τ1 , τ2 ) where l and k identify the l th and kth generation offspring of the kernel τ10 and τ20 , respectively, according to equation (6). Let denote the complete set of kernel and offspring by {τ } = (τ1l , τ2k ) = (τ1 +
(2)
i=1
This will transform the domain of the system from a multisheeted Riemann surface into the complex plane, where the poles can be easier calculated. 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| ≤ (3) 2α with v ∈ C , which the stable region has been displayed by hatched lines in Fig. 1.
v
n
n
℘(τ1 , τ2 ) = ℘0 (τ1 , τ2 ) + ∑ + ∑ ℘lk (τ1 , τ2 )
v
(7)
l=0 k=0 l+k>0
Fig. 1.
The kernel and the offspring constitute the complete (and exhaustive) distribution of (τ1 , τ2 ) points for which the characteristic equation C(s, τ1 , τ2 ) has root sets with at least one imaginary pair. Outside the set of curves ℘(τ1 , τ2 ) there cannot be a point, which results in an imaginary characteristic root of (1). Thus in mathematical formalism, the complete imaginary root set of (1)
The v-stability region for fractional systems
Note that under this transformation, the imaginary axis in the s-domain is mapped into the lines π (4) ∠v = ± 2α
Ω[s|C(s, τ1 , τ2 ) = 0, (τ1 , τ2 ) ∈ R2+ ] ∩C0
≡ Ω[s|C(s, τ1 , τ2 ) = 0, (τ1 , τ2 ) ∈ ℘(τ1 , τ2 )] ∩C0
s = ω e± jπ /2
Let us assume that ; s = ± jω or in other words, are the roots of characteristic equation (1) for a (τ1 , τ2 ) ∈ R2+ . Then for the auxiliary variable, the roots are defined as follows: √ √ (5) v = α s = α ω e j±π /2α √ Therefore, with the auxiliary variable v = α s , there is a direct relation between the roots on the imaginary axis for the s-domain with the ones having argument ±π /2α in the v-domain.
(8)
is generated only by a small number of contours in (τ1 , τ2 ) ∈ R2+ . These are the only locations in the (τ1 , τ2 ) space where the system (1) could transit from stable to unstable posture (or vice versa). These contours ℘(τ1 , τ2 ) must be determined exhaustively. Since ℘(τ1 , τ2 ) is completely generated from the kernel via equation (6), it is sufficient to determine the kernel itself exhaustively. The determination of the complete set of kernel and offspring is, mathematically speaking, a very challenging problem. In order to achieve this we deploy a transformation called the Rekasius substitution [11]. 1 − Ti s , Ti ∈ R, i = 1, 2 (9) e −τ i s = 1 + Ti s It is important to note that, this substitution is an exact expression of e−τi s for purely imaginary roots s = ± jω .
III. METHODOLOGY In this section, we propose a method that can extract the 2-D PSSC for fractional-order systems with multiple time delay. In other words, we perform the stability analysis of (1) in 2-D delay parameter space. 171
and
Moreover, transformation (9) is different from the first-order Pade’ approximation of e−τi s ≈ (1 − 0.5τi s)/(1 + 0.5τi s). Transformation (9) can also be written as: α
e−τi v =
1 − Ti vα , 1 + Ti vα
Ti ∈ R,
i = 1, 2
2 [tan−1 (ωc Ti ) + rπ ], ωc
r = 0, 1, ...
Note that all ai and bi are real polynomials in T1 . hℜ and hℑ , which have positive degrees in terms of T2 , are assumed to have no common factors. Such common factors, if they exist, can be separately studied. Definition 3: The resultant RT2 with respect to ω and T1 is the resultant of hℜ and hℑ by eliminating T2 [12]. The resultant of RT2 and ∂ RT2 /∂ T1 with respect to ω is called the discriminant of RT2 by eliminating T1 [13],[14]. Theorem 2: Minimum and maximum positive real roots of the discriminant of ressultant of hℜ and hℑ with respect to ω , that correspond to (T1 , T2 ) ∈ R2 solutions in (14) yield the exact lower and upper bounds of the crossing frequency set. Proof: For the delay-dependent case, finite lower bound ω and upper bound ω of Ω are known to exist [2]. To find the global maximum ω and the global minimum ω , we start studying the extrema of ω via ∂ ω /∂ τ1 = 0, which is identical to studying ∂ω ∂ ω ∂ T1 = =0 (17) ∂ τ1 ∂ T1 ∂ τ1
(10)
(11)
This equation describes an asymmetric mapping in which one Ti is mapped into countably infinite τi which are distributed with a periodicity of 2π /ωc . √ Since (10) is exact when v = α ωc e jπ /2α , it is convenient to use it for solving s = va roots of (2). By inserting equation (10) into equation (2), we have: 1 − T v α i 1 − T v α ℓ m 2 1 + q (v) ∑ i 1 + T2 vα α 1 + T v 1 i=1 ℓ=0 (12) where T1 ∈ R and T2 ∈ R. By multiplying equation (12) by (1 + T1 vα )n (1 + T2 vα )m the polynomial form of the characteristic equation is reached: n
Cv (v, T1 , T2 ) =
∑ pℓ (v)
Where ∂ T1 /∂ τ1 = 0.5(1 + ω 2 T12 ) 6= 0 as per (11). Since ∂ T1 /∂ τ1 6= 0, we can study ∂ ω /∂ τ1 = 0 alternatively on ∂ ω /∂ T1 At this point, the differentiability of ω with respect to T1 is essential as established above, and holds for the regular points of RT2 (hℜ , hℑ ). Under this condition, we can write ∂ RT2 (hℜ , hℑ ) ∂ ω ∂ RT2 (hℜ , hℑ ) =0 (18) + ∂ T1 ∂ T1 ∂ω
(1 + T1 vα )n (1 + T2 vα )mCv (v, T1 , T2 ) = h(v, T1 , T2 ) n
=
∑ pℓ (v) (1 − T1 vα )ℓ (1 + T1 vα )n−ℓ (1 + T2 vα )m
ℓ=0 m
(13)
+ ∑ qi (v) (1 − T2 vα )i (1 + T1 vα )n (1 + T2 vα )m−i i=1
This expression is a polynomial in v of which the coefficients are parametric functions of T1 and T2 . As is observed, characteristic equation (2), which had transcendental terms, has been converted into algebraic equation (13). Theorem 1: A systems with charactristic equation (1) has finite crossing points for all {τ }. Proof: Let assume s √= jωc be a pair of roots √ for C(s, τ1 , τ2 ) then v = α s = α ωc e± jπ /2α would be a pair of roots for h(v, T1 , T2 ). since h(v, T1 , T2 ) is a finite degree polynomial with maximum degree of max(deg(pℓ (v)), deg(qi (v))) + α (m + n) then the number of crossing points of (1) that are the real roots of h(v, T1 , T2 ) is finite. To find the crossing frequencies set in equation (1), relation (5) should be inserted into relation (13) and then the real and imaginary parts of the resulting equation should be separated as follows: h(v, T1 , T2 )|v= α√ω
Which leads to ∂ RT2 (hℜ , hℑ )/∂ T1 = 0, assuming RT2 (hℜ , hℑ )/∂ ω 6= 0, Thus, one has two equations, RT2 and ∂ RT2 /∂ T1 , and they should be simultaneously zero. This requires to study the zeros of the resultant of these two equations, particularly by eliminating T1 . The resultant of RT2 and ∂ RT2 /∂ T1 becomes only a function of ω . Z(ω ) = RT1 (RT2 , ∂ RT2 /∂ T1 )
jπ /2α
m
hℜ = ∑ ai (ω , T1 )T2i = 0
(19)
Which is the discriminant by Definition 3. The minimum and maximum positive real zeros of Z(ω ) that correspond to (T1 , T2 ) ∈ R2 solutions in (14) are the exact lower and upper bounds of the crossing frequency set, respectively. In the sequel, the methodology ACFS is presented step by step. Notice that ACFS methodology only requires frequency sweeping from the precise lower bound ω to the precise upper bound ω that ACFS identifies via Theorem 2. For each ω ∈ [ω , ω ] with an appropriately chosen step size, perform the following steps: 1) Solve the polynomial equation RT2 (hℜ , hℑ ) for T1 ∈ R values. 2) For each T1 ∈ R found from above, if T2 ∈ R values exist satisfying hℜ = 0 and hℑ = 0 then proceed to the next step, otherwise increase ω by an amount of the step size, and restart from the step above.
= hℜ (ω , T1 , T2 ) + jhℑ (ω , T1 , T2 ) (14) In the above relations, hℜ = ℜ(h), hℑ = ℑ(h). For ω to be a zero of (14), hℜ and hℑ should be concurrently zero for some (T1 , T2 ). Let us investigate those (T1 , T2 ) solutions from ce
(16)
i=0
By examining the amplitude and phase of equation (8), the relationship between Ti and τi will be given next.
τi =
m
hℑ = ∑ bi (ω , T1 )T2i = 0
(15)
i=0
172
α
Since e−ℓτ1 v |v=vc , τ = τ10 + (2π /ωc )l, l = 0, 1, 2, ... expressions remain the same for all the l values, we can state that α e−τ1 v |v=vc is dependent only on τ10 and independent of l, or the actual value τ1 itself. Therefore,the root tendency is invariant for all τ1l , l = 0, 1, 2, .... In other words,the imaginary root always crosses either to C+ (for RT = +1) or to C− (for RT = −1), when one of the delays is kept fixed, independent of the actual values of the second delay. This theorem helps identifying certain sections of ℘0 and the offspring curves to be marked as stabilizing transitions along the τ1 (or τ2 ) axis or vice versa. Once we detect completely ℘0 (τ1 , τ2 ), ℘lk (τ1 , τ2 ) offspring curves and the invariant root sensitivities, we can determine all possible stability regions in the parametric space of time delays {τ } using the well-known D-Subdivision methodology [15]. This implies the exhaustiveness of our methodology because it covers the complete set of stability regions in the entire semiinfinite time delay space entirely. Over example cases we will also show the exactness of the methodology; determining the precise boundaries of time delays where stability switches (i.e., from stable to unstable or vice versa) occur.
3) Via (6) and (11), calculate the delay values (τ1 , τ2 ) corresponding to (T1 , T2 ) ∈ R2 pairs, and restart from step 1 increasing ω by an amount of the step size. B. Direction of Crossing The root sensitivities associated with each purely imaginary characteristic root crossing jω with respect to one of the time delay, τi , is defined as ∂ C/∂ τi ds = − i = 1, 2 (20) Sτsi s= jω = ∂ C/∂ s c d τi s= jωc s= jωc
And the corresponding root tendency with respect to one of the delays is given as: s i Root Tendency = RT |τs= jωc = sgn ℜ Sτi s= jωc (21) π π dv dv = sgn sin( )ℜ − cos( )ℑ v=vc 2α dτ 2α dτ τ =τi
In case the result is 0, a higher order analysis is needed, since this might be the case where the root just touches the imaginary axis and returns to its original half-plane. Theorem 3: Take an imaginary root, s = jωc (or v = √ α ω e jπ /2α ) caused by any one of the infinitely many grid points in (τ1 , τ2 ) defined by (τ10 + 2π l/ωc , τ20 + 2π k/ωc ), 1 l = 0, 1, 2, ..., k = 0, 1, 2, ... .The root tendency RT |τs= jωc (or τ2 RT |s= jωc ) remains invariant so long as the grid points on different offspring curves. are selected keeping τ2k (or τ1l ) fixed. Proof: characteristic equation (1) can be written in terms of the auxiliary parameter v as (2) , the simple roots of (2) are continuously differentiable with respect to {τ } ∈ R2+ . Then one can find dv/d τ for simple roots of (2) as follows: ∂ Cv/∂ τ1 dv = − ∂C v/∂ v d τ1
IV. ILLUSTRATIVE EXAMPLE We present three two cases, which display all the features discussed in the text. Example 1: Consider a multiple time delayed fractional order system with the transcendental characteristic function: √ √ √ C1 (s, τ1 , τ2 ) = ( s)3 − 1.5( s)2 + 4( s) + 8 (25) √ √ − 1.3( s)2 e−τ1 s − 0.2( s)2 e−τ2 s C1 (s, τ1 , τ2 ) is equivalent to C2 (s, τ ) when τ1 = τ2 = τ . √ √ √ √ C2 (s, τ ) = ( s)3 − 1.5( s)2 + 4( s) + 8 − 1.5( s)2 e−τ s (24) C2 (s, τ ) has a pair of poles (s = ±8 j) on the imaginary axis for τ = 0. A very involved calculation scheme based on Cauchy’s integral has been used in [16] to show that this system is unstable for τ = 0.99 and stable for τ = 1. Also it is shown in [17], [18] that the system C2 (s, τ ) is inputoutput stable for 0.0499 < τ < 0.7854 , 0.9983 < τ < 1.5708 , 1.9486 < τ < 2.3562 , 2.8953 < τ < 3.1416 and 3.8437 < τ < 3.9270. Our objective in this example is to find all the stability map for C2 (s, τ1 , τ2 ) based on the method described in this article. √ Using auxiliary variable v = s into (25) the characteristic equation is obtained as:
(22)
Since τ1 and τ2 are interchangeable we prove the theorem for τ1 and claim that it holds for τ2 also. First, we keep τ2 fixed and look at the root tendency given by (21) with respect to τ1 at grid points corresponding to ωc , i.e., (τ10 + 2π l/ωc , l = 0, 1, 2, .... It is given as h π π dv i τ1 RT |s= jωc = sgn ℜ sin( ) + jcos( ) 2α 2α d τ1 v=vc τ =τ1
= sgn ℜ −
∂ Cv/∂ v
sin( 2πα ) + jcos( 2πα ) α
∂ Cv/∂ τ1
−1 v=vc
τ =τ1 α
P (v) qi (v) α −1 q (v))e−iτ2 v 2 2 ∑nℓ=0 ℓdv e−ℓτ1 v + ∑m i i=0 ( dv − iτ2 v Cv1 (v, τ1 , τ2 ) = v3 −1.5v2 +4v+8−1.3v2 e−τ1 v −0.2v2 e−τ2 v α π π n v −ℓ α τ 1 [sin( 2α ) + jcos( 2α )] ∑ℓ=0 Pℓ (v) ℓ v e (26) −1 ατ1 v=vc By applying the criterion expressed in the previous section, − τ =τ1 (sin( 2πα ) + jcos( 2πα ))v we can eliminate exponential term from (26) as follows: Pℓ (v) −ℓτ1 vα n ∑ℓ=0 dv e h(v, T1 , T2 ) = (v3 − 1.5v2 + 4v + 8)(1 + T1 v2 )(1 + T2 v2 ) = sgn ℜ α [sin( 2πα ) + jcos( 2πα )] ∑nℓ=0 Pℓ (v) ℓ vα e−ℓτ1 v − 1.3v2 (1 − T1 v2 )(1 + T2 v2 ) − 0.2v2 (1 + T1 v2 )(1 − T2 v2 ) qi (v) α −1 q (v))e−iτ vα − i ( τ v ∑m (27) i 2 i=0 dv v=vc pω √ + α j π4 = π π n v α τ −ℓ 1 τ = τ ω e (1+ j) = a(1+ j) By inserting expression v = [sin( 2α ) + jcos( 2α )] ∑ℓ=0 Pℓ (v) ℓ v e 1 2 (23) in the above equation and equating the real and imaginary
= sgn ℜ
173
parts of the obtained relation to zero, we get: hℜ =[(−8a7 + 16a5 + 32a4 )T1 + (4a5 − 10.4a4 + 8a3 )]T2 + (4a5 − 1.6a4 + 8a3 )T1 + 2a3 − 4a − 8 = 0
(28)
and hℑ =[(8a7 + 16a5 )T1 + (4a5 − 8a3 − 16a2 )]T2
+ (4a5 − 8a3 − 16a2 )T1 − 2a3 + 6a2 − 4a = 0
(29)
Using homomorphism resultant algorithm [12] eliminate T2 from hℜ and hℑ , RT2 (hℜ , hℑ ) = (a12 − 0.2a11 − 4.4a9 + 4a8 + 8a7 + 8a6 )T12 + 1.3(a9 − 2a7 − 4a6 )T1 + 0.25a8 − 0.7a7 + 0.975a6 − 2.4a5 + a4 + 2a3 + 2a2 = 0
Fig. 2.
(30) Discriminant of the resultant of hℜ and hℑ in Theorem 2 is the resultant of RT2 and ∂ RT2 /∂ T1 with eliminating T1 ,
5
∂ RT2 ) = a14 (a6 − 0.2a5 − 4.4a3 + 4a2 + 8a + 8)× RT1 (RT2 , ∂ T1 (a12 − 3a11 + 2.77a10 − 14.78a9 + 29a8 + 0.36a7 + 43.08a6 4
3
8
4.5 7.8 4 7.6
3.5
2
− 75.84a − 91.84a − 48a + 128a + 128a + 64) = 0
3 7.4
(31)
τ2
5
Stability map of the system in example 1
The real solutions of (31) for a is
2.5 7.2
2
a = 1.82 , a=2 ,
ω = 6.6248 ω =8
(32)
1.5
7
1
p where a = ω /2 . The minimum and maximum positive real roots of Z(ω ) = RT1 (RT2 , ∂ RT2 /∂ T1 ) are computed as 6.6248 and 8 . From Theorem 2, it is concluded that [ω , ω ] is [6.6248, 8]. One can now use this ω range and the ACFS method in previous section, in order to extract the stability maps on τ1 − τ2 domain by sweeping the frequency from 6.6248 to 8. The potential stability switching curves of the system is extracted in Fig. 2, where the kernel curve is shown in red, and the offspring curves are given in blue when viewed in color. The shaded zone in Fig.2 show the complete map of stability for the given system. The number of unstable roots in each region, NU, is also shown sparingly. Obviously, in the stable regions, NU = 0. To get a better understanding of the properties of this system, the kernel and offspring curves have been plotted as a function of frequency distribution in Fig. 3. The color spectrum in the “color bar” indicates the selected ω ; dark blue designates ω = 0.6.6248, and dark red is for ω = 8. Example 2: Consider the following fractional-order system with two time delays [16],[18]: √ √ √ (33) C(s) = ( 6 s)5 + ( 6 s)3 + ( 6 s)2 e−0.5s + e−s
6.8
0.5
6.6
0 0
1
2
3
4
5
ω
τ1
Fig. 3.
Color-coded ωc of the kernel and offspring curves in example 1
Now we evaluate the stability of this system for all the values of (τ1 , τ2 ), and also study the stability for τ1 = 0.5 and τ2 = 1 This system has no unstable pole for (τ1 , τ2 ) = (0, 0) (in fact, it has no pole in the physical Riemann sheet). √ By using of auxiliary variable v = 6 s in (34), we will have: 6
6
Cv2 (v, τ1 , τ2 ) = v5 + (v3 + v2 )e−τ1 v + e−τ2 v
(35)
By eliminating exponential term in (35) , we have the following equation: h(v, T1 , T2 ) = v5 (1 + T1 v6 )(1 + T2 v6 ) + (v3 + v2 )(1 − T1 v6 )(1 + T2 v6 ) + (1 + T1 v6 )(1 − T2 v6 ) (36) √ jπ √ And , since v = 6 s = 6 ω e 12 we will have: √ √ √ √ v = a(( 6 + 2) + j( 6 − 2)) (37) √ where a = 6 ω /4. By inserting equation (37) into equation (36) and making the real and imaginary parts of h(v, T1 , T2 ) to zero, and Following Theorem 2, corresponding frequency
It has been demonstrated in [16] that this system is stable. In order to apply our procedure, we transformed C(s) into the following form. √ √ √ C3 (s, τ1 , τ2 ) = ( 6 s)5 + ( 6 s)3 + ( 6 s)2 e−τ1 s + e−τ2 s (34) 174
ranges are found as: [ω , ω ] = [0.1191, 8.3958]
for studying the asymptotic stability of multiple time-delay fractional-order systems in the parameter space of delays. It is evident from the literature that the stability assessment of this class of dynamics remains unsolved. The method √ commenced by deploying the auxiliary variable v = α s and Rekasius substitution for the transcendental terms in the characteristic equation, reducing it into a finite dimensional algebraic equation. Moreover, a single-variable function is derived to precisely calculate the lower and upper bounds of the crossing frequency set. These bounds are critical to the ACFS implementation. By means of ACFS, potential stability switching curves (PSSC) on any 2-D delay domain are extracted exhaustively. The proposed methodology detects all the stability regions precisely, in the space of the time delays. This map, in fact, is the exact display of robustness against delay uncertainties. Two examples presented to highlight the proposed approach.
(38)
Upon sweeping ω in this range, the potential stability switching curves of the system is extracted in Fig.4. The kernel curve is shown in gray and the offspring curves are given in black (red and blue, respectively, when viewed in color) and the shaded regions are stable operating points (NU = 0) and distribution of the number of unstable roots, NU, is also marked sparingly. As (τ1 , τ2 ) = (0.5, 1) is included in the shaded regions , we can ensure that the original system C(s) is stable. In Fig. 5, the kernel and offspring curves of this system for the changes of ω from ω = 0.1191 (dark blue) to ω = 8.3958 (dark red) has been presented for a better perception of the system.
R EFERENCES
Fig. 4.
Stability map of the system in example 2
5
8.3958 8
4.5 7
τ2
4 3.5
6
3
5
2.5 4 2 3 1.5 2
1
1
0.5
0.1191
0 0
0.5
1
1.5
2
2.5
ω
τ1
Fig. 5.
Colo-coded ωc of the kernel and offspring curves in example 2
V. CONCLUSIONS In this paper, a novel methodology, based on Advanced Clustering with Frequency Sweeping (ACFS), is proposed 175
[1] I.I., Delice, R. Sipahi, Advanced Clustering with Frequency Sweeping (ACFS) Methodology for the Stability Analysis of Multiple Time Delay Systems, American Control Conference, Baltimore, 50125017,2010. [2] R. Sipahi, I.I., Delice, Advanced Clustering with Frequency Sweeping Methodology for the Stability Analysis of Multiple Time-Delay Systems, IEEE Trans. on Automatic Control, vol.56, no.2, 467-472, 2011 [3] K. Gu, S.-I. Niculescu, and J. Chen, On stability crossing curves for general systems with two delays, J. Math. Anal. Appl., vol. 311, no. 1, pp. 231-253, 2005. [4] R. Sipahi and N. Olgac, Complete stability robustness of third-order LTI multiple time-delay systems, Automatica, vol. 41, no. 8, pp. 14131422, 2005. [5] H. Fazelinia, R. Sipahi, and N. Olgac, Stability analysis of multiple time delayed systems using ‘building block’ concept, IEEE Trans. Autom. Control, vol. 52, no. 5, pp. 799-810, May 2007. [6] S. Pakzad, and M. A. Pakzad, Stability condition for discrete systems with multiple state delays, WSEAS Trans. on Systems and Control, vol.6, no.11, pp. 417-426, 2011. ¨ urk and A. Uraz, An analysis stability test for a certain class [7] N. Ozt¨ of distributed parameter systems with delays, IEEE Trans. on Circuits and Systems, vol. 34, no. 4, pp. 393-396, 1985. [8] E.I. Jury and E. Zeheb, On a stability test for a class of distributed system with delays, IEEE Trans. on Circuits and Systems, vol. 37, no.10, pp. 1027-1028, 1986. [9] C. Bonnet and J.R. Partington, Analysis of fractional delay systems of retarded and neutral type, Automatica, vol. 38, no. 7, pp. 1133-1138, 2002. [10] C. Bonnet and J.R. Partington, Stabilization of some fractional delay systems of neutral type, Automatica, vol. 43, no. 12, pp. 2047-2053, 2007. [11] Z. V. Rekasius, A stability test for systems with delays, in Proc. Joint Autom. Control Conf., San Francisco, CA, 1980. [12] G. E. Collins, The calculation of multivariate polynomial resultants, J. Assoc. Comput. Mach., vol. 18, no. 4, pp. 515532, 1971. [13] I. M. Gelfand, M. M. Kapranov, and A. V. Zelevinsky, Discriminants, Resultants, and Multidimensional Determinants. Boston, MA: Birkhuser, 1994. [14] D. Wang, Elimination Methods. New York: Springer-Verlag, 2001. [15] V.B. Kolmanovskii, V.R. Nosov, Stability of Functional Differential Equations, Academic Press, London, 1989. [16] C. Hwang and Y.C. Cheng, A numerical algorithm for stability testing of fractional delay systems, Automatica, vol. 42, no. 5, pp. 825-831, 2006. [17] M. A. Pakzad, S. Pakzad, Stability map of fractional order time-delay systems, WSEAS Trans. on Systems, vol.10, no.11, pp. 541-550, 2012. [18] M. A. Pakzad, M.A. Nekoui and S. Pakzad, Stability analysis of time-delayed linear fractional-order systems, International Journal of Control, Automation and Systems, vol.11, no.3, 2013, doi: 10.1007/s12555-012-0164-4.
Stability Analysis of Multiple Time Delayed Fractional Order Systems Mohammad Ali Pakzad1 , Sara Pakzad2 , and Mohammad Ali Nekoui3
Obviously, the accuracy and completeness of the analysis strongly depends on finding all the existing PSSC without any approximations [1]. There has been a large effort to deal with this problem, for the standard case (integer order systems); see [1]-[6], and others. The researchers of [7] and [8] are probably the pioneer to consider stability analysis of the fractional delay system with single-delay. They have extended the Ruth-Hurwitz criteria for analyzing the stability of some √ special delay systems to those involve fractional power s. For single delay case; In [9], necessary and sufficient conditions for BIBO stability of the retarded fractional-order delay systems and sufficient conditions for some neutral types have been introduced; and also from the analitical analysis point of view, the effective analitical approaches to check BIBO stability of fractionalorder delay systems have been discussed in [17] and [18]. The original idea in this strategy is derived based on the method (advanced clustering with frequency sweeping) reported in Reference [1] to achieve the stability map of integer order systems with multiple time delays. In this paper, we extend the approach of [1] to fractionalorder systems with multiple delays. Necessary and sufficient conditions which yield the exact lower and upper bounds of the crossing frequency set (CFS) can be computed via an automated sequential formula. These bounds are crucial as they determine the sweeping range of the only parameter, the frequency that ACFS sweeps.
Abstract— A new methodology, based on advanced clustering with frequency sweeping (ACFS), is presented for the stability analysis of fractional-order systems with multiple time delays against delay uncertainties. The problem is known to be notoriously complex, primarily because the systems are infinite dimensional due to delays. Multiplicity of the delays in this study complicates the analysis even further. And “fractionalorder” feature of the systems makes the problem much more challenging compared to integer order systems. ACFS does not impose any restrictions in the number of delays, and it can directly extract the 2-D cross sections of the stability views in any two delay domain. We show that this procedure analytically reveals all possible stability regions exclusively in the space of the delays. As an added strength, it does not require the delay-free system under consideration to be stable. The main contribution of this document is that we demonstrate for the first time that the stability maps of a fractional-order system with multiple time delays can be obtained efficiently. Two illustrative examples are presented to confirm the proposed method results.
I. INTRODUCTION In a generic sense, the fractional-order systems are identified by non-integer powers of the Laplace variable s. When time delays and fractional-order derivative are involved in dynamical systems, we have fractional-delay systems. The characteristic equation of a fractional-delay system involves exponential type transcendental terms, so a fractional delay system has infinitely many roots. This makes the stability analysis of fractional-delay systems a challenging task. In this paper, one of the most important and unresolved problems of fractional delay systems is studied: the asymptotic stability of a general class of fractional-order systems with multiple time delays, against delay uncertainties. We make two qualitative comments first. (a) There is no methodology presently to the knowledge of the authors, which can reveal the aimed “stability map” of this system. (b) As a natural consequence of (a) there is no methodology to present this stability map exhaustively. Exhaustiveness means including all the stable regions in a given segment of interest in (τ1 , τ2 ) space attached to or detached from the origin (0, 0). The main objective in 2-D stability analysis is to construct all the potential stability switching curves (PSSC) which partition the delay space into stable and unstable regions.
II. PRELIMINARIES AND DEFINITIONS Standard notation has been used throughout the article: N , Z : sets of natural and integer numbers. R(R+ , R− ) : set of real (positive real, negative real) numbers. C:√ set of complex numbers. j = −1 : the imaginary unit. z : complex conjugate of a complex number z ℜ(z) : real part of a complex number z ℑ(z) : imaginary part of a complex number z |z|, ∠z:magnitude and argument of a complex number z RTℓ (p1 , p2 ): denotes the resultant of bivariate polynomials p1 (T1 , T2 ) and p2 (T1 , T2 ) with eliminating Tℓ , where ℓ = 1, 2 Consider a fractional-order system with the following characteristic equation:
1 M.A. Pakzad is with Department of Electrical Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
n
C(s, τ1 , τ2 ) =
[email protected] 2 S.
South
ℓ=0
Pakzad is with the Departments of Electrical Engineering, Tehran Branch, Islamic Azad University, Tehran, Iran
(1)
i=1
where parameters τ1 and τ2 are non-negative, such that (τ1 , τ2 ) ∈√ R2+ and pℓ and qi are real polynomials in complex variable α s with arbitrary order (where α ∈ N). Note that the zeros of characteristic equation (1) are in fact the poles
[email protected] 3 M.A. Nekoui is with the Faculty of Electrical and Computer Engineering, K.N.Toosi University of Technology, Tehran, Iran
[email protected] 978-1-4799-0176-0/$31.00 ©2013 AACC
m √ √ −ℓτ1 s α s) e + p ( ∑ qi ( α s) e−iτ2 s ∑ ℓ
170
of the system under investigation. We find out from [10] that the transfer function of a system with a characteristic equation in the form of (1) will be H∞ stable if, and only if, it doesn’t 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. By applying this auxiliary variable in characteristic equation (1), the following relation is obtained: n
Cv (v, τ1 , τ2 ) =
α
∑ pℓ (v) e−ℓτ1 v
ℓ=0
m
α
+ ∑ qi (v) e−iτ2 v
A. ACFS Methodology If there exists an imaginary root of equation (1) at s = ± jωc (c for crossing) for a given set of time delays {τ } = (τ1 , τ2 ) the same imaginary root will also exist at all the countably infinite grid points of 2π 2π l, τ2 + k) ωc ωc l = 0, 1, 2, ... , k = 0, 1, 2, ... (6) 2π 2π ≤ 0, τ20 − ≤0 τ10 − ωc ωc This signifies that τi0 is the smallest positive τil , τi0 = min(τil ), i = 1, 2, l = 0, 1, 2, ..., (τil > 0). Notice that for a fixed ωc the distinct points of (6) generate a grid in {τ } ∈ R2+ space with equal grid size, 2π /ωc in both dimensions. Definition 1 (“kernel curves”): Assume that the set of (τ10 , τ20 )|ωc is determined exhaustively in {τ } = (τ1 , τ2 ) space for all possible ωc values satisfying (1) and (6).These curves as a group are called the “kernel curves” of system described by the characteristic equation (1). We denote these curves by ℘0 (τ1 , τ2 ). Definition 2 (“offspring curves”): The trajectories of (τ1 , τ2 ) grid points in (6) for l = 1, 2, ..., k = 1, 2, ... corresponding to the kernel are called the ”offspring curves” or “offspring” in short. They are represented by ℘lk (τ1 , τ2 ) where l and k identify the l th and kth generation offspring of the kernel τ10 and τ20 , respectively, according to equation (6). Let denote the complete set of kernel and offspring by {τ } = (τ1l , τ2k ) = (τ1 +
(2)
i=1
This will transform the domain of the system from a multisheeted Riemann surface into the complex plane, where the poles can be easier calculated. 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| ≤ (3) 2α with v ∈ C , which the stable region has been displayed by hatched lines in Fig. 1.
v
n
n
℘(τ1 , τ2 ) = ℘0 (τ1 , τ2 ) + ∑ + ∑ ℘lk (τ1 , τ2 )
v
(7)
l=0 k=0 l+k>0
Fig. 1.
The kernel and the offspring constitute the complete (and exhaustive) distribution of (τ1 , τ2 ) points for which the characteristic equation C(s, τ1 , τ2 ) has root sets with at least one imaginary pair. Outside the set of curves ℘(τ1 , τ2 ) there cannot be a point, which results in an imaginary characteristic root of (1). Thus in mathematical formalism, the complete imaginary root set of (1)
The v-stability region for fractional systems
Note that under this transformation, the imaginary axis in the s-domain is mapped into the lines π (4) ∠v = ± 2α
Ω[s|C(s, τ1 , τ2 ) = 0, (τ1 , τ2 ) ∈ R2+ ] ∩C0
≡ Ω[s|C(s, τ1 , τ2 ) = 0, (τ1 , τ2 ) ∈ ℘(τ1 , τ2 )] ∩C0
s = ω e± jπ /2
Let us assume that ; s = ± jω or in other words, are the roots of characteristic equation (1) for a (τ1 , τ2 ) ∈ R2+ . Then for the auxiliary variable, the roots are defined as follows: √ √ (5) v = α s = α ω e j±π /2α √ Therefore, with the auxiliary variable v = α s , there is a direct relation between the roots on the imaginary axis for the s-domain with the ones having argument ±π /2α in the v-domain.
(8)
is generated only by a small number of contours in (τ1 , τ2 ) ∈ R2+ . These are the only locations in the (τ1 , τ2 ) space where the system (1) could transit from stable to unstable posture (or vice versa). These contours ℘(τ1 , τ2 ) must be determined exhaustively. Since ℘(τ1 , τ2 ) is completely generated from the kernel via equation (6), it is sufficient to determine the kernel itself exhaustively. The determination of the complete set of kernel and offspring is, mathematically speaking, a very challenging problem. In order to achieve this we deploy a transformation called the Rekasius substitution [11]. 1 − Ti s , Ti ∈ R, i = 1, 2 (9) e −τ i s = 1 + Ti s It is important to note that, this substitution is an exact expression of e−τi s for purely imaginary roots s = ± jω .
III. METHODOLOGY In this section, we propose a method that can extract the 2-D PSSC for fractional-order systems with multiple time delay. In other words, we perform the stability analysis of (1) in 2-D delay parameter space. 171
and
Moreover, transformation (9) is different from the first-order Pade’ approximation of e−τi s ≈ (1 − 0.5τi s)/(1 + 0.5τi s). Transformation (9) can also be written as: α
e−τi v =
1 − Ti vα , 1 + Ti vα
Ti ∈ R,
i = 1, 2
2 [tan−1 (ωc Ti ) + rπ ], ωc
r = 0, 1, ...
Note that all ai and bi are real polynomials in T1 . hℜ and hℑ , which have positive degrees in terms of T2 , are assumed to have no common factors. Such common factors, if they exist, can be separately studied. Definition 3: The resultant RT2 with respect to ω and T1 is the resultant of hℜ and hℑ by eliminating T2 [12]. The resultant of RT2 and ∂ RT2 /∂ T1 with respect to ω is called the discriminant of RT2 by eliminating T1 [13],[14]. Theorem 2: Minimum and maximum positive real roots of the discriminant of ressultant of hℜ and hℑ with respect to ω , that correspond to (T1 , T2 ) ∈ R2 solutions in (14) yield the exact lower and upper bounds of the crossing frequency set. Proof: For the delay-dependent case, finite lower bound ω and upper bound ω of Ω are known to exist [2]. To find the global maximum ω and the global minimum ω , we start studying the extrema of ω via ∂ ω /∂ τ1 = 0, which is identical to studying ∂ω ∂ ω ∂ T1 = =0 (17) ∂ τ1 ∂ T1 ∂ τ1
(10)
(11)
This equation describes an asymmetric mapping in which one Ti is mapped into countably infinite τi which are distributed with a periodicity of 2π /ωc . √ Since (10) is exact when v = α ωc e jπ /2α , it is convenient to use it for solving s = va roots of (2). By inserting equation (10) into equation (2), we have: 1 − T v α i 1 − T v α ℓ m 2 1 + q (v) ∑ i 1 + T2 vα α 1 + T v 1 i=1 ℓ=0 (12) where T1 ∈ R and T2 ∈ R. By multiplying equation (12) by (1 + T1 vα )n (1 + T2 vα )m the polynomial form of the characteristic equation is reached: n
Cv (v, T1 , T2 ) =
∑ pℓ (v)
Where ∂ T1 /∂ τ1 = 0.5(1 + ω 2 T12 ) 6= 0 as per (11). Since ∂ T1 /∂ τ1 6= 0, we can study ∂ ω /∂ τ1 = 0 alternatively on ∂ ω /∂ T1 At this point, the differentiability of ω with respect to T1 is essential as established above, and holds for the regular points of RT2 (hℜ , hℑ ). Under this condition, we can write ∂ RT2 (hℜ , hℑ ) ∂ ω ∂ RT2 (hℜ , hℑ ) =0 (18) + ∂ T1 ∂ T1 ∂ω
(1 + T1 vα )n (1 + T2 vα )mCv (v, T1 , T2 ) = h(v, T1 , T2 ) n
=
∑ pℓ (v) (1 − T1 vα )ℓ (1 + T1 vα )n−ℓ (1 + T2 vα )m
ℓ=0 m
(13)
+ ∑ qi (v) (1 − T2 vα )i (1 + T1 vα )n (1 + T2 vα )m−i i=1
This expression is a polynomial in v of which the coefficients are parametric functions of T1 and T2 . As is observed, characteristic equation (2), which had transcendental terms, has been converted into algebraic equation (13). Theorem 1: A systems with charactristic equation (1) has finite crossing points for all {τ }. Proof: Let assume s √= jωc be a pair of roots √ for C(s, τ1 , τ2 ) then v = α s = α ωc e± jπ /2α would be a pair of roots for h(v, T1 , T2 ). since h(v, T1 , T2 ) is a finite degree polynomial with maximum degree of max(deg(pℓ (v)), deg(qi (v))) + α (m + n) then the number of crossing points of (1) that are the real roots of h(v, T1 , T2 ) is finite. To find the crossing frequencies set in equation (1), relation (5) should be inserted into relation (13) and then the real and imaginary parts of the resulting equation should be separated as follows: h(v, T1 , T2 )|v= α√ω
Which leads to ∂ RT2 (hℜ , hℑ )/∂ T1 = 0, assuming RT2 (hℜ , hℑ )/∂ ω 6= 0, Thus, one has two equations, RT2 and ∂ RT2 /∂ T1 , and they should be simultaneously zero. This requires to study the zeros of the resultant of these two equations, particularly by eliminating T1 . The resultant of RT2 and ∂ RT2 /∂ T1 becomes only a function of ω . Z(ω ) = RT1 (RT2 , ∂ RT2 /∂ T1 )
jπ /2α
m
hℜ = ∑ ai (ω , T1 )T2i = 0
(19)
Which is the discriminant by Definition 3. The minimum and maximum positive real zeros of Z(ω ) that correspond to (T1 , T2 ) ∈ R2 solutions in (14) are the exact lower and upper bounds of the crossing frequency set, respectively. In the sequel, the methodology ACFS is presented step by step. Notice that ACFS methodology only requires frequency sweeping from the precise lower bound ω to the precise upper bound ω that ACFS identifies via Theorem 2. For each ω ∈ [ω , ω ] with an appropriately chosen step size, perform the following steps: 1) Solve the polynomial equation RT2 (hℜ , hℑ ) for T1 ∈ R values. 2) For each T1 ∈ R found from above, if T2 ∈ R values exist satisfying hℜ = 0 and hℑ = 0 then proceed to the next step, otherwise increase ω by an amount of the step size, and restart from the step above.
= hℜ (ω , T1 , T2 ) + jhℑ (ω , T1 , T2 ) (14) In the above relations, hℜ = ℜ(h), hℑ = ℑ(h). For ω to be a zero of (14), hℜ and hℑ should be concurrently zero for some (T1 , T2 ). Let us investigate those (T1 , T2 ) solutions from ce
(16)
i=0
By examining the amplitude and phase of equation (8), the relationship between Ti and τi will be given next.
τi =
m
hℑ = ∑ bi (ω , T1 )T2i = 0
(15)
i=0
172
α
Since e−ℓτ1 v |v=vc , τ = τ10 + (2π /ωc )l, l = 0, 1, 2, ... expressions remain the same for all the l values, we can state that α e−τ1 v |v=vc is dependent only on τ10 and independent of l, or the actual value τ1 itself. Therefore,the root tendency is invariant for all τ1l , l = 0, 1, 2, .... In other words,the imaginary root always crosses either to C+ (for RT = +1) or to C− (for RT = −1), when one of the delays is kept fixed, independent of the actual values of the second delay. This theorem helps identifying certain sections of ℘0 and the offspring curves to be marked as stabilizing transitions along the τ1 (or τ2 ) axis or vice versa. Once we detect completely ℘0 (τ1 , τ2 ), ℘lk (τ1 , τ2 ) offspring curves and the invariant root sensitivities, we can determine all possible stability regions in the parametric space of time delays {τ } using the well-known D-Subdivision methodology [15]. This implies the exhaustiveness of our methodology because it covers the complete set of stability regions in the entire semiinfinite time delay space entirely. Over example cases we will also show the exactness of the methodology; determining the precise boundaries of time delays where stability switches (i.e., from stable to unstable or vice versa) occur.
3) Via (6) and (11), calculate the delay values (τ1 , τ2 ) corresponding to (T1 , T2 ) ∈ R2 pairs, and restart from step 1 increasing ω by an amount of the step size. B. Direction of Crossing The root sensitivities associated with each purely imaginary characteristic root crossing jω with respect to one of the time delay, τi , is defined as ∂ C/∂ τi ds = − i = 1, 2 (20) Sτsi s= jω = ∂ C/∂ s c d τi s= jωc s= jωc
And the corresponding root tendency with respect to one of the delays is given as: s i Root Tendency = RT |τs= jωc = sgn ℜ Sτi s= jωc (21) π π dv dv = sgn sin( )ℜ − cos( )ℑ v=vc 2α dτ 2α dτ τ =τi
In case the result is 0, a higher order analysis is needed, since this might be the case where the root just touches the imaginary axis and returns to its original half-plane. Theorem 3: Take an imaginary root, s = jωc (or v = √ α ω e jπ /2α ) caused by any one of the infinitely many grid points in (τ1 , τ2 ) defined by (τ10 + 2π l/ωc , τ20 + 2π k/ωc ), 1 l = 0, 1, 2, ..., k = 0, 1, 2, ... .The root tendency RT |τs= jωc (or τ2 RT |s= jωc ) remains invariant so long as the grid points on different offspring curves. are selected keeping τ2k (or τ1l ) fixed. Proof: characteristic equation (1) can be written in terms of the auxiliary parameter v as (2) , the simple roots of (2) are continuously differentiable with respect to {τ } ∈ R2+ . Then one can find dv/d τ for simple roots of (2) as follows: ∂ Cv/∂ τ1 dv = − ∂C v/∂ v d τ1
IV. ILLUSTRATIVE EXAMPLE We present three two cases, which display all the features discussed in the text. Example 1: Consider a multiple time delayed fractional order system with the transcendental characteristic function: √ √ √ C1 (s, τ1 , τ2 ) = ( s)3 − 1.5( s)2 + 4( s) + 8 (25) √ √ − 1.3( s)2 e−τ1 s − 0.2( s)2 e−τ2 s C1 (s, τ1 , τ2 ) is equivalent to C2 (s, τ ) when τ1 = τ2 = τ . √ √ √ √ C2 (s, τ ) = ( s)3 − 1.5( s)2 + 4( s) + 8 − 1.5( s)2 e−τ s (24) C2 (s, τ ) has a pair of poles (s = ±8 j) on the imaginary axis for τ = 0. A very involved calculation scheme based on Cauchy’s integral has been used in [16] to show that this system is unstable for τ = 0.99 and stable for τ = 1. Also it is shown in [17], [18] that the system C2 (s, τ ) is inputoutput stable for 0.0499 < τ < 0.7854 , 0.9983 < τ < 1.5708 , 1.9486 < τ < 2.3562 , 2.8953 < τ < 3.1416 and 3.8437 < τ < 3.9270. Our objective in this example is to find all the stability map for C2 (s, τ1 , τ2 ) based on the method described in this article. √ Using auxiliary variable v = s into (25) the characteristic equation is obtained as:
(22)
Since τ1 and τ2 are interchangeable we prove the theorem for τ1 and claim that it holds for τ2 also. First, we keep τ2 fixed and look at the root tendency given by (21) with respect to τ1 at grid points corresponding to ωc , i.e., (τ10 + 2π l/ωc , l = 0, 1, 2, .... It is given as h π π dv i τ1 RT |s= jωc = sgn ℜ sin( ) + jcos( ) 2α 2α d τ1 v=vc τ =τ1
= sgn ℜ −
∂ Cv/∂ v
sin( 2πα ) + jcos( 2πα ) α
∂ Cv/∂ τ1
−1 v=vc
τ =τ1 α
P (v) qi (v) α −1 q (v))e−iτ2 v 2 2 ∑nℓ=0 ℓdv e−ℓτ1 v + ∑m i i=0 ( dv − iτ2 v Cv1 (v, τ1 , τ2 ) = v3 −1.5v2 +4v+8−1.3v2 e−τ1 v −0.2v2 e−τ2 v α π π n v −ℓ α τ 1 [sin( 2α ) + jcos( 2α )] ∑ℓ=0 Pℓ (v) ℓ v e (26) −1 ατ1 v=vc By applying the criterion expressed in the previous section, − τ =τ1 (sin( 2πα ) + jcos( 2πα ))v we can eliminate exponential term from (26) as follows: Pℓ (v) −ℓτ1 vα n ∑ℓ=0 dv e h(v, T1 , T2 ) = (v3 − 1.5v2 + 4v + 8)(1 + T1 v2 )(1 + T2 v2 ) = sgn ℜ α [sin( 2πα ) + jcos( 2πα )] ∑nℓ=0 Pℓ (v) ℓ vα e−ℓτ1 v − 1.3v2 (1 − T1 v2 )(1 + T2 v2 ) − 0.2v2 (1 + T1 v2 )(1 − T2 v2 ) qi (v) α −1 q (v))e−iτ vα − i ( τ v ∑m (27) i 2 i=0 dv v=vc pω √ + α j π4 = π π n v α τ −ℓ 1 τ = τ ω e (1+ j) = a(1+ j) By inserting expression v = [sin( 2α ) + jcos( 2α )] ∑ℓ=0 Pℓ (v) ℓ v e 1 2 (23) in the above equation and equating the real and imaginary
= sgn ℜ
173
parts of the obtained relation to zero, we get: hℜ =[(−8a7 + 16a5 + 32a4 )T1 + (4a5 − 10.4a4 + 8a3 )]T2 + (4a5 − 1.6a4 + 8a3 )T1 + 2a3 − 4a − 8 = 0
(28)
and hℑ =[(8a7 + 16a5 )T1 + (4a5 − 8a3 − 16a2 )]T2
+ (4a5 − 8a3 − 16a2 )T1 − 2a3 + 6a2 − 4a = 0
(29)
Using homomorphism resultant algorithm [12] eliminate T2 from hℜ and hℑ , RT2 (hℜ , hℑ ) = (a12 − 0.2a11 − 4.4a9 + 4a8 + 8a7 + 8a6 )T12 + 1.3(a9 − 2a7 − 4a6 )T1 + 0.25a8 − 0.7a7 + 0.975a6 − 2.4a5 + a4 + 2a3 + 2a2 = 0
Fig. 2.
(30) Discriminant of the resultant of hℜ and hℑ in Theorem 2 is the resultant of RT2 and ∂ RT2 /∂ T1 with eliminating T1 ,
5
∂ RT2 ) = a14 (a6 − 0.2a5 − 4.4a3 + 4a2 + 8a + 8)× RT1 (RT2 , ∂ T1 (a12 − 3a11 + 2.77a10 − 14.78a9 + 29a8 + 0.36a7 + 43.08a6 4
3
8
4.5 7.8 4 7.6
3.5
2
− 75.84a − 91.84a − 48a + 128a + 128a + 64) = 0
3 7.4
(31)
τ2
5
Stability map of the system in example 1
The real solutions of (31) for a is
2.5 7.2
2
a = 1.82 , a=2 ,
ω = 6.6248 ω =8
(32)
1.5
7
1
p where a = ω /2 . The minimum and maximum positive real roots of Z(ω ) = RT1 (RT2 , ∂ RT2 /∂ T1 ) are computed as 6.6248 and 8 . From Theorem 2, it is concluded that [ω , ω ] is [6.6248, 8]. One can now use this ω range and the ACFS method in previous section, in order to extract the stability maps on τ1 − τ2 domain by sweeping the frequency from 6.6248 to 8. The potential stability switching curves of the system is extracted in Fig. 2, where the kernel curve is shown in red, and the offspring curves are given in blue when viewed in color. The shaded zone in Fig.2 show the complete map of stability for the given system. The number of unstable roots in each region, NU, is also shown sparingly. Obviously, in the stable regions, NU = 0. To get a better understanding of the properties of this system, the kernel and offspring curves have been plotted as a function of frequency distribution in Fig. 3. The color spectrum in the “color bar” indicates the selected ω ; dark blue designates ω = 0.6.6248, and dark red is for ω = 8. Example 2: Consider the following fractional-order system with two time delays [16],[18]: √ √ √ (33) C(s) = ( 6 s)5 + ( 6 s)3 + ( 6 s)2 e−0.5s + e−s
6.8
0.5
6.6
0 0
1
2
3
4
5
ω
τ1
Fig. 3.
Color-coded ωc of the kernel and offspring curves in example 1
Now we evaluate the stability of this system for all the values of (τ1 , τ2 ), and also study the stability for τ1 = 0.5 and τ2 = 1 This system has no unstable pole for (τ1 , τ2 ) = (0, 0) (in fact, it has no pole in the physical Riemann sheet). √ By using of auxiliary variable v = 6 s in (34), we will have: 6
6
Cv2 (v, τ1 , τ2 ) = v5 + (v3 + v2 )e−τ1 v + e−τ2 v
(35)
By eliminating exponential term in (35) , we have the following equation: h(v, T1 , T2 ) = v5 (1 + T1 v6 )(1 + T2 v6 ) + (v3 + v2 )(1 − T1 v6 )(1 + T2 v6 ) + (1 + T1 v6 )(1 − T2 v6 ) (36) √ jπ √ And , since v = 6 s = 6 ω e 12 we will have: √ √ √ √ v = a(( 6 + 2) + j( 6 − 2)) (37) √ where a = 6 ω /4. By inserting equation (37) into equation (36) and making the real and imaginary parts of h(v, T1 , T2 ) to zero, and Following Theorem 2, corresponding frequency
It has been demonstrated in [16] that this system is stable. In order to apply our procedure, we transformed C(s) into the following form. √ √ √ C3 (s, τ1 , τ2 ) = ( 6 s)5 + ( 6 s)3 + ( 6 s)2 e−τ1 s + e−τ2 s (34) 174
ranges are found as: [ω , ω ] = [0.1191, 8.3958]
for studying the asymptotic stability of multiple time-delay fractional-order systems in the parameter space of delays. It is evident from the literature that the stability assessment of this class of dynamics remains unsolved. The method √ commenced by deploying the auxiliary variable v = α s and Rekasius substitution for the transcendental terms in the characteristic equation, reducing it into a finite dimensional algebraic equation. Moreover, a single-variable function is derived to precisely calculate the lower and upper bounds of the crossing frequency set. These bounds are critical to the ACFS implementation. By means of ACFS, potential stability switching curves (PSSC) on any 2-D delay domain are extracted exhaustively. The proposed methodology detects all the stability regions precisely, in the space of the time delays. This map, in fact, is the exact display of robustness against delay uncertainties. Two examples presented to highlight the proposed approach.
(38)
Upon sweeping ω in this range, the potential stability switching curves of the system is extracted in Fig.4. The kernel curve is shown in gray and the offspring curves are given in black (red and blue, respectively, when viewed in color) and the shaded regions are stable operating points (NU = 0) and distribution of the number of unstable roots, NU, is also marked sparingly. As (τ1 , τ2 ) = (0.5, 1) is included in the shaded regions , we can ensure that the original system C(s) is stable. In Fig. 5, the kernel and offspring curves of this system for the changes of ω from ω = 0.1191 (dark blue) to ω = 8.3958 (dark red) has been presented for a better perception of the system.
R EFERENCES
Fig. 4.
Stability map of the system in example 2
5
8.3958 8
4.5 7
τ2
4 3.5
6
3
5
2.5 4 2 3 1.5 2
1
1
0.5
0.1191
0 0
0.5
1
1.5
2
2.5
ω
τ1
Fig. 5.
Colo-coded ωc of the kernel and offspring curves in example 2
V. CONCLUSIONS In this paper, a novel methodology, based on Advanced Clustering with Frequency Sweeping (ACFS), is proposed 175
[1] I.I., Delice, R. Sipahi, Advanced Clustering with Frequency Sweeping (ACFS) Methodology for the Stability Analysis of Multiple Time Delay Systems, American Control Conference, Baltimore, 50125017,2010. [2] R. Sipahi, I.I., Delice, Advanced Clustering with Frequency Sweeping Methodology for the Stability Analysis of Multiple Time-Delay Systems, IEEE Trans. on Automatic Control, vol.56, no.2, 467-472, 2011 [3] K. Gu, S.-I. Niculescu, and J. Chen, On stability crossing curves for general systems with two delays, J. Math. Anal. Appl., vol. 311, no. 1, pp. 231-253, 2005. [4] R. Sipahi and N. Olgac, Complete stability robustness of third-order LTI multiple time-delay systems, Automatica, vol. 41, no. 8, pp. 14131422, 2005. [5] H. Fazelinia, R. Sipahi, and N. Olgac, Stability analysis of multiple time delayed systems using ‘building block’ concept, IEEE Trans. Autom. Control, vol. 52, no. 5, pp. 799-810, May 2007. [6] S. Pakzad, and M. A. Pakzad, Stability condition for discrete systems with multiple state delays, WSEAS Trans. on Systems and Control, vol.6, no.11, pp. 417-426, 2011. ¨ urk and A. Uraz, An analysis stability test for a certain class [7] N. Ozt¨ of distributed parameter systems with delays, IEEE Trans. on Circuits and Systems, vol. 34, no. 4, pp. 393-396, 1985. [8] E.I. Jury and E. Zeheb, On a stability test for a class of distributed system with delays, IEEE Trans. on Circuits and Systems, vol. 37, no.10, pp. 1027-1028, 1986. [9] C. Bonnet and J.R. Partington, Analysis of fractional delay systems of retarded and neutral type, Automatica, vol. 38, no. 7, pp. 1133-1138, 2002. [10] C. Bonnet and J.R. Partington, Stabilization of some fractional delay systems of neutral type, Automatica, vol. 43, no. 12, pp. 2047-2053, 2007. [11] Z. V. Rekasius, A stability test for systems with delays, in Proc. Joint Autom. Control Conf., San Francisco, CA, 1980. [12] G. E. Collins, The calculation of multivariate polynomial resultants, J. Assoc. Comput. Mach., vol. 18, no. 4, pp. 515532, 1971. [13] I. M. Gelfand, M. M. Kapranov, and A. V. Zelevinsky, Discriminants, Resultants, and Multidimensional Determinants. Boston, MA: Birkhuser, 1994. [14] D. Wang, Elimination Methods. New York: Springer-Verlag, 2001. [15] V.B. Kolmanovskii, V.R. Nosov, Stability of Functional Differential Equations, Academic Press, London, 1989. [16] C. Hwang and Y.C. Cheng, A numerical algorithm for stability testing of fractional delay systems, Automatica, vol. 42, no. 5, pp. 825-831, 2006. [17] M. A. Pakzad, S. Pakzad, Stability map of fractional order time-delay systems, WSEAS Trans. on Systems, vol.10, no.11, pp. 541-550, 2012. [18] M. A. Pakzad, M.A. Nekoui and S. Pakzad, Stability analysis of time-delayed linear fractional-order systems, International Journal of Control, Automation and Systems, vol.11, no.3, 2013, doi: 10.1007/s12555-012-0164-4.
