Tr 5

ISA Transactions 48 (2009) 32–37

Contents lists available at ScienceDirect

ISA Transactions journal homepage: www.elsevier.com/locate/isatrans

An efficient numerical algorithm for stability testing of fractional-delay systems Farshad Merrikh-Bayat ∗ , Masoud Karimi-Ghartemani Sharif University of Technology, Tehran, Iran

article

info

Article history: Received 11 June 2008 Received in revised form 13 September 2008 Accepted 13 October 2008 Available online 11 November 2008 Keywords: Fractional delay equation Stability analysis Numerical algorithm Fractional-order system

a b s t r a c t This paper presents a numerical algorithm for BIBO stability testing of a certain class of the so-called fractional-delay systems. The characteristic function of the systems under consideration is a multi-valued function of the Laplace variable s which is defined on a Riemann surface with finite number of Riemann sheets where the origin is a branch point. The stability analysis of such systems is not straightforward because there is no universally applicable analytical method to find the roots of the characteristic equation on the right half-plane of the first Riemann sheet. The proposed method is based on the Rouche’s theorem which provides the number of the zeros of a given function in a given simple closed contour. One advantage of the proposed method over previous works is that it gives the number and the location of the unstable poles. The algorithm has a reliable result which is illustrated by several examples. © 2008 ISA. Published by Elsevier Ltd. All rights reserved.

1. Introduction

those described by the transfer function

A large number of circuits and systems have distributed parameters and/or delay elements. Thermal processes, hole diffusion of transistors, electromagnetic devices, and transmission lines are typical examples. The mathematical models in the Laplace transform domain for √ these systems commonly appear as a ratio of two polynomials in s in combination with the fractional-delay √ term exp(− s) [1,2]. For example, H (s) =

tanh

√

√

s

,

s

(1)

appears in a boundary controlled and observed diffusion process in a bounded domain [2]. As another example, the transfer function cosh

H (s) = √

√

sx0

s sinh




√ , s

0 < x0 < 1,

(2)

corresponds to the heat equation with Neumann boundary control [2]. Other examples of a similar type can be found in Curtain and Zwart [2]. The presence of a fractional-order time delay in a feedback control system leads to a characteristic equation with infinity many isolated roots which is often called the fractionaldelay equation. In the literature, the fractional-delay systems are

q0 (s) + G(s) = p0 (s) +

n2 P k=1 n1

P

r

qk (s)e−βk s

, pk (s)e−γk s

k =1

where r is a real number such that 0 < r < 1, the pk are of the Pl Pm form jk=0 aj sαj with αj ∈ R+ and the qk are of the form j=k0 bj sδj with δj ∈ R+ . The function G as defined in (3) is a multi-valued function of the Laplace variable s the domain of definition of which is, in general, a Riemann surface with infinite number of Riemann sheets [3]. In the rest of this paper, we restrict our studies to a special case of (3) in which the powers of the Laplace variable s are positive rational numbers. This is not a considerable loss of generality as it covers almost all practical studies [2]. It is a well recognized fact that a fractional-delay system has, in general, infinite number of poles and zeros unless it can be expressed in the special form of G(s) =

b0 smα + b1 s(m−1)α + · · · + bm snα + a1 s(n−1)α + · · · + an

Corresponding author. Tel.: +98 21 66165982. E-mail addresses: [email protected] (F. Merrikh-Bayat), [email protected] (M. Karimi-Ghartemani).

,

(4)

where m, n ∈ N, α ∈ (0, 1] and ai , bj ∈ R for i = 1, . . . , n and j = 1, . . . , m. The function G as defined in (4) is often called the fractional system of commensurate order. It is a well-known fact [4] that (4) is BIBO stable if and only if all roots of the equation

w n + a1 w n − 1 + · · · + an = 0 , ∗

(3)

r

(5)

lie in the sector

| arg(w)| >

0019-0578/$ – see front matter © 2008 ISA. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.isatra.2008.10.003

π 2

α,

(6)

F. Merrikh-Bayat, M. Karimi-Ghartemani / ISA Transactions 48 (2009) 32–37

