Cent. Eur. J. Phys.

DOI: 10.2478/s11534-013-0226-0

VERS ITA

Central European Journal of Physics

General formula for stability testing of linear systems with fractional-delay characteristic equation

Research Article

Farshad Merrikh-Bayat*

Faculty of Electrical Engineering, University of Zanjan Zanjan, Iran, P.O.Box313

Received 03 February 2013; accepted 29 March 2013

Abstract: In some applications (especially in the filed of control theory) the characteristic equation of system contains

fractional powers of the Laplace variable s possibly in combination with exponentials of fractional powers of s. The aim of this paper is to propose an easy-to-use and effective formula for bounded-input bounded-output (BIBO) stability testing of a linear time-invariant system with fractional-delay characteristic equation

in the general form of A(s) = P0(s) + YLNl-i Pi(s)exp ( - Zispi) = 0, where Pt(s) (i = 0.....N) are the

so-called fractional-order polynomials and Zi and ft are positive real constants. The proposed formula determines the number of the roots of such a characteristic equation in the right half-plane of the first Riemann sheet by applying Rouche's theorem. Numerical simulations are also presented to confirm the efficiency of the proposed formula.

PACS C200B): 02.60.-x, 02.70.-c, 89.75.-k

Keywords: fractional-delay system • stability testing • Riemann surface • Rouche's theorem • simulation

© Versitasp. zo.o.

1. Introduction

As It is also summarized In [1], fractional-order differential equations can model many real-world physical systems better than classical integer-order differential equations. For example, it is a well-known fact that fractional-order differential equations are more effective for modeling materials having memory [2], mass diffusion and heat conduction [3], viscous dampers [4], viscoelastic relaxation [5], and so on.

In dealing with fractional-order differential equations the

^E-mail: f.bayat@znu.ac.ir

generalized operator 0Dt, where a and t are the limits and a is the order of operation, is often defined by the Riemann-Liouville definition (a > 0):

™ , > 1 d" ft x(t)

aDax (t) = —-T— --V-7<dr,

a t () r(n - a) dt" Ja (t — r)a—n+1 , (1)

n — 1 < a < ",

where T(.) is the well-known Euler's gamma function. Other definitions for fractional differential operators can also be found in the literature [6]. Fortunately all of these definitions lead to the same result in the Laplace transform domain assuming zero initial conditions for variables. More precisely, it is shown that the Laplace transform (L) of a fractional derivative of a signal x(t) considering null

Springer

initial conditions is given by:

L {Dax(t)} = saX(s),

where X(s) = L{x(t)}.

In engineering it is a common practice to describe linear time-invariant (LTI) systems by their transfer function; which is equal to the Laplace transform of the impulse response of system. Many physical phenomena including certain types of dielectric polarization, the metal-liquid interface polarization, and the distributed constant system in a transmission line can usually be represented by a transfer function in fractional powers of s [7]-[9]. In some cases the transfer function of system may even consist of fractional powers of s in combination with exponentials of fractional powers of s. For example, the transfer functions:

G (s) =

G (s) =

cosh (xovS) VSsinh ( vS),

sinh (xovS) sinh ( vS) '

0 < xo < 1,

0 < xo < 1,

appear in boundary control of one-dimensional heat equation with Neumann and Dirichlet boundary conditions [10]. Other examples of this type can be found in [10]-[12]. Note that a branch point of both of the transfer functions given in (3) and (4) is located at the origin. In practice the branch point of the transfer function of a system may be at a point rather than at the origin. As an example of this type, consider the asymmetric transmission line shown in Fig. 1. In this figure R, C, L, and G show the resistance, capacitance, inductance, and conductance in the unit of length, respectively, and the length of line is assumed equal to L. According to the Krishoff's law equations governing this system are as the following:

As it can be observed, the transfer function given in (8) has branch points at -R/L and -G/C. It is understood that analysis and control of LTI systems with nontrivial transfer functions such as those given in (3), (4), and (8) are much more difficult than systems with classical transfer functions. See also [13]-[16] to gain more information about some of the recently developed methods and tools in this area.

