Stability Map of Fractional Delay Systems in the ...
2013 American Control Conference (ACC) Washington, DC, USA, June 17-19, 2013
Stability Map of Fractional Delay Systems in the Parametric Space of Delays and Coefficient Sara Pakzad1 , Mohammad Ali Pakzad2 , and Mohammad Ali Nekoui3 has been investigated widely for the standard case (integer order systems).
Abstract— In this paper, we examine an open problem: The stability analysis of fractional order systems with uncertain parameters in both time-delay space and coefficient space. It is evident from the literature that the stability assessment of this class of dynamics remains unsolved. The Rekasius transformation is used as a connection between time-delay space and coefficient space. For generating potential stability switching curves (PSSC) an efficient procedure is proposed to extend the paradigm called, advanced clustering with frequency sweeping (ACFS), to the determination of stability regions in coefficientdelay space. We show that this methodology analytically reveals all possible stability regions exclusively in the space of the delay and coefficient. Two examples are presented to highlight the proposed approach.
For systems with coefficient uncertainty, the Ddecomposition method is used in studying the total set of parameters (controller coefficients or plant characteristics) which provide stability of the system. The technique is often exploited for the design of low-order controllers, see [6]-[8]. Parallel to the stability analysis of time delay system with uncertain coefficients, many papers devoted to the stability analysis of time-delay system with uncertain time delays. There has been a large effort to deal with this problem, see [9]-[15], and others.
I. INTRODUCTION
Though there are a lot of works focusing on the stability analysis of time-delayed system with uncertain coefficients or delays, there is comparatively little work on the stability analysis of time-delayed system with both uncertain coefficients and time delays. A recent paper [15] extended the CTCR method to time delay systems with both coefficient uncertainty and time-delay uncertainty for the standard case.
Time delay is an inherent part of many dynamical and physical systems such as communication systems and biological and chemical processes; and this time delay could lead to inadequate performance and even the instability of the mentioned systems. Therefore, time-delayed systems play significant roles in theoretical as well as practical fields; and this influence can be observed in numerous research articles written on various problems that involve this class of systems [1], [2], [3], [4], [5]. In a generic sense, the fractional-order systems are recognized by non-integer powers of the Laplace variable s. When time delays and fractional-order derivative are appeared in dynamical systems, we have time-delayed fractional-order 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 time delays, against delay and coefficient uncertainties. Usually, the stability analysis of time delayed system with uncertain coefficients or uncertain time delays are treated separately. In the last two decades, the determination of stability region of uncertain parameters in time delay space or coefficient space is a main stream research problem and
¨ urk and Uraz [16] may be For fractional order systems, Ozt¨ the pioneer to consider stability of the fractional delay system with delay and coefficient uncertainty. They have developed the Ruth-Hurwitz criteria for analyzing the stability of some √ special delay systems to those involve fractional power s. From the analitical analysis point of view, the effective analitical approaches to check BIBO stability of fractionalorder delay systems have been discussed in [23] and [25]. The stability robustness against both coefficient uncertainty and delay uncertainty of time-delayed fractional order system is the main question in this article. The main contribution of this document is that we demonstrate for the first time that the stability maps of a fractional order system with uncertain parameters in both time-delay space and coefficient space can be obtained efficiently. The main objective in the stability analysis is to construct all the potential stability switching curves (PSSC), which partition the coefficientdelay space into stable and unstable regions. Obviously, the accuracy and completeness of the analysis strongly depend on finding all the PSSC without any approximations. The method starts with the determination of all possible values of uncertain parameters which result in purely imaginary characteristic roots. The Rekasius transformation is used as a connection between time-delay space and coefficient space. The original idea in this strategy is derived from [9] to achieve the stability map of integer order systems with multiple time delays.
1 S. Pakzad is with the Departments of Electrical Engineering, South Tehran Branch, Islamic Azad University, Tehran, Iran
[email protected] 2 M.A. Pakzad is with Department of Electrical Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
[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
176
II. PRELIMINARIES AND DEFINITIONS
v
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: n
CE(s, K, τ ) =
√
∑ Pℓ ( α s, K) e−ℓτ s
v
Fig. 1.
III. METHODOLOGY
(1)
ℓ=0
τ is non-negative, such that τ ∈ R+ , K ∈ where parameter √ α for ℓ ∈ NN is a real polynomial in the R and Pℓ ( s, K) √ complex variable α s (where α ∈ N). Note that the zeros of characteristic equation (1) are in fact the poles of the system under investigation. We find out from [17] 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 doesnt have any pole at ℜ(s) ≥ 0 (in particular, no poles of fractional order at s=0). For √ fractional order systems, if a 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
CEv (v, K, τ ) =
∑ Pℓ (v, K) e
−ℓτ vα
(2)
ℓ=0
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: π (3) |∠v| ≤ 2α with v ∈ C , which the stable region has been displayed by hatched lines in Fig. 1. Note that under this transformation, the imaginary axis in the s-domain is mapped into the lines π ∠v = ± (4) 2α
The v-stability region for fractional systems
In this section, we propose a method that can extract the 2-D PSSC for fractional order systems with both coefficient uncertainty and delay uncertainty. In other words, we perform the stability analysis of (1) in K-τ parameter space. 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 {τ } the same imaginary root will also exist at all the countably infinite grid points of {τ } = (τ1 +
2π l 2π ), l = 0, 1, 2, ... , τ1 − ≤0 ωc ωc
(6)
This signifies that τ1 is the smallest positive τℓ , τ1 = min(τℓ ), l = 0, 1, 2, ..., (τℓ > 0). Notice that for a fixed ωc the distinct points of (6) generate a grid in {τ } ∈ R+ space with equal grid size, 2π /ωc . Definition 1 (“kernel curves”): Assume that the set of (K, τ1 )|ωc is determined exhaustively in R2 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 (K, τ ). Definition 2 (“offspring curves”): The trajectories of (K, τl ) grid points in (6) for l = 1, 2, ... corresponding to the kernel are called the “offspring curves” or “offspring” in short. They are represented by ℘l (K, τ ) where l identify the l th generation offspring of the kernel τ1 and K, respectively, according to equation (6). Let’s denote the complete set of kernel and offspring by
Let us assume that ; s = ± jω or in other words, s = ω e± jπ /2 are the roots of characteristic equation (1) for a (τ1 , τ2 ) ∈ R2+ . Then for the auxiliary variable, the roots are defined as follows: √ √ v = α s = α ω e j±π /2α (5) √ 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.
n
℘(K, τ ) = ℘0 (K, τ ) + ∑ ℘l (K, τ ) l=1
(7)
The kernel and the offspring constitute the complete (and exhaustive) distribution of (K, τ ) points for which the characteristic equation CE(s, K, τ ) has root sets with at least one imaginary pair. Outside the set of curves ℘(K, τ ) 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) 177
Ω[s|C(s, K, τ ) = 0, τ ∈ R+ , K ∈ R] ∩C0
≡ Ω[s|C(s, K, τ ) = 0, (τ , K) ∈ ℘(K, τ )] ∩C0
imaginary parts of the resulting equation should be separated as follows:
(8)
h(v, K, T )|v= α√ω
is generated only by a small number of contours in (K, τ ). These are the only locations in the (K, τ ) space where the system (1) could transit from stable to unstable posture (or vice versa). These contours ℘(K, τ ) must be determined exhaustively. The determination of the complete set of kernel and offspring is, mathematically speaking, a very challenging problem. In order to achieve this we replace the exponential type transcendental terms in (1) with Rekasius substitution [18]. 1−Ts , T ∈R (9) e−τ s = 1+Ts It is important to note that, this substitution is an exact expression of e−τ s for purely imaginary roots s = ± jω . Moreover, transformation (9) is different from the first-order Pade’ approximation of e−τ s which is
m
and
Not that all ai and bi are real polynomials in T . hℜ and hℑ , which have positive degrees in terms of K, are assumed to have no common factors. Such common factors, if they exist, can be separately studied. Definition 3: The resultant RK with respect to ω and T is the resultant of hℜ and hℑ by eliminating K [19]. The resultant of RK and ∂ RK /∂ T with respect to ω is called the discriminant of RK by eliminating T [20],[21]. Theorem 2: Minimum and maximum positive real roots of the discriminant of ressultant of hℜ and hℑ with respect to ω , that correspond to (K, T ) ∈ 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 [10]. To find the global maximum ω and the global minimum ω , maximaize and minimize ω with respect to τ by studying ∂ ω /∂ τ = 0, which is identical to studying
(11)
This equation describes an asymmetric mapping in which one T is mapped into countably infinite τl which are distributed with a periodicity of 2π /ωc . By inserting equation (10) into equation (2), we have: 1 − T v α ℓ n (12) CEv (v, K, T ) = ∑ Pℓ (v, K) 1 + T vα ℓ=0
∂ω ∂ω ∂T = =0 ∂τ ∂T ∂τ
(17)
Where ∂ T /∂ τ = 0.5(1 + ω 2 T 2 ) 6= 0 as per (11). Since ∂ T /∂ τ 6= 0, we can study ∂ ω /∂ τ = 0 alternatively on ∂ ω /∂ T . At this point, the differentiability of ω with respect to T is essential as established above, and holds for the regular points of RK (hℜ , hℑ ). Under this condition, we can write ∂ω ∂ RK (hℜ , hℑ )/∂ T =0 (18) =− ∂T ∂ RK (hℜ , hℑ )/∂ ω
By multiplying equation (12) by (1 + T vα )n the polynomial form of the characteristic equation is reached: (1 + T vα )nCEv (v, K, T ) = h(v, K, T ) n
(16)
i=0
1 − T vα , T ∈ R, (10) 1 + T vα By examining the amplitude and phase of equation (9), the relationship between T and τ can be obtained as follows:
∑ Pℓ (v, K) (1 − T vα )ℓ (1 + T vα )n−ℓ
m
hℑ = ∑ bi (ω , T )K i = 0
α
=
(15)
i=0
e−τ v =
l = 0, 1, ...
= hℜ (ω , K, T ) + jhℑ (ω , K, T ) (14)
hℜ = ∑ ai (ω , T )K i = 0
1 − 0.5τ s 1 + 0.5τ s Transformation (9) can also be written as:
2 [tan−1 (ωc T ) + l π ], ωc
jπ /2α
In the above relations, hℜ = ℜ(h), hℑ = ℑ(h). For ω to be a zero of (14), hℜ and hℑ should be concurrently zero for some (K, T ). Let us investigate those (K, T ) solutions from
e−τ s ≈
τl =
ce
(13)
ℓ=0
Which leads to ∂ RK (hℜ , hℑ )/∂ T = 0, assuming RK (hℜ , hℑ )/∂ ω 6= 0, Thus, one has two equations, RK and ∂ RK /∂ T , and they should be simultaneously zero. This requires to study the zeros of the resultant of these two equations, particularly by eliminating T . The resultant of RK and ∂ RK /∂ T becomes only a function of ω .
This expression is a polynomial in v of which the coefficients are parametric functions of T and K. 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 √ 1 C(s, K, τ ) then v = s /α = α ωc e± jπ /2α would be a pair of roots for h(v, K, T ). since h(v, K, T ) is a finite degree polynomial with maximum degree of max(deg(pℓ (v)) + α n) then the number of crossing points of (1) that are the real roots of h(v, K, T ) is finite. To find the crossing frequencies set in equation (1), relation (5) should be inserted into (13) and then the real and
Z(ω ) = RT (RK , ∂ RK /∂ T )
(19)
Which is the discriminant by Definition 3. The minimum and maximum positive real zeros of Z(ω ) that correspond to (K, T ) ∈ R2 solotions in (14) are the exact lower and upper bounds of the crossing frequency set, respectively. In the sequel, our new methodology based on advanced clustering with frequency sweeping presented step by step. 178
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 [9]:
Based on (22) and definition given in (21), the root tendency of each time delay τl is written as follows: h π π dv i RT = sgn ℜ sin( ) + jcos( ) 2α 2α d τ v=vc τ =τl
α
[sin( 2πα ) + jcos( 2πα )] ∑nℓ=0 Pℓ (v)ℓvα e−ℓτ v v=vc = sgn ℜ α α P (v) τ =τl ∑nℓ=0 ℓdv e−ℓτ v − Pℓ (v)ℓα vα −1 e−ℓτ v
1) Solve the polynomial equation RK (hℜ , hℑ ) for T ∈ R values. 2) For each T ∈ R found from above, if K ∈ 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. 3) Via (11), calculate the delay values τ corresponding to (K, T ) ∈ R2 pairs, and restart from step 1 increasing ω by an amount of the step size.
α
P (v) ∑nℓ=0 ℓdv e−ℓτ v = sgn ℜ α π π [sin( 2α ) + jcos( 2α )] ∑nℓ=0 Pℓ (v) ℓ vα e−ℓτ v −1 ατ v=vc − π π τ =τ l (sin( 2α ) + jcos( 2α ))v α
P (v) ∑nℓ=0 ℓdv e−ℓτ v v=vc = sgn ℜ α π π n −ℓ α τ v τ =τl [sin( 2α ) + jcos( 2α )] ∑ℓ=0 Pℓ (v)ℓv e (23) The root tendency in each time delay τl is independent from the time delay itself and constant for each crossing frequency, α because e−τ v and Pℓ (v) do not depend on τl . This theorem helps identifying certain sections of ℘0 and the ‘offspring curves’ to be marked as stabilizing transitions along the τ1 (or K) axis or vice versa. Once we detect completely ℘0 (K, τ ), ℘l (K, τ ) ‘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 [22]. 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 K − τ hyperplane where stability switches (i.e., from stable to unstable or vice versa) occur.
B. Direction of Crossing With fixed K to K0 , the root sensitivities associated with each purely imaginary characteristic root crossing jω with respect to , τ , is defined as ∂ C/∂ τ ds s Sτ |s= jωc = = − ∂ C (20) d τ s= jωc /∂ s s= jωc And the corresponding root tendency with respect to delay is given as: Root Tendency = RT |τs= jωc = sgn ℜ Sτs |s= jωc ds dv × = sgn ℜ s= jωc dv d τ √ ±j π v=vc = α ω c e 2α − jπ dv dv (21) = sign ℜ je 2α = sgn ℜ vα −1 dτ dτ v=vc π π dv = sgn ℜ [sin( ) + jcos( )] 2α 2α d τ v=vc π π dv dv − cos( )ℑ = sgn sin( )ℜ v=vc 2α dτ 2α dτ τ =τ
IV. ILLUSTRATIVE EXAMPLE We present three two cases, which display all the features discussed in the text. Example 1: Consider the fractional order system with two time delays in [24], which has the following transfer function
If it is positive, then it is a destabilizing crossing, whereas if it is negative, this means a stabilizing crossing. 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. In other words, the imaginary root always crosses either to C+ (for RT = +1) or to C− (for RT = −1), when the K is kept fixed. If RT = +1, number of unstable roors (NU) increases by 2. If RT = −1, NU decreases by 2. Theorem 3: The root tendency at a crossing, jωc is invariant with respect to time delay τl . 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 τ ∈ R+ . Then one can find dv/d τ for simple roots of (2) as follows:
C1 (s, τ ) =
1/2
s + s +s
1/3
(24)
e−τ s + e−2τ s
A schematic of the closed loop system with a controller is given in Fig.2. The closed-loop transfer function of the system (24) with feedback controller is: G1 (s, K, τ ) =
1 s + K + s + s1/3 e−τ s + e−2τ s 5/6
(25)
1/2
It has been demonstrated in [24] that this system is stable for K = 0 and τ = 0.5. Also, in [25] for K = 0, it has been shown that the system is BIBO stable for τ ∈ (0, 2.3562). Now we evaluate the stability of this system for all the values of K, τ . The characteristic equation of (25) is: 5 1 1 (26) CE1 (s, K, τ ) = s /6 + K + s /2 + s /3 e−τ s + e−2τ s √ By using of auxiliary variable v = 6 s in (26), we will have:
α
dv ∑nℓ=0 Pℓ (v) ℓ vα e−ℓτ v = α α P (v) dτ ∑nℓ=0 ℓdv e−ℓτ v − ∑nℓ=0 Pℓ (v)ℓατ vα −1 e−ℓτ v
1 5/6
(22)
6
6
CEv1 (v, τ1 , τ2 ) = v5 + K + (v3 + v2 )e−τ v + e−2τ v 179
(27)
By eliminating exponential term in (27) , we have the following equation: h(v, K, T ) = (v5 + K)(1 + T v6 )2
(28) + (v3 + v2 )(1 − T v6 )(1 + T v6 ) + (1 − T v6 )2 √ jπ √ And , since v = 6 s = 6 ω e 12 we will have: √ √ √ √ (29) v = a(( 6 + 2) + j( 6 − 2)) √ Where a = 6 ω /4. By inserting equation (29) into equation (28) and making the real and imaginary parts of h(v, K, T ) to zero, and Following Theorem 2, corresponding frequency ranges are found as: [ω , ω ] = [0, 3.1845]
(30)
Upon sweeping ω in this range, the potential stability switching curves of the system is extracted in Fig.3. Set s = 0 in (26) yields K +1 = 0
Fig. 3.
(31) system is unstable for τ = 0.99 and stable for τ = 1 . Also it is showen in [23] that the system CE2 (s, K, τ ) 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 described in this article. for CE2 (s, K, τ ) based on the method √ Using auxiliary variable v = s into (33) the characteristic equation is obtained as:
Solving (31), K=-1. Thus, when K = −1, (27) has one standing root at origin. Root tendency (RT) of the crossing across the origin is investigated for a complete stability picture. For K < −1, one roots is real and lies in the righthalf plane on the positive real axis. In Fig. 3, 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. From Fig. 3, it is clear that when K > 1.57, the example system in (26) is delay independent stable.
r(t)
e(t)
C(s)
2
CEv2 (v, K, τ ) = v3 − 1.5v2 + 4v + 8 − 1.5Kv2 e−τ v
(34)
Replace the exponential term in (34) with Rekasius transformation and refer to (13)
y(t)
h(v, K, T ) = (v3 − 1.5v2 + 4v + 8)(1 + T v2 )
+_
− 1.5Kv2 (1 − T v2 )
(35)
√ π By inserting expression v = ω e j 4 = a(1 + j) in the above equation and equating the real and imaginary parts of the obtained relation to zero, we get:
K Fig. 2.
Stability boundaries and stability regions (shaded) in K − τ plane
The block diagram of the closed loop control system.
hℜ = (6a4 T )K + (4a5 − 6a4 + 8a3 )T + (2a3 − 4a − 8) (36)
Example 2: Consider a fractional delay system that the dynamic equation in transfer-function form is the following: √ −1.5( s)2 e−τ s √ √ √ C2 (s, τ ) = (31) ( s)3 − 1.5( s)2 + 4( s) + 8
and
The closed-loop transfer function of the above system with feedback controller in Fig.2 is: √ −1.5( s)2 e−τ s √ 2 √ √ G2 (s, K, τ ) = √ 3 ( s) − 1.5( s) + 4( s) + 8 − 1.5K( s)2 e−τ s (32) The closed-loop transcendental characteristic equation of this system is √ √ √ √ CE2 (s, K, τ ) = ( s)3 −1.5( s)2 +4( s)+8−1.5K( s)2 e−τ s (33) For K = 1, a very involved calculation scheme based on Cauchy’s integral has been used in [24] to show that this 180
hℑ = (3a2 )K + (4a5 − 8a3 − 16a2 )T + (−2a3 + 3a2 − 4a) (37) Using homomorphism resultant algorithm [19] eliminate K from hℜ and hℑ , RK (hℜ , hℑ ) = (−24a9 + 48a7 + 96a6 )T 2 + (24a7 − 36a6 + 48a5 )T + 6a5 − 12a3 − 24a2 = 0
(38)
Discriminant of the resultant of hℜ and hℑ in Theorem 2 is the resultant of RK and ∂ RK /∂ T with eliminating T ,
∂ RK )=0 (39) ∂T One can now use the efficient method in previous section, in order to extract the stability maps on K − τ domain by RT (RK ,
R EFERENCES
Fig. 4.
[1] M. A. Pakzad, Kalman filter design for time delay systems, WSEAS Trans. on Systems, vol. 11, no. 10, pp. 551-560, 2012. [2] 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. [3] Y. S. Suh, Y. S. Ro, H. J. Kang, and H. H. Lee, Necessary and sufficient stability condition of discrete state delay systems, International Journal of Control, Automation, and Systems, vol. 2, no. 4, pp. 501-508, 2004. [4] M. A. Pakzad, and B. Moaveni, Delay-Dependent state estimation for time delay systems, WSEAS Trans. on Systems and Control, vol. 8 , no. 1, pp. 1-10, 2013. ¨ urk and A. Uraz, An analysis stability test for a certain class [5] N. Ozt¨ of distributed parameter systems with delays, IEEE Trans. on Circuits and Systems, vol. 34, no. 4, pp. 393-396, 1985. [6] J. Ackermann, and D. Kaesbauer, Stable polyhedra in parameter space. Automatica vol.39,pp. 937943, 2003. [7] S.P. Bhattacharyya, R.N. Tantaris, and L.H. Keel, Stabilization of discrete-time systems by firstorder controllers. IEEE Trans Autom Control vol.48,pp. 858861, 2003. [8] C. Hwang, and J.H. Hwang, Stabilization of first-order plus dead-time unstable processes using PID controllers. IEE Proc., Control Theory Appl vol.15,pp. 8994, 2004. [9] 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. [10] 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 [11] 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. 231253, 2005. [12] R. Sipahi and N. Olgac, Complete stability robustness of third-order LTI multiple time-delay systems, Automatica, vol. 41, no. 8, pp. 14131422, 2005. [13] 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. 799810, May 2007. [14] 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. [15] D. Zheng, Z. Ren, J. Fang, Stability analysis of time delayed system with coefficient uncertainty and time delay uncertainty ,European Journal of Control, vol.16, no.1, pp. 5-13, 2010. ¨ urk and A. Uraz, An analytic stability test for a certain class [16] N. Ozt¨ of distributed parameter systems with a distributed lag, IEEE Trans. Autom. Control, vol. 29, no. 4, pp. 368 - 370, 1985. [17] 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. [18] Z. V. Rekasius, A stability test for systems with delays, in Proc. Joint Autom. Control Conf., San Francisco, CA, 1980. [19] G. E. Collins, The calculation of multivariate polynomial resultants, J. Assoc. Comput. Mach., vol. 18, no. 4, pp. 515532, 1971. [20] I. M. Gelfand, M. M. Kapranov, and A. V. Zelevinsky, Discriminants, Resultants, and Multidimensional Determinants. Boston, MA: Birkhuser, 1994. [21] D. Wang, Elimination Methods. New York: Springer-Verlag, 2001. [22] V.B. Kolmanovskii, V.R. Nosov, Stability of Functional Differential Equations, Academic Press, London, 1989. [23] 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. [24] 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. [25] 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 boundaries and stability regions (shaded) in K − τ plane
sweeping the frequency from 0 to 8. The potential stability switching curves of the system is extracted in Fig. 4, 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.4 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. From Fig. 4, it is evident that for −0.99 < K < 0.99, the example system is delay independent stable. V. CONCLUSIONS In this article, we have presented a new procedure based on, Advanced Clustering with Frequency Sweeping (ACFS), method to analyze the stability of linear fractional order systems with both coefficient uncertainty and time-delay uncertainty. 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. The proposed methodology detects all the stability regions precisely, in the space of coefficient and the time delays. This map, in fact, is the exact display of robustness against delay and coefficient uncertainties. After finding stability boundaries in parametric space, the number of unstable roots in each region is calculated with the definition of root tendency on the boundaries of each region. two illustrative examples presented to confirm the proposed method results. ACKNOWLEDGMENT This work was supported by the Science and Research Branch, Islamic Azad University, Tehran, Iran. 181
Stability Map of Fractional Delay Systems in the Parametric Space of Delays and Coefficient Sara Pakzad1 , Mohammad Ali Pakzad2 , and Mohammad Ali Nekoui3 has been investigated widely for the standard case (integer order systems).
Abstract— In this paper, we examine an open problem: The stability analysis of fractional order systems with uncertain parameters in both time-delay space and coefficient space. It is evident from the literature that the stability assessment of this class of dynamics remains unsolved. The Rekasius transformation is used as a connection between time-delay space and coefficient space. For generating potential stability switching curves (PSSC) an efficient procedure is proposed to extend the paradigm called, advanced clustering with frequency sweeping (ACFS), to the determination of stability regions in coefficientdelay space. We show that this methodology analytically reveals all possible stability regions exclusively in the space of the delay and coefficient. Two examples are presented to highlight the proposed approach.
For systems with coefficient uncertainty, the Ddecomposition method is used in studying the total set of parameters (controller coefficients or plant characteristics) which provide stability of the system. The technique is often exploited for the design of low-order controllers, see [6]-[8]. Parallel to the stability analysis of time delay system with uncertain coefficients, many papers devoted to the stability analysis of time-delay system with uncertain time delays. There has been a large effort to deal with this problem, see [9]-[15], and others.
I. INTRODUCTION
Though there are a lot of works focusing on the stability analysis of time-delayed system with uncertain coefficients or delays, there is comparatively little work on the stability analysis of time-delayed system with both uncertain coefficients and time delays. A recent paper [15] extended the CTCR method to time delay systems with both coefficient uncertainty and time-delay uncertainty for the standard case.
Time delay is an inherent part of many dynamical and physical systems such as communication systems and biological and chemical processes; and this time delay could lead to inadequate performance and even the instability of the mentioned systems. Therefore, time-delayed systems play significant roles in theoretical as well as practical fields; and this influence can be observed in numerous research articles written on various problems that involve this class of systems [1], [2], [3], [4], [5]. In a generic sense, the fractional-order systems are recognized by non-integer powers of the Laplace variable s. When time delays and fractional-order derivative are appeared in dynamical systems, we have time-delayed fractional-order 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 time delays, against delay and coefficient uncertainties. Usually, the stability analysis of time delayed system with uncertain coefficients or uncertain time delays are treated separately. In the last two decades, the determination of stability region of uncertain parameters in time delay space or coefficient space is a main stream research problem and
¨ urk and Uraz [16] may be For fractional order systems, Ozt¨ the pioneer to consider stability of the fractional delay system with delay and coefficient uncertainty. They have developed the Ruth-Hurwitz criteria for analyzing the stability of some √ special delay systems to those involve fractional power s. From the analitical analysis point of view, the effective analitical approaches to check BIBO stability of fractionalorder delay systems have been discussed in [23] and [25]. The stability robustness against both coefficient uncertainty and delay uncertainty of time-delayed fractional order system is the main question in this article. The main contribution of this document is that we demonstrate for the first time that the stability maps of a fractional order system with uncertain parameters in both time-delay space and coefficient space can be obtained efficiently. The main objective in the stability analysis is to construct all the potential stability switching curves (PSSC), which partition the coefficientdelay space into stable and unstable regions. Obviously, the accuracy and completeness of the analysis strongly depend on finding all the PSSC without any approximations. The method starts with the determination of all possible values of uncertain parameters which result in purely imaginary characteristic roots. The Rekasius transformation is used as a connection between time-delay space and coefficient space. The original idea in this strategy is derived from [9] to achieve the stability map of integer order systems with multiple time delays.
1 S. Pakzad is with the Departments of Electrical Engineering, South Tehran Branch, Islamic Azad University, Tehran, Iran
[email protected] 2 M.A. Pakzad is with Department of Electrical Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
[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
176
II. PRELIMINARIES AND DEFINITIONS
v
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: n
CE(s, K, τ ) =
√
∑ Pℓ ( α s, K) e−ℓτ s
v
Fig. 1.
III. METHODOLOGY
(1)
ℓ=0
τ is non-negative, such that τ ∈ R+ , K ∈ where parameter √ α for ℓ ∈ NN is a real polynomial in the R and Pℓ ( s, K) √ complex variable α s (where α ∈ N). Note that the zeros of characteristic equation (1) are in fact the poles of the system under investigation. We find out from [17] 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 doesnt have any pole at ℜ(s) ≥ 0 (in particular, no poles of fractional order at s=0). For √ fractional order systems, if a 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
CEv (v, K, τ ) =
∑ Pℓ (v, K) e
−ℓτ vα
(2)
ℓ=0
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: π (3) |∠v| ≤ 2α with v ∈ C , which the stable region has been displayed by hatched lines in Fig. 1. Note that under this transformation, the imaginary axis in the s-domain is mapped into the lines π ∠v = ± (4) 2α
The v-stability region for fractional systems
In this section, we propose a method that can extract the 2-D PSSC for fractional order systems with both coefficient uncertainty and delay uncertainty. In other words, we perform the stability analysis of (1) in K-τ parameter space. 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 {τ } the same imaginary root will also exist at all the countably infinite grid points of {τ } = (τ1 +
2π l 2π ), l = 0, 1, 2, ... , τ1 − ≤0 ωc ωc
(6)
This signifies that τ1 is the smallest positive τℓ , τ1 = min(τℓ ), l = 0, 1, 2, ..., (τℓ > 0). Notice that for a fixed ωc the distinct points of (6) generate a grid in {τ } ∈ R+ space with equal grid size, 2π /ωc . Definition 1 (“kernel curves”): Assume that the set of (K, τ1 )|ωc is determined exhaustively in R2 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 (K, τ ). Definition 2 (“offspring curves”): The trajectories of (K, τl ) grid points in (6) for l = 1, 2, ... corresponding to the kernel are called the “offspring curves” or “offspring” in short. They are represented by ℘l (K, τ ) where l identify the l th generation offspring of the kernel τ1 and K, respectively, according to equation (6). Let’s denote the complete set of kernel and offspring by
Let us assume that ; s = ± jω or in other words, s = ω e± jπ /2 are the roots of characteristic equation (1) for a (τ1 , τ2 ) ∈ R2+ . Then for the auxiliary variable, the roots are defined as follows: √ √ v = α s = α ω e j±π /2α (5) √ 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.
n
℘(K, τ ) = ℘0 (K, τ ) + ∑ ℘l (K, τ ) l=1
(7)
The kernel and the offspring constitute the complete (and exhaustive) distribution of (K, τ ) points for which the characteristic equation CE(s, K, τ ) has root sets with at least one imaginary pair. Outside the set of curves ℘(K, τ ) 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) 177
Ω[s|C(s, K, τ ) = 0, τ ∈ R+ , K ∈ R] ∩C0
≡ Ω[s|C(s, K, τ ) = 0, (τ , K) ∈ ℘(K, τ )] ∩C0
imaginary parts of the resulting equation should be separated as follows:
(8)
h(v, K, T )|v= α√ω
is generated only by a small number of contours in (K, τ ). These are the only locations in the (K, τ ) space where the system (1) could transit from stable to unstable posture (or vice versa). These contours ℘(K, τ ) must be determined exhaustively. The determination of the complete set of kernel and offspring is, mathematically speaking, a very challenging problem. In order to achieve this we replace the exponential type transcendental terms in (1) with Rekasius substitution [18]. 1−Ts , T ∈R (9) e−τ s = 1+Ts It is important to note that, this substitution is an exact expression of e−τ s for purely imaginary roots s = ± jω . Moreover, transformation (9) is different from the first-order Pade’ approximation of e−τ s which is
m
and
Not that all ai and bi are real polynomials in T . hℜ and hℑ , which have positive degrees in terms of K, are assumed to have no common factors. Such common factors, if they exist, can be separately studied. Definition 3: The resultant RK with respect to ω and T is the resultant of hℜ and hℑ by eliminating K [19]. The resultant of RK and ∂ RK /∂ T with respect to ω is called the discriminant of RK by eliminating T [20],[21]. Theorem 2: Minimum and maximum positive real roots of the discriminant of ressultant of hℜ and hℑ with respect to ω , that correspond to (K, T ) ∈ 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 [10]. To find the global maximum ω and the global minimum ω , maximaize and minimize ω with respect to τ by studying ∂ ω /∂ τ = 0, which is identical to studying
(11)
This equation describes an asymmetric mapping in which one T is mapped into countably infinite τl which are distributed with a periodicity of 2π /ωc . By inserting equation (10) into equation (2), we have: 1 − T v α ℓ n (12) CEv (v, K, T ) = ∑ Pℓ (v, K) 1 + T vα ℓ=0
∂ω ∂ω ∂T = =0 ∂τ ∂T ∂τ
(17)
Where ∂ T /∂ τ = 0.5(1 + ω 2 T 2 ) 6= 0 as per (11). Since ∂ T /∂ τ 6= 0, we can study ∂ ω /∂ τ = 0 alternatively on ∂ ω /∂ T . At this point, the differentiability of ω with respect to T is essential as established above, and holds for the regular points of RK (hℜ , hℑ ). Under this condition, we can write ∂ω ∂ RK (hℜ , hℑ )/∂ T =0 (18) =− ∂T ∂ RK (hℜ , hℑ )/∂ ω
By multiplying equation (12) by (1 + T vα )n the polynomial form of the characteristic equation is reached: (1 + T vα )nCEv (v, K, T ) = h(v, K, T ) n
(16)
i=0
1 − T vα , T ∈ R, (10) 1 + T vα By examining the amplitude and phase of equation (9), the relationship between T and τ can be obtained as follows:
∑ Pℓ (v, K) (1 − T vα )ℓ (1 + T vα )n−ℓ
m
hℑ = ∑ bi (ω , T )K i = 0
α
=
(15)
i=0
e−τ v =
l = 0, 1, ...
= hℜ (ω , K, T ) + jhℑ (ω , K, T ) (14)
hℜ = ∑ ai (ω , T )K i = 0
1 − 0.5τ s 1 + 0.5τ s Transformation (9) can also be written as:
2 [tan−1 (ωc T ) + l π ], ωc
jπ /2α
In the above relations, hℜ = ℜ(h), hℑ = ℑ(h). For ω to be a zero of (14), hℜ and hℑ should be concurrently zero for some (K, T ). Let us investigate those (K, T ) solutions from
e−τ s ≈
τl =
ce
(13)
ℓ=0
Which leads to ∂ RK (hℜ , hℑ )/∂ T = 0, assuming RK (hℜ , hℑ )/∂ ω 6= 0, Thus, one has two equations, RK and ∂ RK /∂ T , and they should be simultaneously zero. This requires to study the zeros of the resultant of these two equations, particularly by eliminating T . The resultant of RK and ∂ RK /∂ T becomes only a function of ω .
This expression is a polynomial in v of which the coefficients are parametric functions of T and K. 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 √ 1 C(s, K, τ ) then v = s /α = α ωc e± jπ /2α would be a pair of roots for h(v, K, T ). since h(v, K, T ) is a finite degree polynomial with maximum degree of max(deg(pℓ (v)) + α n) then the number of crossing points of (1) that are the real roots of h(v, K, T ) is finite. To find the crossing frequencies set in equation (1), relation (5) should be inserted into (13) and then the real and
Z(ω ) = RT (RK , ∂ RK /∂ T )
(19)
Which is the discriminant by Definition 3. The minimum and maximum positive real zeros of Z(ω ) that correspond to (K, T ) ∈ R2 solotions in (14) are the exact lower and upper bounds of the crossing frequency set, respectively. In the sequel, our new methodology based on advanced clustering with frequency sweeping presented step by step. 178
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 [9]:
Based on (22) and definition given in (21), the root tendency of each time delay τl is written as follows: h π π dv i RT = sgn ℜ sin( ) + jcos( ) 2α 2α d τ v=vc τ =τl
α
[sin( 2πα ) + jcos( 2πα )] ∑nℓ=0 Pℓ (v)ℓvα e−ℓτ v v=vc = sgn ℜ α α P (v) τ =τl ∑nℓ=0 ℓdv e−ℓτ v − Pℓ (v)ℓα vα −1 e−ℓτ v
1) Solve the polynomial equation RK (hℜ , hℑ ) for T ∈ R values. 2) For each T ∈ R found from above, if K ∈ 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. 3) Via (11), calculate the delay values τ corresponding to (K, T ) ∈ R2 pairs, and restart from step 1 increasing ω by an amount of the step size.
α
P (v) ∑nℓ=0 ℓdv e−ℓτ v = sgn ℜ α π π [sin( 2α ) + jcos( 2α )] ∑nℓ=0 Pℓ (v) ℓ vα e−ℓτ v −1 ατ v=vc − π π τ =τ l (sin( 2α ) + jcos( 2α ))v α
P (v) ∑nℓ=0 ℓdv e−ℓτ v v=vc = sgn ℜ α π π n −ℓ α τ v τ =τl [sin( 2α ) + jcos( 2α )] ∑ℓ=0 Pℓ (v)ℓv e (23) The root tendency in each time delay τl is independent from the time delay itself and constant for each crossing frequency, α because e−τ v and Pℓ (v) do not depend on τl . This theorem helps identifying certain sections of ℘0 and the ‘offspring curves’ to be marked as stabilizing transitions along the τ1 (or K) axis or vice versa. Once we detect completely ℘0 (K, τ ), ℘l (K, τ ) ‘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 [22]. 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 K − τ hyperplane where stability switches (i.e., from stable to unstable or vice versa) occur.
B. Direction of Crossing With fixed K to K0 , the root sensitivities associated with each purely imaginary characteristic root crossing jω with respect to , τ , is defined as ∂ C/∂ τ ds s Sτ |s= jωc = = − ∂ C (20) d τ s= jωc /∂ s s= jωc And the corresponding root tendency with respect to delay is given as: Root Tendency = RT |τs= jωc = sgn ℜ Sτs |s= jωc ds dv × = sgn ℜ s= jωc dv d τ √ ±j π v=vc = α ω c e 2α − jπ dv dv (21) = sign ℜ je 2α = sgn ℜ vα −1 dτ dτ v=vc π π dv = sgn ℜ [sin( ) + jcos( )] 2α 2α d τ v=vc π π dv dv − cos( )ℑ = sgn sin( )ℜ v=vc 2α dτ 2α dτ τ =τ
IV. ILLUSTRATIVE EXAMPLE We present three two cases, which display all the features discussed in the text. Example 1: Consider the fractional order system with two time delays in [24], which has the following transfer function
If it is positive, then it is a destabilizing crossing, whereas if it is negative, this means a stabilizing crossing. 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. In other words, the imaginary root always crosses either to C+ (for RT = +1) or to C− (for RT = −1), when the K is kept fixed. If RT = +1, number of unstable roors (NU) increases by 2. If RT = −1, NU decreases by 2. Theorem 3: The root tendency at a crossing, jωc is invariant with respect to time delay τl . 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 τ ∈ R+ . Then one can find dv/d τ for simple roots of (2) as follows:
C1 (s, τ ) =
1/2
s + s +s
1/3
(24)
e−τ s + e−2τ s
A schematic of the closed loop system with a controller is given in Fig.2. The closed-loop transfer function of the system (24) with feedback controller is: G1 (s, K, τ ) =
1 s + K + s + s1/3 e−τ s + e−2τ s 5/6
(25)
1/2
It has been demonstrated in [24] that this system is stable for K = 0 and τ = 0.5. Also, in [25] for K = 0, it has been shown that the system is BIBO stable for τ ∈ (0, 2.3562). Now we evaluate the stability of this system for all the values of K, τ . The characteristic equation of (25) is: 5 1 1 (26) CE1 (s, K, τ ) = s /6 + K + s /2 + s /3 e−τ s + e−2τ s √ By using of auxiliary variable v = 6 s in (26), we will have:
α
dv ∑nℓ=0 Pℓ (v) ℓ vα e−ℓτ v = α α P (v) dτ ∑nℓ=0 ℓdv e−ℓτ v − ∑nℓ=0 Pℓ (v)ℓατ vα −1 e−ℓτ v
1 5/6
(22)
6
6
CEv1 (v, τ1 , τ2 ) = v5 + K + (v3 + v2 )e−τ v + e−2τ v 179
(27)
By eliminating exponential term in (27) , we have the following equation: h(v, K, T ) = (v5 + K)(1 + T v6 )2
(28) + (v3 + v2 )(1 − T v6 )(1 + T v6 ) + (1 − T v6 )2 √ jπ √ And , since v = 6 s = 6 ω e 12 we will have: √ √ √ √ (29) v = a(( 6 + 2) + j( 6 − 2)) √ Where a = 6 ω /4. By inserting equation (29) into equation (28) and making the real and imaginary parts of h(v, K, T ) to zero, and Following Theorem 2, corresponding frequency ranges are found as: [ω , ω ] = [0, 3.1845]
(30)
Upon sweeping ω in this range, the potential stability switching curves of the system is extracted in Fig.3. Set s = 0 in (26) yields K +1 = 0
Fig. 3.
(31) system is unstable for τ = 0.99 and stable for τ = 1 . Also it is showen in [23] that the system CE2 (s, K, τ ) 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 described in this article. for CE2 (s, K, τ ) based on the method √ Using auxiliary variable v = s into (33) the characteristic equation is obtained as:
Solving (31), K=-1. Thus, when K = −1, (27) has one standing root at origin. Root tendency (RT) of the crossing across the origin is investigated for a complete stability picture. For K < −1, one roots is real and lies in the righthalf plane on the positive real axis. In Fig. 3, 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. From Fig. 3, it is clear that when K > 1.57, the example system in (26) is delay independent stable.
r(t)
e(t)
C(s)
2
CEv2 (v, K, τ ) = v3 − 1.5v2 + 4v + 8 − 1.5Kv2 e−τ v
(34)
Replace the exponential term in (34) with Rekasius transformation and refer to (13)
y(t)
h(v, K, T ) = (v3 − 1.5v2 + 4v + 8)(1 + T v2 )
+_
− 1.5Kv2 (1 − T v2 )
(35)
√ π By inserting expression v = ω e j 4 = a(1 + j) in the above equation and equating the real and imaginary parts of the obtained relation to zero, we get:
K Fig. 2.
Stability boundaries and stability regions (shaded) in K − τ plane
The block diagram of the closed loop control system.
hℜ = (6a4 T )K + (4a5 − 6a4 + 8a3 )T + (2a3 − 4a − 8) (36)
Example 2: Consider a fractional delay system that the dynamic equation in transfer-function form is the following: √ −1.5( s)2 e−τ s √ √ √ C2 (s, τ ) = (31) ( s)3 − 1.5( s)2 + 4( s) + 8
and
The closed-loop transfer function of the above system with feedback controller in Fig.2 is: √ −1.5( s)2 e−τ s √ 2 √ √ G2 (s, K, τ ) = √ 3 ( s) − 1.5( s) + 4( s) + 8 − 1.5K( s)2 e−τ s (32) The closed-loop transcendental characteristic equation of this system is √ √ √ √ CE2 (s, K, τ ) = ( s)3 −1.5( s)2 +4( s)+8−1.5K( s)2 e−τ s (33) For K = 1, a very involved calculation scheme based on Cauchy’s integral has been used in [24] to show that this 180
hℑ = (3a2 )K + (4a5 − 8a3 − 16a2 )T + (−2a3 + 3a2 − 4a) (37) Using homomorphism resultant algorithm [19] eliminate K from hℜ and hℑ , RK (hℜ , hℑ ) = (−24a9 + 48a7 + 96a6 )T 2 + (24a7 − 36a6 + 48a5 )T + 6a5 − 12a3 − 24a2 = 0
(38)
Discriminant of the resultant of hℜ and hℑ in Theorem 2 is the resultant of RK and ∂ RK /∂ T with eliminating T ,
∂ RK )=0 (39) ∂T One can now use the efficient method in previous section, in order to extract the stability maps on K − τ domain by RT (RK ,
R EFERENCES
Fig. 4.
[1] M. A. Pakzad, Kalman filter design for time delay systems, WSEAS Trans. on Systems, vol. 11, no. 10, pp. 551-560, 2012. [2] 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. [3] Y. S. Suh, Y. S. Ro, H. J. Kang, and H. H. Lee, Necessary and sufficient stability condition of discrete state delay systems, International Journal of Control, Automation, and Systems, vol. 2, no. 4, pp. 501-508, 2004. [4] M. A. Pakzad, and B. Moaveni, Delay-Dependent state estimation for time delay systems, WSEAS Trans. on Systems and Control, vol. 8 , no. 1, pp. 1-10, 2013. ¨ urk and A. Uraz, An analysis stability test for a certain class [5] N. Ozt¨ of distributed parameter systems with delays, IEEE Trans. on Circuits and Systems, vol. 34, no. 4, pp. 393-396, 1985. [6] J. Ackermann, and D. Kaesbauer, Stable polyhedra in parameter space. Automatica vol.39,pp. 937943, 2003. [7] S.P. Bhattacharyya, R.N. Tantaris, and L.H. Keel, Stabilization of discrete-time systems by firstorder controllers. IEEE Trans Autom Control vol.48,pp. 858861, 2003. [8] C. Hwang, and J.H. Hwang, Stabilization of first-order plus dead-time unstable processes using PID controllers. IEE Proc., Control Theory Appl vol.15,pp. 8994, 2004. [9] 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. [10] 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 [11] 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. 231253, 2005. [12] R. Sipahi and N. Olgac, Complete stability robustness of third-order LTI multiple time-delay systems, Automatica, vol. 41, no. 8, pp. 14131422, 2005. [13] 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. 799810, May 2007. [14] 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. [15] D. Zheng, Z. Ren, J. Fang, Stability analysis of time delayed system with coefficient uncertainty and time delay uncertainty ,European Journal of Control, vol.16, no.1, pp. 5-13, 2010. ¨ urk and A. Uraz, An analytic stability test for a certain class [16] N. Ozt¨ of distributed parameter systems with a distributed lag, IEEE Trans. Autom. Control, vol. 29, no. 4, pp. 368 - 370, 1985. [17] 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. [18] Z. V. Rekasius, A stability test for systems with delays, in Proc. Joint Autom. Control Conf., San Francisco, CA, 1980. [19] G. E. Collins, The calculation of multivariate polynomial resultants, J. Assoc. Comput. Mach., vol. 18, no. 4, pp. 515532, 1971. [20] I. M. Gelfand, M. M. Kapranov, and A. V. Zelevinsky, Discriminants, Resultants, and Multidimensional Determinants. Boston, MA: Birkhuser, 1994. [21] D. Wang, Elimination Methods. New York: Springer-Verlag, 2001. [22] V.B. Kolmanovskii, V.R. Nosov, Stability of Functional Differential Equations, Academic Press, London, 1989. [23] 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. [24] 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. [25] 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 boundaries and stability regions (shaded) in K − τ plane
sweeping the frequency from 0 to 8. The potential stability switching curves of the system is extracted in Fig. 4, 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.4 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. From Fig. 4, it is evident that for −0.99 < K < 0.99, the example system is delay independent stable. V. CONCLUSIONS In this article, we have presented a new procedure based on, Advanced Clustering with Frequency Sweeping (ACFS), method to analyze the stability of linear fractional order systems with both coefficient uncertainty and time-delay uncertainty. 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. The proposed methodology detects all the stability regions precisely, in the space of coefficient and the time delays. This map, in fact, is the exact display of robustness against delay and coefficient uncertainties. After finding stability boundaries in parametric space, the number of unstable roots in each region is calculated with the definition of root tendency on the boundaries of each region. two illustrative examples presented to confirm the proposed method results. ACKNOWLEDGMENT This work was supported by the Science and Research Branch, Islamic Azad University, Tehran, Iran. 181