in w plane where w = sα . In the field of infinite-dimensional fractional-delay systems most studies are concerned about the stability of √ a class of distributed systems whose transfer functions involve s and/or √ e− s [5,6]. Recently, several authors have studied the stability of a class of more general fractional-delay systems. For example, Hotzel [7] presented the stability conditions for fractional-delay systems with the characteristic equation (asα + b) + (csα + d)e−ρ s = 0, and Chen and Moore [8] analyzed the stability of a class of fractional-delay systems whose characteristic function can be represented as the product of factors of the form (as + b)r ecs + d, where the parameters a, b, c, d, and r are real numbers. See Hwang and Cheng [9] for a survey on this subject. In the literature few theorems are available for stability testing of fractional-delay systems. Almost all of these theorems are based on the locations of the transfer function poles [10,11] and since there is no universally applicable analytical method for solving fractional-delay equations in s domain, the numerical approach is commonly used. The most famous numerical approach for studying the stability of a certain class of the fractional-delay systems is the so-called Lambert W function which provides the solutions of the equation ses = x where x ∈ C is an arbitrarily chosen constant [12]. There are also few other numerical methods that are based on the Cauchy integral theorem [3] which states that if f (s) is nonsingular on and inside the contour Γ , then the complex integral I =

I

1 2π i

Γ

f (s)ds,

(7)

vanishes. Otherwise, it is evaluated to be the sum of the residues of the function f (s) at the poles in Γ . So, in order to examine the stability of a system with the characteristic function f (s) one may try to evaluate (7) and then compare its value with zero, where Γ is the border of the region of instability. But since it may occur that the sum of the residues of f (s) at all poles in Γ is zero, the I = 0 is only a necessary condition for f (s) to be nonsingular in Γ . To overcome this difficulty, Hwang and Cheng [9] proposed a numerical algorithm which uses the modified explanation of (7) in the form of

Z

i∞

Jk = −i∞

f (s)

(s + h1 + ih2 )k f (ih2 )

ds.

(8)

In (8), h1 > 0 and h2 are randomly chosen real constants lying in a specified interval and k is a positive integer. The randomness of the parameters h1 and h2 makes the probability of the zero sum of the residues of all poles of the integrand being practically zero. Hence, the stability of a given fractional-delay system can be achieved by evaluating the integral Jk and comparing its value with zero. As stated in [9], this approach should be carried out for different values of h1 , h2 , and k in order to obtain a reliable result. Anyway, the algorithm proposed by Hwang and Cheng [9] provides no idea about the number and the location of unstable poles. In this paper, we will introduce another numerical algorithm for stability testing of the fractional-delay systems which does not suffer from the above drawbacks. The proposed method is applicable to all systems the transfer function of which are in the form of (3), assuming that all powers of s are positive rational numbers. This method is based on Rouche’s theorem which computes the number of the zeros of a function in a given contour. One advantage of the proposed method is that it enables one to determine the number of unstable poles. It also computes the location of unstable poles with an arbitrary precision. The rest of this paper is arranged as follows. Rouche’s theorem is reviewed in Section 2. The proposed numerical algorithm and numerical examples are presented in Sections 3 and 4, respectively. Finally, Section 5 concludes the paper.

33

Fig. 1. The suitable contour for stability testing in classical case.

2. Review of Rouche’s theorem Consider the function f : C → C which has zeros of orders m1 , . . . , mk respectively at α1 , . . . , αk and poles of orders p1 , . . . , pn respectively at β1 , . . . , βn . Obviously the function f can be written as f (s) = g (s)

(s − α1 )m1 . . . (s − αk )mk , (s − β1 )p1 . . . (s − βn )pn

(9)

where g (s) has neither (finite) pole nor (finite) zero. By taking the logarithm from both sides of (9) we obtain ln f (s) = ln g (s) + m1 ln(s − α1 ) + · · · + mk ln(s − αk )

− p1 ln(s − β1 ) − · · · − pn ln(s − βn ).

(10)

Derivation with respect to s yields f 0 ( s) f ( s)

=

g 0 (s) g (s)

−

+

m1 s − α1

p1 s − β1

+ ··· +

− ··· −

mk s − αk

pn s − βn

.

(11)

Now let γ be a simple, closed, counterclockwise curve such that f (s) has neither pole nor zero on it. Then it is concluded from the Residue theorem that 1

f 0 ( s)

I

2π i

γ

f (s)

ds =

k X j =1

mj −

n X

pj = M − P ,

(12)

j=1

considering the fact that g (s) is analytic on and inside γ . Eq. (12) is an explanation of Rouche’s theorem. If f (s) has no pole inside γ then it is concluded from (12) that 1 2π i