In the literature, a transfer function in fractional powers of s in combination with exponentials of fractional powers of s is commonly referred to as a fractional-delay transfer function. The transfer functions given in (3) and (4) are two examples of fractional-delay transfer functions. The fractional-delay characteristic equation is obtained by equating the denominator of a fractional-delay transfer function to zero. More formally, the fractional-delay characteristic equation under consideration in this paper can be expressed in the general form of

A(s) = Po(s) + ^ PL(s) exp (-ZLsp/) = 0, (9)

where Z and pi are positive real constants and Pi(s) (i = 1,..., N) are fractional-order polynomials in the form of

Pl(s) = £ alksak,

where aik and aik are real and positive real constants, respectively, and

Po(s) = sa" + an_i sa-1 + ... + ai sai + ao, (11)

where, without any loss of generality, it is assumed that the powers of s in (11) satisfy the following relations:

di(x,t >= Gv (x,t ) + C-V M )

dV M >= Ri(x,t ) + L^ )

subject to the boundary conditions:

an > an-i > ••• > ai > 0.

Note that although the transfer function given in (8) cannot be written in the form (9) the following discussions are general and still valid for it.

v(0, t) = ^s(t), i(i, t) = 0. (7)

Taking the Laplace transform from both sides of (5), (6) and (7) (assuming that the initial condition of system is equal to zero), and after algebraic calculations the transfer function of system is obtained as follows:

r (s) = VM =_1_ (8)

(s) Vs(s) cosh (i(R + is)1/2(r + Cs)1/2). (8)

i(x,t)

Figure 1. The model of an asymmetric transmission line. The transfer function V0(s)/Vs(s) has branch points at -R/i and

In control engineering, fractional-delay characteristic equations commonly appear in feedback systems where an LTI process with dead time is controlled by a fractional-order controller. For example, consider the unity-feedback system shown in Fig. 2 where a process with transfer function

is controlled with the so-ca lled P\AD" controller with transfer function [17]:

C(s) = Kp (1 + + Tds

where KP, Ti, Td, A, and " are unknown parameters of the controller to be determined. It can be verified that the characteristic equation of this feedback system (i.e., the equation obtained by equating the denominator of closed-loop transfer function from R to Y to zero) is as the following

A(s) = TsA(1 + sT) + KpK ( TsA + 1 + T<TdsA+") e-sL = 0,

which can be considered as a special case of (9) with N = 1,

Po(s)= TisA(1 + sT), (16)

Pi(s) = KpK (TS + 1 + TJdsA+),

Z1 = L, and ft = 1.

Stability analysis of feedback systems with fractional-delay characteristic equation is of crucial importance in control engineering. For example, if in the above-mentioned feedback system one tries to find the optimal values of KP, Ti, Td, A, and " by means of meta-heuristic optimization algorithms such that a certain cost function (e.g., the ISE performance index corresponding to the tracking of unit step command) is minimized, a method to check the feasibility of the solutions generated by meta-heuristic optimization algorithm from the stability point of

C(s) G(s)

Controller

Process

Figure 2. Unity feedback system. Stability analysis of this system when both the delay (either of classical or fractional type) and fractional powers of s exist in the loop is a challenging task.

view is needed. It should be noted that since in such optimization problems the cost function is usually expressed in the frequency domain (by applying Parseval's theorem), the resulted optimal controller may destabilize the feedback system while it minimizes the cost function in frequency domain [18]. In some practical cases one may even face more complicated problems. For example, it may be desired to design a PIAD" controller for a process which is modeled by a fractional-order transfer function possibly in combination with delay term (see, for example, [6] for modeling a heating furnace by fractional-order transfer function which is proved to be more accurate than the integer-order one). Stability analysis of the feedback system becomes more difficult in dealing with such cases.

In general stability analysis of the feedback system when both the time delay (either of integer or fractional type) and fractional-order terms exist in the loop is a challenging task. Even the stability analysis of a feedback system which consists of both the PIAD" controller and a process with dead-time is not straightforward. So far many researchers have tried to develop analytical or numerical methods for stability testing of systems with fractional-delay characteristic equations (see, for example, [19]-[23] and the references therein for more information on this subject). Probably the most famous analytical method for stability testing of fractional-order systems (as a special case of fractional-delay systems) is the sector stability test of Matignon [24], which is initially reported in the work of Ikeda and Takahashi [25]. Application of this method is limited to the case where the sigma term does not exist in (9) and P0(s) is of commensurate order. Few numerical algorithms for stability testing of (9) can also be found in the literature (see, for example, [26] and [19] and the references therein for more information). As far as the author knows, all of these methods suffer from the limitation that they can be applied only to a certain class of fractional-delay systems [26], or the results are of probabilistic nature [19].

The aim of this paper is to propose a formula for determining the number of unstable roots of (9). The proposed formula is actually a generalization of the method already proposed by author in [26]. However, the formula developed in this paper has the advantage of being much simpler compared to the one presented in [26], and moreover, it can be easily applied to a more general form of fractional-delay characteristic equations.

The rest of this paper is organized as follows. The proposed formula for stability testing of fractional-delay characteristic equations is presented in Section 2. Four numerical examples are studied in Section 3, and finally, Section 4 concludes the paper.

2. Proposed formula for stability testing of fractional-delay characteristic equations

The first step In dealing with multi-valued complex functions (such as the one presented in (9)) is to construct the domain of definition of the function appropriately. The domain of definition of the characteristic function given in (9) is, in general, in the form of a Riemann surface with infinite number of Riemann sheets where the origin is a branch point and the branch cut is considered (arbitrarily) at R-. Equation A(s) = 0 as defined in (9) has, in general, infinite number of roots which are distributed on this Riemann surface. It is understood that a system with characteristic equation (9) is stable if and only if it does not have any roots in the right half-plane of the first Riemann sheet [19, 27]. Hence, stability analysis of a system with characteristic equation (9) is equivalent to investigation for the roots of A(s) = 0 in the right half-plane of the first Riemann sheet. The following uses Rouche's theorem for this purpose.

First, let us briefly review Rouche's theorem. Consider the complex function f : C ^ C which has zeros of orders mj,...,mk respectively at and does not have

any poles. This function can be written as

f (s) = g(s)(s - Zi)mi x (s - z2f2 x-.-x (s - z* p, (18)

where g(s) has neither pole nor zero. Taking the natural logarithm from both sides of (18) leads to

In f (s) = In g(s) + mj ln(s — z-j) + m2 ln(s — z2) + ... ... + mk ln(s - z*).

Derivation with respect to s from both sides of (19) yields

f£? = + ^ + ^ + ... + A. (20)

f (s) g(s) s - z- s - z2 s - zk

Now let y be a simple, closed, counterclockwise contour such that f(s) has no zeros (or singularities like branch point and branch cut in dealing with multi-valued functions) on it. Then it is concluded from the Residue theorem that

^ f^Jds = M, (21)

2ni JY f (s)

where M is equal to the total number of the roots of f(s) = 0 inside y. Clearly, if the contour y is considered such that all zeros of f(s) lie inside it then we have M = mj. Equation (21) can be used to calculate the

Figure 3. The contour y considered on the first Riemann sheet for stability testing of the fractional-delay system under consideration. A system with characteristic equation A(s) = 0 (as defined in (9)) is stable if and only if it does not have any roots inside y.

number of zeros of the given function f(s) inside the desired contour y (which, of course, should have the above-mentioned properties). For this purpose, we can simply use a numerical integration technique to evaluate the integral in the right hand side of (21) for the given contour Y and function f.

According to the above discussions, by setting f(s) equal to A(s) and y equal to the border of the region of instability (which is equal to the closed right half-plane of the first Riemann sheet) the value obtained for M from (21) will be equal to the number of unstable roots of the characteristic equation. In the following we consider the contour y as shown in Fig. 3 and f(s) = A(s) (where A(s) is as defined in (9)) and then simplify the integral in the left hand side of (21) to arrive at a more effective formula for stability testing of the fractional-delay system under consideration (clearly, the system is stable if and only if M = 0). Note that the very small semicircle in Fig. 3 is used to avoid the branch-point located at the origin.

According to (21) and Fig. 3 we can write

M ^f^ds =±([ +¡ + 1 ). (22) 2nljY A(s) 2m\JCj+C3 JC2 JCJ ()

In (22), the Integral fc +C Is calculated as

I — f w + rA'(iw):

JC1+C3 J™ A(iw) J-

lci+c3

f™ A'(LU).. f™ A'(-LU). ...

i aUidw + i A—J)(-i)dw (24)

■L r - L i™ I )' dw, (25)

I A(Lw) J£ \ A(Lw) j ( )

which yields

A'(Lw) A(Lw)

JC1+C3 Je

The integral C in (22) is calculated as

[ = im f2 eie^de

Jc2 e^jn ■

A(eeie

i - ^\e Are

Jn A(eeie

Le de.

In the above equation Ume^0 A(eeiS) is equal to a nonzero constant (else, the characteristic function has a strong singularity at the origin and the corresponding system is unstable) and A'(eeiS) ~ Ken as e ^ 0, where K and n > —1 are two constants. Hence, assuming A(0) = 0, f tends to zero as e ^ 0. Finally, f in (22) is calculated as

f — Um f

Jc4 R^™J-.

-—i:

A'(Reie A(Reie

RLeie de

— iann.

^ IT^R

A(Reie)

Le de (30) (31)

the integrand at different points on the w axis. Hence, the numerical integration algorithm may halt if the integrand becomes singular at the origin (which is the case if, for example, we have 0 < a1 < 1 in (11)). In practice, in order to determine the number of unstable poles of the given fractional-delay transfer function we can consider the lower and upper bound of the integral in (32) equal to sufficiently small and big positive numbers, respectively. Numerical simulations performed by author showed that in most cases the MATLAB function quadl (as well as quadgk) can evaluate the integral in (32) with a high degree of accuracy. Note that quadl implements a high order method using an adaptive Gauss/Lobatto quadrature rule

[28] and quadgk implements adaptive quadrature based on a Gauss-Kronrod pair (15th and 7th order formulas)

[29]. Some numerical examples are presented in the next section.

3. Numerical examples

In the following we study the application of (32) for stability testing of some fractional-delay systems. In each case, the impulse response of the corresponding system is also plotted to verify the correctness of the result. The method used in this paper to calculate the impulse response of the given fractional-order system is based on the formula proposed in [30] for numerical inversion of Laplace transforms. In this method the impulse response of the given fractional-order system is approximated by numerical inversion of its transfer function. The MATLAB code of this method, invlap.m, can freely be downloaded from http://www.mathworks.com/matlabcentral/fileexchange/. Most of the following examples have already been studied by author in [26].

Example 1. Consider a system with characteristic equation

(See (9) and (11) for the definition of an.) Substitution of (26) and (31) in (22) and considering the fact that f = 0 results in

M — H -

1 f™ ^

n Je—0+

A'(Lw) A(Lw)

A1 (s) = (sn/2+1)(sn/3 + 1) (33)

= s5n/e + sn/2 + sn/3 +1=0. (34)

The roots of this equation can be calculated analytically, which are as the following:

where M is equal to the number of unstable poles of a system with characteristic equation A(s) = 0 as defined in (9).

Equation (32) is the main result of this paper. It should be noted that the value of e in (32) cannot in general be considered exactly equal to zero. This is because of the fact that the numerical integration technique used to evaluate the integral in (32) performs this task by evaluating

sk1— ei2(2k1+1\ k1 G Z,

s — e'3(2k2+1)i k2 G Z.

As it is observed, the characteristic equation given in (34) has Infinite many roots which are distributed on a Riemann surface with infinite number of Riemann sheets. It

is concluded from (35) and (36) that (34) has four roots on the first Riemann sheet which are e±j2 and e±j3, and none of them are located in the right half-plane (recall that all roots whose phase angle lies in the range [-n, n) belong to the first Riemann sheet).

Comparing (34) with (9) and (11) yields an = 5n/6 (note that (34) has no delay terms). Application of (32) assuming that the lower and upper bound of integral in (32) are equal to 0 and 1000, respectively, yields M = 2.3300 x 10-4 which is consistent with the above-mentioned analytical result. Figure 4 shows the impulse response of a system with transfer function

Hi (s) =

A1 (s) s5n/6 + sn/2 + sn/3 + 1 .

It can be observed in this figure that the impulse response of the system is absolutely summable, as expected.

Example 2. Stability of a system with fractional-delay characteristic equation:

A2(s) = s + K (y/s + 1).

is studied in [31 ] and it is especially shown that it is stable for K < 21.51 and unstable for K > 21.51. Application of (32) assuming that K = 21, an = 1, and the lower and upper bound of integral are equal to 0 and 500, respectively, yields M = 3.4227 x 10-9, which implies the stability of system as it is expected. Figure 5 shows the impulse response of a system with transfer function

H2(s) =

A2(s) s + 21 (+ 1),

As it can be observed, the impulse response is absolutely summable, as it is expected. Repeating the above procedure with K = 22 yields M = 2.0174, which means that in this case the system has two unstable poles. This result is also consistent with the one presented in [26].

Example 3. It is shown in [32] that a system with characteristic equation

A3(s) = s15 - 1.5s + 4s0 5 + 8 - 1.5se-Ts = 0, (40)

is stable for the values of t e (0.99830,1.57079) and unstable for other values of t. It is also shown by author in [26] that this system has two unstable poles for t = 0.99. Application of (32) assuming t = 1 and considering the fact that here we have an = 1.5 yields M = 0.0082 (the lower and upper bound of integral are considered equal to 0 and 500, respectively). As it is observed, the result obtained by using the proposed method is fairly close to zero. Figure 6 shows the impulse response of a system with transfer function:

H3(s) =

A3(s) = s15 - 1.5s + 4s0 5 +8 - 1,5se-s'

As it can be observed in this figure, the impulse response of the system is absolutely summable, and consequently, the corresponding system is stable. In this example the application of (32) assuming t = 0.99 yields 1.9994 which is consistent with the result presented in [26].

Example 4. It is shown in [19] (by applying Lambert W function) that a system with the following characteristic equation

A4(s) = s5/6 + (s1/2 + s1/3) e-05s + e-s = 0, (42)

10 t(s)

is stable. Clearly, here we have an = 5/6. Application of (32) (assuming that the lower and upper bound of integral are equal to 0 and 100, respectively) leads to M = 0.0290, which implies the stability of system. Figure 7 shows the impulse response of a system with transfer function

H4(s) =

A4(s) s5/6 + (s1/2 + s1/3) e-0 5s + e-

It can be observed that the impulse response is absolutely summable and consequently, the system is stable, as it is expected.

4. Conclusion

Figure 6. Impulse response of a system with transfer function (41).

Figure 7. Impulse response of a system with transfer function (43).

In this paper, bounded-input bounded-output stability of a very general class of linear time-invariant systems is studied in the frequency domain. The main result of this paper is the formula proposed for stability testing of LTI systems whose characteristic equation contain, in general, both the fractional powers of s and the exponentials of fractional powers of s. The proposed formula determines the number of unstable poles of the given fractional-delay transfer function by applying Rouche's theorem. More precisely, it is shown in this paper that the Rouche's theorem can be used to arrive at a simple formula for counting the number of roots of a multi-valued characteristic equation in the right half-plane of the first Riemann sheet. Four numerical examples are presented which support the efficiency of the proposed formula.

Appendix A: THE MATLAB CODE OF NUMERICAL EXAMPLES

The following are the MATLAB functions used in examples of this paper to calculate the number of unstable poles. Note that the numbers in the left hand-side of lines are to hint and should not be typed in practice. Example 1:

1: function y = f(w) 2: s = i*w;

3: y=real(((5*pi/6)*s.~(5*pi/6-1)+(pi/2)*s.~(pi/2-1)+(pi/3)*s.~(pi/3-1)) ./(s.~(5*pi/6)+s.~(pi/2)+s.~(pi/3)+1));

Example 2:

1: function y = f(w) 2: s = i*w;

3: y=real((1-22*0.5*exp(-s.~0.5))./(s+22*(s.~0.5).*exp(-s.~0.5) +21*exp(-s.~0.5)));

Example 3:

1: function y = f(w)

2: s = i*w; 3: tau = 0.99;

4: y=real((1.5*s.~0.5-1.5+2*s.~(-0.5)-1.5*(exp(-tau*s)-tau*s.*exp(-tau*s))) ./(s.~1.5-1.5*s+4*s.~0.5+8-1.5*s.*exp(-tau*s)));

Example 4:

1: function y = f(w) 2: s = i*w;

3: y=real(((5/6)*s.~(-1/6)+((1/2)*s.~(-1/2)+(1/3)*s.~(-2/3)).*exp(-.5*s)+

(s.~(1/2)+s.~(1/3))*(-.5).*exp(-.5*s)-exp(-s))

./(s.~(5/6)+(s.~(1/2)+s.~(1/3)).*exp(-.5*s)+exp(-s)));

In each case the number of unstable roots can easily be obtained by using the quadl or quadgk function of MATLAB. For example, since in Example 1 we have an = 5n/6 we can simply use the MATLAB command M=5*pi/6/2-quadl(@f,0,1000)/pi to determine the number of its unstable poles.

References

[1] R. HUfer, Applications of Fractional Calculus In Physics (World Scientific, New Jersey, 2000)

[2] R. L. Bagley, P. Torvik, J. Appl. Mech. 51, 294 (1984)

[3] V. G. Jenson, G. V. Jeffreys, Mathematical Methods in Chemical Engineering, 2nd edition (Academic Press, New York, 1977)

[4] N. Nakris, M. C. Constantinous, J. Struct. Eng. 117, 2708 (1991)

[5] G. Haneczok, M. Weller, Mat. Sci. Eng. A-Struct. 370, 209 (2004)

[6] I. Podlubny, Fractional Differential Equations (Academic Press, New York, 1999)

[7] B. Onaral, H. H. Sun, H. P. Schwan, In: Proceedigs of the 10th Northeast Bioengineering Conference, Mar. 1982, Hanover, NH, 46

[8] D. W. Davidson, R. H. Cole, Chem. Phys. 18, 1414 (1950)

[9] R. J. Schwartz, B. Friedland, Linear Systems (McGraw-Hill, New York, 1965)

[10] R. F. Curtain, H.J. Zwart, An Introduction to Infinite Dimensional Linear Systems Theory (Springer, Berlin, 1995)

[11] H. Zwart, Syst. Control Lett. 52, 247 (2004)

[12] T. Helie, D. Matignon, Signal Process. 86, 2516 (2006)

[13] D. Baleanu, K. Diethelm, E. Scalas, J.J. Trujillo, Fractional Calculus: Models and Numerical Methods (Series on Complexity, Nonlinearity and Chaos), (World Scientific, Singapore, 2012)

[14] K. Balachandran, J. Kokila, J. J. Trujillo, Comput.

Math. Appl. 64, 3037 (2012)

[15] S. Abbas, M. Benchohra, J. Nieto, Adv. Diff. Equ. 2011, 379876 (2011)

[16] A. Babakhani, D. Baleanu, Abstr. Appl. Anal. 2011, 391971 (2011)

[17] I. Podlubny, IEEE T. Automat. Contr. 44, 208 (1999)

[18] F. Merrikh-Bayat, Can. J. Chem. Eng. 90,1400 (2012)

[19] C. Hwang, Y.-C. Cheng, Automatica 42, 825 (2006)

[20] F. Jarad, T. Abdeljawad, D. Baleanu, Nonlinear AnalReal 14, 780 (2013)

[21] F. Jarad, T. Abdeljawad, D. Baleanu, K. Bigen, Abstr. Appl. Anal. 2012, 476581 (2012)

[22] H. Delavari, D. Baleanu, J. Sadati, Nonlinear Dynam. 67, 2433 (2012)

[23] E. Kaslik, S. Sivasundaram, J. Comput. Appl. Math. 236, 4027 (2012)

[24] D. Matignon, ESAIM: Proceedings 5, 145 (1998)

[25] M. Ikeda, S. Takahashi, Electron. Comm. Jpn. 160, 41 (1977)

[26] F. Merrikh-Bayat, M. Karimi-Ghartemani, ISA T. 48, 32 (2009)

[27] F. Merrikh-Bayat, M. Karimi-Ghartemani, Math. Probl. Eng. 2008, D0I:10.1155/2008/419046 (2008)

[28] W. Gander, W. Gautschi, BIT 40, 84 (2000)

[29] L.F. Shampine, J. Comput. Appl. Math. 211,131 (2008)

[30] J. Valsa, L. Brancik, Int. J. Numer. Model El. 11, 153 (1998)

[31] N. Ozturk, A. Uraz, IEEE T. Automat. Contr. 29, 368

(1984)

[32] N. Ozturk, A. Uraz, IEEE T. Circuits Syst. 32, 393

(1985)