f 0 ( s)

I γ

f (s)

ds = M .

(13)

In classical case (i.e. in dealing with single-valued transfer functions), the BIBO stability of a given system can be examined by choosing f (s) equal to the characteristic function where γ is the semicircle shown in Fig. 1. In the next section, a modified explanation of (13) will be used to find the number of the zeros of the multi-valued characteristic function ∆(s) inside the contour γ which is chosen equal to the border of instability. 3. The proposed numerical algorithm Let ∆(s) be the characteristic function of a (multi-valued) fractional-delay system which can be transformed to the single1

e(w) with the change of variable w = s v . Then, valued function ∆

34

F. Merrikh-Bayat, M. Karimi-Ghartemani / ISA Transactions 48 (2009) 32–37

Note that in the above equation the order of integration and limitation can be changed since the integral is uniformly convergent. In order to evaluate the second integral in (15) when R → ∞ and α = 2πv , the change of variable r = tan x is applied. By using this change of variable the improper integral is transformed to a proper one, which can be numerically evaluated by software. It then follows that

(

∞

Z

= 0

e0 (rei 2πv ) i π ∆ e 2v e(rei 2πv ) ∆

)

π

Z

(

2

dr =

= 0

e0 (tan xei 2πv ) i π ∆ e 2v e(tan xei 2πv ) ∆

× (1 + tan2 x)dx.

)

(20)

Substituting (19) and (20) in (15) yields Fig. 2. The contour used for stability testing of fractional delay systems.

e(w) = 0 inside the contour γ (as the number of the roots of ∆ shown in Fig. 2) is calculated from (13) as M (R, α) =

e0 (re−iα ) −iα ∆ e dr + e −i α ) 2π i r =0 ∆(re  Z 0 e 0 iα ∆ (re ) iα e dr , + e iα r =R ∆(re ) Z

1

R

Z

α θ=−α

e0 (Reiθ ) iθ ∆ Rie dθ e(Reiθ ) ∆ (14)

which after some algebra yields α

) e0 (Reiθ ) ∆ M (R, α) = R R ei θ d θ e(Reiθ ) π ∆ 0 ) ! Z R ( e0 iα ∆ (re ) iα − = e dr , e(reiα ) ∆ 0 Z

1

(

(15)

where R{.} and ={.} stand for the real and the imaginary part, respectively. Note that for reasons of simplicity it is assumed that e(w) = 0 has no roots on γ . It is a well-known fact that the ∆ characteristic equation ∆(s) = 0 corresponds to a BIBO stable e(w) = 0 has no roots in the sector system [10,11] if and only if ∆

| arg(w)| ≤

π 2v

.

(16)

√

√

For example, a system with the characteristic equation s + e s = 0 is BIBO stable if and only if all roots of the equation w + ew = 0 lie in the sector | arg(w)| > π /4. Hence, in order to find the number of unstable roots of the (multi-valued) characteristic function ∆(s), first we transform it e(w) with the change of variable to the single-valued function ∆

w = s v and then calculate M (∞, 2πv ) from (15). To proceed, we e(w) in the general form of consider ∆ 1

e(w) = e ∆ P0 (w) +

n X

βi e Pi (w)e−αi w ,

(17)

i=1

where αi ∈ R , βi ∈ {1, . . . , v}, and Pi are polynomials in w , in particular, +

e P0 (w) = a0 w p + · · · + ap−1 w + ap ,

(18)

where p ∈ N and ai ∈ R for i = 0, . . . , p. Note that the e(w) (as defined in (17)) covers almost all characteristic function ∆ real-world applications [2]. Now, for R → ∞ and α = 2πv , the first integral in (15) is evaluated as π

 e0 (Reiθ ) iθ ∆ Re dθ = e(Reiθ ) R→∞ 0 ∆   Z π 2v a0 p(Reiθ )p−1 + · · · + ap−1 iθ pπ lim R Re dθ = . iθ )p + · · · + a R →∞ a ( Re 2v 0 p 0 Z

lim

2v

) π Z π ( e0  2 p π 1 ∆ (tan xei 2v ) i π 2 v = M ∞, − = e e(tan xei 2πv ) 2v 2v π 0 ∆ × (1 + tan2 x)dx.

Eq. (21) is the key formula of this paper which can be used to find the number of unstable zeros of the characteristic function (17). In some applications the exact location of the roots of the (multi-valued) equation ∆(s) = 0 on the first Riemann sheet is e(w) = 0 in required. In order to find the modulus of the roots of ∆ the sector | arg w| < π /v (which corresponds to the first Riemann sheet), first we let α = π /v and then calculate M (R, π /v) for several values of R from (15). Then, the distance of the roots from the origin can be investigated by plotting M (R, π /v) versus R. For instance, let R1 and R2 (> R1 ) be two successive values of R such that M (R2 ) ≈ M (R1 ) + k,

e(w) = 0 using the of ∆(s) = 0 can be obtained from the roots of ∆ map s = w v . There is another point that should be noted here and that is the proposed method can also be used to examine the stability of more complicated systems. For example, consider a system with √ 3 the characteristic function ∆(s) = 1 + s2 + s + 1 + e−s . In order to find the number of unstable poles first we should apply a change of variable to make the characteristic function single valued. For √ 3 2 example, we can use w = s + s + 1 + e−s which results in e(w) = 1 + w. Now, according to (13), the number of unstable ∆ roots, M, can be calculated as M =

(19)

(22)

for some k ∈ N (this means that there is an integer jump in the plot of M versus R). In this case, (22) implies that there are k roots in the region defined by R1 < |w| < R2 and | arg w| < π /v . In the proposed method, the values of R1 and R2 must be chosen to be close so that (22) holds for k = 1 (or k = 2, if k = 1 is not possible, i.e. if there is a pair of complex conjugate roots in the contour). In brief, in order to determine the modulus of the roots e(w) = 0, plot M (R, π /v) versus R. The values of R, for which of ∆ there is an integer jump in this plot, correspond to the modulus of the roots. After determining the distance of the roots from the origin, their phase angle can also be determined. Let R1 , . . . , Rn be the distance e(w) = 0 from the origin where by assumption of the roots of ∆ e(w) = 0 such that R1 < · · · < Rn . Assume that wi is a root of ∆ |wi | = Ri . In order to determine arg(wi ), plot M (R0 , α) versus α where R0 ∈ (Ri , Ri+1 ) is an arbitrarily chosen constant and 0 < α < π /v . The values of α for which there is an integer jump in this plot correspond to the phase angle of wi . Finally, the roots



R

(21)

1 2π i

I Γ

I e0 (w) ∆ 1 dw = e(w) 2π i ∆

1 Γ

w+1

dw,

where Γ is the projection of the semicircle shown in Fig. 1 under √ 3 2 the map s 7→ w = s + s + 1 + e−s which is shown in Fig. 3.

F. Merrikh-Bayat, M. Karimi-Ghartemani / ISA Transactions 48 (2009) 32–37

35

Fig. 3. The contour of integration√for studying the stability of a system with 3 characteristic function ∆(s) = 1 + s2 + s + 1 + e−s .

e0 (w)/∆ e(w) In this example it is easily seen that M = 0 since ∆ has no pole inside Γ . In general, it is not straightforward to apply the above method to find the number of unstable poles because the contour Γ is not easily parameterized and consequently, the evaluation of

H

Γ

Fig. 4. The plot of M (R, π/100) versus R corresponding to Example 1.

e0 (w)/∆ e(w)dw is a tricky task. ∆

4. Numerical examples In this section, several examples are provided to show the effectiveness of the proposed algorithm. A convenient approach Rb to evaluate the proper integral a g (x)dx is to define I (b) :=

Rb a

g (x)dx and then solve the initial-value problem

dI (x) dx

= g (x),

I (a) = 0.

(23)

In the following examples, this approach is used to evaluate the integrals in (15) and (21). The Matlab function ode23 which is an implementation of an explicit Runge–Kutta (2,3) pair of Bogacki and Shampine [13] is found to be accurate enough for this purpose. Example 1. Consider the fractional-order transfer function proposed by Podlubny et al. [14] for a heating furnace: P ( s) =

1 14994s1.31 + 6009.5s0.97 + 1.69

,

(24)

which is intuitively stable. In the following, we will determine the location of the (stable) poles of P (s) on the first Riemann sheet using the proposed algorithm. The transfer function P (s) can be represented in the equivalent form P ( s) =

1 14994s

131 100

97

+ 6009.5s 100 + 1.69

,

(25)

which corresponds to

e P (w) =

1 14994w 131 + 6009.5w 97 + 1.69

.

(26)

It is easily seen that P (s) has 131 poles distributed on 100 Riemann sheets. In order to find the location of the poles of e P (w) on the first e(w) = 14994w131 + 6009.5w 97 + 1.69 Riemann sheet we let ∆ and then calculate M (R, π /100) from (15) for several values of e(w) with (17) and (18) leads to p = R (note that comparing ∆ 131). Fig. 4 shows the plot of M (R, π /100) versus R. As it is seen, M (0.918, π /100) ≈ 0 and M (0.919, π /100) ≈ 2 which implies

Fig. 5. The plot of M (0.92, α) versus α corresponding to Example 1.

that there are two complex conjugate roots in the region defined by

| arg w| < π /100,

0.918 < |w| < 0.919,

(27)

i.e. the modulus of the roots is between 0.918 and 0.919. By plotting M (R, π /100) versus R ∈ [0.918, 0.919] the exact modulus of the roots is found to be R1 ≈ 0.91848. The phase angle of the roots can be obtained by plotting M (0.92, α) versus α , which is shown in Fig. 5. It is concluded from this figure that the absolute value of the phase angles is between 0.03122 rad and 0.03125 rad. Hence, the roots (in the complex w plane) are calculated as w1,2 ≈ 0.91848e±i0.0312 . The location of the roots on the first Riemann sheet of the s Riemann surface can be obtained using the map s = w 100 , which results in s1,2 ≈ 2.0275 × 10−4 e±i3.12 . As expected, both poles are located on the left half-plane of the first Riemann sheet. Note that the number of unstable poles can also be calculated from (21) which leads to M (∞, π /200) = −0.0001 and can be safely regarded as zero. Example 2. Consider the fractional-delay characteristic function √ √ √ ∆(s) = ( s)2 + K ( s + 1)e− s ,

(28)

where K is a positive real constant. It has been shown by Ozturk and Uraz [5] that ∆(s) corresponds to a stable system if K < 21.51

36

F. Merrikh-Bayat, M. Karimi-Ghartemani / ISA Transactions 48 (2009) 32–37

Fig. 6. The plot of M (R, π/4) versus R corresponding to Example 2 (K = 21).

Fig. 7. The plot of M (R, π/4) versus R corresponding to Example 2 (K = 22).

and has some roots in the sector of instability | arg(w)| ≤ π /4 if K > 21.51. In order to deal with the stability of the characteristic function (28) with the proposed method, the change of variable w = s1/v = s1/2 is applied which results in

e(w) = w 2 + K (w + 1)e−w . ∆

(29)

Comparing (29) with (17) and (18) leads to p = 2. Numerical evaluation of (21) for K = 22 and K = 21 results in M (∞, π /4) = 2.0000 and M (∞, π /4) = 2.1412 × 10−5 , respectively which means that (28) has two unstable poles for K = 22 and is stable for K = 21. Fig. 6 shows the plot of M (R, π /4) versus R for k = 21. As it is expected, there is no integer jump in this figure. In order e(w) = 0 when to find the modulus of the unstable roots of ∆ K = 22, M (R, π /4) must be calculated for several values of R. Fig. 7 shows the plot of M (R, π /4) versus R in this case. As viewed, M (3, π/4) ≈ 0 and M (3.1, π /4) ≈ 2 which implies that there are two (complex conjugate roots) in the region defined by 3 < |w| < 3.1,

| arg(w)| <

π 4

.

Fig. 8. The plot of M (3.5, α) versus α corresponding to Example 2 (K = 22).

(30)

In the same manner, by inspecting the values of R between 3 and 3.1 the exact modulus of the roots is found as R1 ≈ 3.085. The e(w) = 0 can be obtained by plotting phase angle of the roots of ∆ M (R0 , α) versus α , where R0 > R1 is an arbitrarily chosen constant. Fig. 8 shows the plot of M (3.5, α) versus α . According to this figure, there are two roots in the region defined by 0 < |w| < 3.5,

0.764 < | arg(w)| < 0.7954.

(31)

By inspecting the values of α in the interval [0.764, 0.7954] the e(w) = 0 is absolute value of the phase angle of the roots of ∆ found to be 0.782 rad. Therefore, the roots of (29) in the complex w plane are w1,2 ≈ 3.085e±i0.782 . The locations of the unstable roots of (28) can be calculated using the map s = w 2 , which results in s1,2 ≈ 9.5172e±i1.564 . As expected, both roots lie in the right halfplane of the first Riemann sheet.

proposed method. By using the change of variable w = s1/2 the corresponding characteristic function in w -plane is obtained as

e(w) = w3 − 1.5w 2 − 1.5w 2 e−τ w2 + 4w + 8. ∆ Now, the numerical evaluation of (21) for τ = 0.99 and τ = 1 results in M (∞, π /4) = 2.0082 and M (∞, π /4) = −0.0047, respectively which means that (32) is stable for τ = 1 and unstable for τ = 0.99. This result is completely consistent with the one presented by Ozturk and Uraz [6]. One can also use the proposed method to find the locations of the unstable poles which is not discussed here. Example 4. Consider the fractional-order characteristic function

∆(s) = s5/6 + (s1/2 + s1/3 )e−0.5s + e−s .

(33)

Example 3. Consider a fractional-delay system with the characteristic function

Hwang and Cheng [9] have obtained the roots of (33) via the Lambert W function and have shown that the corresponding system is stable. Clearly, in order to apply the proposed numerical algorithm, the change of variable w = s1/6 should be used which yields

√ √ √ √ ∆(s) = ( s)3 − 1.5( s)2 − 1.5( s)2 e−τ s + 4( s) + 8.

e(w) = w5 + (w 3 + w2 )e−0.5w6 + e−w6 . ∆

(32)

Ozturk and Uraz [6] have shown that this system has no pole in the sector of instability | arg(w)| ≤ π /4 for τ ∈ (0.99830, 1.57079). In the following, we study the stability properties of (32) using the

(34)

Putting (34) in (21) yields M (∞, π /12) = −0.0190 which can be safely regarded as zero. Hence, the characteristic function (33) corresponds to a stable system as expected.

F. Merrikh-Bayat, M. Karimi-Ghartemani / ISA Transactions 48 (2009) 32–37

5. Conclusions An effective numerical algorithm for stability testing of a large class of fractional delay systems is presented. The proposed numerical method is based on Rouche’s theorem which can be used to find the number of zeros of a given function in a simple closed contour. One advantage of the proposed method is that it also provides the number and the exact location of the stable and unstable poles. Four numerical examples are presented to confirm the effectiveness of the proposed method. References [1] Chen CF, Chiu RF. Evaluation of irrational and transcendental transfer functions via the fast fourier transform. International Journal of Electronics 1973;35: 267–76. [2] Curtain RF, Zwart HJ. An introduction to infinite-dimensional linear systems thoery. Text in applied mathematics, vol. 21. Berlin: Springer; 1995. [3] LePage WR. Complex variables and the Laplace transform for engineers. McGraw-Hill; 1961.

37

[4] Matignon D. Stability properties for generalized fractional differential systems, In: ESAIM: Proc. 5, 1998, p. 145–58. [5] Ozturk N, Uraz A. An analytic stability test for a certain class of distributed parameter systems with a distributed lag. IEEE Transactions on Automatic Control 1984;29(4):368–70. [6] Ozturk N, Uraz A. An analytic stability test for a certain class of distributed parameter systems with delay. IEEE Transactions on CAS 1985;32(4): 393–396. [7] Hotzel R. Some stability conditions for fractional delay systems. Journal of Mathematical Systems, Estimation, and Control 1998;8(4):1–19. [8] Chen YQ, Moore KL. Analytical stability bound for a class of delayed fractionalorder dynamical systems. Nonlinear Dynamics 2002;29(1–4):191–200. [9] Hwang C, Cheng Y-C. A numerical algorithm for stability testing of fractional delay systems. Automatica 2006;42:825–31. [10] Bonnet C, Partington JR. Stabilization of fractional exponential systems including delays. Kybernetika (Prague) 2001;37(3):345–53. [11] Bonnet C, Partington JR. Analysis of fractional delay systems of retarded and neutral type. Automatica 2002;38(7):1133–8. [12] Corless RM, Gonnet GH, Hare DEG, Jeffrey DJ, Knuth DE. On the Lambert W function. Advances in Computational Mathematics 1996;5:329–59. [13] Bogacki P, Shampine LF. A 3(2) pair of Runge–Kutta formulas. Applied Mathematics Letters 1989;2:1–9. [14] Podlubny I, Dorcak L, Kostial I. On fractional derivatives, fractional-order dynamic system and PIλ Dµ -controllers, In: Proc. of the 36th IEEE CDC, 1999.