Synthesis of fractional Laguerre basis for system approximation

Automatica 43 (2007) 1640 – 1648 www.elsevier.com/locate/automatica

Brief paper

Synthesis of fractional Laguerre basis for system approximation夡 M. Aoun a,b , R. Malti a,∗ , F. Levron c , A. Oustaloup a a IMS, UMR 5218 CNRS, Département LAPS, Université Bordeaux 1, 351 cours de la Libération, F 33405 Talence cedex, France b MACS, École Nationale d’Ingénieurs de Gabès, Rue Omar Ibn El Khattab, T 6029 Gabès, Tunisia c IMB, UMR 5251 CNRS, Université Bordeaux 1, 351 cours de la Libération, F 33405 Talence cedex, France

Received 6 April 2004; received in revised form 17 July 2006; accepted 19 February 2007

Abstract Fractional differentiation systems are characterized by the presence of non-exponential aperiodic multimodes. Although rational orthogonal bases can be used to model any L2 [0, ∞[ system, they fail to quickly capture the aperiodic multimode behavior with a limited number of terms. Hence, fractional orthogonal bases are expected to better approximate fractional models with fewer parameters. Intuitive reasoning could lead to simply extending the differentiation order of existing bases from integer to any positive real number. However, classical Laguerre, and by extension Kautz and generalized orthogonal basis functions, are divergent as soon as their differentiation order is non-integer. In this paper, the first fractional orthogonal basis is synthesized, extrapolating the definition of Laguerre functions to any fractional order derivative. Completeness of the new basis is demonstrated. Hence, a new class of fixed denominator models is provided for fractional system approximation and identification. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Orthonormal basis; Fractional differentiation; Laguerre function; System approximation; Identification

1. Introduction Over the last 15 years, orthogonal functions have been widely used in identification and control of linear systems; see for instance Malti, Ekongolo, and Ragot (1998), Ninness and Gustafsson (1997), Van den Hof, Heuberger, and Bokor (1995) and Wahlberg (1991) and their own references. The most popular orthogonal functions used in control engineering are Laguerre functions with a single pole, Kautz functions with two complex conjugate poles, and the generalized orthogonal basis (GOB) functions. The latter extend the first two to any number of real or complex conjugate poles. All these bases span Lebesgue space of squared integrable functions, provided 夡 This paper was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor Brett Ninness under the direction of Editor Torsten Söderström. ∗ Corresponding author. E-mail addresses: [email protected] (M. Aoun), [email protected] (R. Malti), [email protected] (F. Levron), [email protected] (A. Oustaloup).

0005-1098/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2007.02.013

the completeness condition is satisfied (Akçay & Ninness, 1999). They can hence approximate any linear system, including those with fractional differentiation, provided its impulse response belongs to L2 [0, ∞[. However, fractional systems, such as thermal (Battaglia, Le Lay, Batsale, Oustaloup, & Cois, 2000) or electrochemical (Darling & Newman, 1997) systems (see the special issue of Signal Processing (Ortigueira & Tenreiro Machado, 2006) for other fields of application), are characterized by the presence of two kinds of modes: exponential modes, as in rational systems, and aperiodic multimodes (Oustaloup, 1995) which decay polynomially. Using rational orthogonal bases to approximate a fractional L2 [0, ∞[ system leads to a high truncation order to capture the non-exponential mode’s behavior. An intuitive approach to extend Laguerre functions to fractional differentiation orders led El Sayed (1999) to simply allow their differentiation orders to be positive real numbers without any precaution. However, Abbott (2000), commenting on El-Sayed’s work, has since proved that classical Laguerre functions are divergent whenever their differentiation order is non-integer.

M. Aoun et al. / Automatica 43 (2007) 1640 – 1648

Nomenclature (.) . nm T

N R R+ R+∗ C

C+ C+ L2 [0, ∞[

transpose of a matrix floor operator Dirac delta function equals 1 if n = m, and 0 otherwise set of non-negative integers field of real numbers set of non-negative real numbers set of strictly positive real numbers field of complex numbers

H2 (C+ )

(m) C

1641

open right half-plane {s ∈ C : Re(s) > 0} closed right half-plane {s ∈ C : Re(s) 0} Lebesgue space of squared integrable functions in [0, ∞[ Hardy space of functions F, analytical on C+ andcontinuous on C+ such that F 2 = ∞ (1/2) −∞ |f (x + jy)|2 dy < ∞ Euler’s function, defined as (m) =  ∞ −x gamma m−1 dx∀m ∈ R\{0} e x 0 Cholesky decomposition

In this paper, a fractional orthonormal basis is synthesized by applying the Gram–Schmidt orthogonalization procedure. The objective is twofold: to extrapolate Laguerre functions to any fractional derivative, and to obtain convergent functions forming a complete orthonormal basis in H2 (C+ ). Hence, a new class of fixed denominator models is provided for fractional differentiation systems. The new orthonormal basis is governed by two tuning parameters: an eigenvalue and a commensurable order  ∈]0, 2[ (all differentiation orders are multiples of ). Laguerre functions correspond to the special case  = 1. This is the first working fractional orthonormal basis ever synthesized.

F (s) is called a commensurable transfer function at order  ∈ R+∗ , iff (i , j ) are commensurable1 with . This  is called a commensurable order and the biggest number is always chosen. Time-domain simulation of fractional transfer functions (3) is explained in Aoun (2005), Aoun, Malti, Levron, and Oustaloup (2004) and Poinot and Trigeassou (2003).

2. Mathematical background

 K  sk ts k sin() ∞ x  e−tx e + dx,   x 2 − 2 x  cos() + 2 0 k=1

2.1. Definition of fractional differentiation Differentiation to not only integer but also non-integer orders was defined in the 19th century by Riemann and Liouville. This is usually termed fractional differentiation. The th fractional derivative of a continuous real function f (t) is defined as (Oldham & Spanier, 1974):  +1 1 d D f (t) = ( −  + 1) dt  t  f () × d . − 0 (t − ) 

(1)

The Laplace transform of D f (t), when f is relaxed at t = 0 (f (t) and all its derivatives equal 0 for all t < 0), is given by (Oldham & Spanier, 1974): L(D f (t)) = s  F (s).

2.2. Aperiodic multimodes as compared to exponential modes Taking the inverse Laplace transform of the following simple fractional transfer function 1/(s  − ) gives (Oustaloup, 1995) (4)

where sk , k=1, . . . , K, are the s-roots of s  − =0. The number of s-roots, K, depends on  as shown in Oustaloup (1995). Expression (4) shows the presence of exponential modes in the left part, as in rational transfer functions. Also aperiodic multimodes, which decay polynomially, are present in the right part (Oustaloup, 1995). When  is integer there is no aperiodic multimode as sin() = 0 is a factor of the right part. Although classical rational orthogonal functions span completely L2 [0, ∞[, they are less appropriate for a good approximation of the polynomially decaying behaviors when using a limited number of terms. Hence, fractional orthogonal bases are expected to better approximate fractional models with fewer parameters. 2.3. Stability condition

(2)

This result is coherent with the classical case when  is an integer. Consequently, it is easy to define a symbolic representation of a fractional dynamic system using the transfer function: mA i i=0 ai s F (s) = , (3) mB 1 + j =1 bj s j where ai and bj belong to R and i and j belong to R+ , for i = 0, 1, . . . , mA , and for j = 1, 2, . . . , mB . Both i and j form strictly increasing sequences.

Matignon (1998, Theorem 2.21 p. 150) has established the stability condition of any commensurable fractional transfer function (3). However, here is a revisited version of his theorem: Stability Theorem. A commensurable -order transfer function F (s) = S(s  ) = T (s  )/R(s  ), where T and R are two coprime polynomials, is BIBO stable iff 0<<2

(5)

1 Exactly divisible by the same unit [number] an integral number of times (The American Heritage䉷, 2000).

1642

M. Aoun et al. / Automatica 43 (2007) 1640 – 1648

and for every s ∈ C such that R(s) = 0

which corresponds to the L2 norm of the approximation error, according to the definition of the scalar product (8):

 | arg(s)| <  . 2

(6)

J = f (t) − fN (t)22 .

(14)

Stability conditions (5) and (6) will be necessary when choosing fractional generating functions for the synthesis of the new basis.

Minimizing J and taking advantage of the orthonormality, Fourier coefficients are obtained by computing the scalar product in time or frequency domains:

2.4. Fractional transfer functions belonging to H2 (C+ )

an = f, ln = F, Ln .

Contrary to rational systems, the stability condition does not guarantee that a fractional transfer function belongs to H2 (C+ ). The H2 norm of fractional systems is extensively studied in Malti, Aoun, Levron, and Oustaloup (2003), where it is demonstrated that a stable fractional transfer function as defined in (3), where conditions (5) and (6) are satisfied, belongs to H2 (C+ ) iff its relative degree is greater than 21 : mB − mA > 21 .

(7)

Condition (7) will be necessary when choosing fractional generating functions for the synthesis of the new basis. 2.5. Scalar product, orthogonality and Laguerre functions Laguerre, Kautz and GOB functions form complete orthonormal bases in L2 [0, ∞[, according to the usual definition of the scalar product:  ∞ ln , lm = ln (t)lm (t) dt = nm , (8) 0

where the reciprocal in the frequency domain is obtained by Plancherel’s theorem:  j∞ 1 Ln , Lm = Ln (j )Lm (j ) d = nm . (9) 2j −j

Any function f (t) ∈ L2 [0, ∞[, thus satisfying: f, f 1/2 = f 2 < ∞

(10)

can be written as a linear combination of these functions: F (s) =

∞ 

an Ln (s).

(11)

F (s) is the Laplace transform of f (t). Usually, (11) is truncated to a given order N which is justified by the convergence of Fourier coefficients as n tends to infinity. F (s) is hence approximated by the finite sum: N 

an Ln (s).

(12)

n=0

The Fourier coefficients are computed by minimizing the leastsquares criterion:  ∞ J= (f (t) − fN (t))2 dt (13) 0

As stated previously ln (t) can represent any set of orthonormal functions in L2 [0, ∞[. In the case where Laguerre functions are used, ln (t) is defined as √

e t dn (t n e−2 t ) . n! dt n It has the following Laplace transform:

ln (t) =

Ln (s) =

2

√

2

(s − )n (s + )n+1

(16)

.

(17)

Laguerre functions are well suited for modelling systems with a dominant time constant, because they have a single pole. Abbott (2000) has shown that Laguerre functions (16) are divergent as soon as n becomes non-integer:  ∞ (ln (t))2 dt = ∞, ∀n ∈ R+∗ \N. (18) 0

Therefore, the generalization of Laguerre functions by simply allowing their differentiation order to be real is not possible. 3. Construction of the orthonormal basis This section contains the main results of the paper. 3.1. Gram–Schmidt orthogonalization procedure Given an arbitrary set of functions {Fm }m=m0 ,m0 +1,...,M , where Fm ∈ H2 (C+ ) ∀m, the orthonormalized functions {Gm }m=m0 ,m0 +1,...,M are obtained, according to the Gram–Schmidt procedure, as a linear combination of functions {Fm }m=m0 ,m0 +1,...,M : G =  × F,

(19)

where

n=0

F (s) ≈ FN (s) =

(15)

G = [Gm0 (s) Gm0 +1 (s) · · · GM (s)]T ,

(20)

F = [Fm0 (s) Fm0 +1 (s) · · · FM (s)] ,

(21)

T

and  is a lower-triangular orthogonalization matrix obtained as follows. Since G is the vector of orthonormal functions: G, GT

⎡ G , G

m0 m0 ⎢ ,G

⎢ G = ⎢ m0 +1. m0 ⎣ .. GM , Gm0

Gm0 , Gm0 +1

.. . ··· GM , Gm0 +1

· · · Gm0 , GM ⎤ .. .. ⎥ . . ⎥ ⎥ =I, .. .. ⎦ . . · · · GM , GM

(22)

M. Aoun et al. / Automatica 43 (2007) 1640 – 1648

1643

where I denotes M × M identity matrix. Thus, using (19): G, GT =  F, FT T = I.

(23)

The solution of the previous quadratic form is given by Cholesky decomposition C of the inverse Gram-matrix F, FT :  = C( F, FT −1 ).

(24)

As a result, the functions of the orthonormal set are given by G = C( F, FT −1 ) × F.

(25)

Evaluating the Gram-matrix F, FT requires, first of all, generating functions to be defined. Fig. 1. Conformal mapping 1/(s  + ) of C+ onto  for 0 <  < 1. When 1 <  < 2,  intrudes in C\C+ .

3.2. Generating functions To extrapolate Laguerre functions to any fractional order derivative, the {Fm (s)}m  m0 m ∈ N set of generating functions is chosen where: 1 Fm (s) =  (26) (s + )m

Proof. To prove the completeness of the generating functions {Fm }m=m0 ,...,∞ in H2 (C+ ), it suffices to prove the completeness in a dense subset of H2 (C+ ), denoted H, since for a given ε > 0 and K ∈ H2 (C+ ), there exists an H ∈ H satisfying:

and

K − H 2 < ε.

∈ R+∗ ,

 ∈]0, 2[.

(27)

Both conditions in (27) stem from the Stability Theorem in 2.3, and with impulse response of (26) set to be real-valued. To guarantee that the generating function (26) belongs to H2 (C+ ), an additional condition is obtained by applying (7) on (26): m > 21 . Keeping in mind that m is integer, yields   1 mm0 = + 1. 2

(28)

(29)

Condition (29) defines the starting value m0 of the subscript m to orthogonalize the set of functions Fm (s) defined in (26). 3.3. Evaluating elements of the Gram-matrix F, FT

Each element Fh , Fm of the matrix F, FT is a scalar product of two generating functions (26):  ∞ 1 Fh (j )Fm (j ) d

Fh , Fm = 2 −∞ h  m  ∞ 1 1 1 = d . (30) 2 −∞ (j ) + (−j ) + Evaluating integral (30) is not an easy task because Fm (s) is a complex multivalued function when  = 1. Consequently, a plane cut is necessary. Integral (30) is computed in Appendix A. 3.4. Completeness of the fractional basis Completeness Theorem. Define Fm (s) as in (26) with conditions (27) satisfied. Then, the linear span of {Fm }m=m0 ,...,∞ , where m0 = 1/2 + 1, is dense in H2 (C+ ).

For a given positive integer N, there exists a dense subset H of H2 (C+ ) satisfying Garnett (2007, Corollary 3.3): (i) H (s) is continuous on C+ and infinitely differentiable on the imaginary axis for all H ∈ H. (ii) lim|s|→∞ |s|N |H (s)| = 0, ∀s ∈ C+ . For a given H ∈ H, let G(s) = (s  + )m0 H (s).

(31)

Due to (i) and by setting N > m0 in (ii), G is analytic on C+ , continuous on C+ and vanishes when |s| → ∞. Let w = F1 (s) = 1/(s  + ). F1 is a bijective conformal mapping from C+ on bounded open domain  limited by Jordan curve , where  is a union of two circular arcs with end points 0 and 1/ , as shown in Fig. 1. Define F(z) = G ◦ F1−1 (z),

z ∈ .

(32)

Thus, F is analytic on  and continuous on its closure . All the assumptions in Mergelyan’s theorem (Rudin, 1987, Theorem are satisfied and there is a polynomial P (z) = N 20.5) n n=0 an z such that |F(z) − P (z)|  ,

∀z ∈ 

(33)

which, replacing z by F1 (s), yields   N     n an (F1 (s))   , ∀s ∈ C+ . G(s) −   n=0

Replacing G(s) by its definition (31), noting that Fn (s) = (F1 (s))n , and integrating the square of the obtained expression

1644

M. Aoun et al. / Automatica 43 (2007) 1640 – 1648

along the imaginary axis, yields the desired result: 

 2  j∞ N      2 an Fn+m0 (s) ds  |Fm0 (s)|2 ds, H (s) −  −j∞  −j∞

G2 (s) =

−5.19 1.73 + , s + 1.5 (s + 1.5)2

G3 (s) =

−10.39 1.73 15.59 + + . s + 1.5 (s + 1.5)2 (s + 1.5)3

j∞

n=0

  N      an Fn+m0 (s) ε. H (s) −   n=0

2

Thus, the linear span of {Fm }m0 ,m0 +1,...,∞ is dense in H2 (C+ ), which completes the proof.  +

Consequently, the linear span {Gm }m  m0 is dense in H2 (C ) too, and the orthonormal functions {Gm }m  m0 can be used to approximate any finite energy system. Important remarks: • The new basis has two tuning parameters: an eigenvalue and a commensurable order . • The only necessary conditions to build the basis are > 0 and 0 <  < 2. • When the commensurable order  equals 1, the built basis corresponds exactly to the classical Laguerre basis. Hence, the fractional Laguerre basis can be considered as an extrapolation of the rational Laguerre basis for any commensurable order  ∈]0, 2[. • For a large number of functions (large M), bad conditioning of F, FT matrix (to be inverted in (24)) may occur. Orthogonal functions can then be computed using a recursive Gram–Schmidt procedure: ⎧ Fm0 (s) ⎪ ⎪ ⎪ Gm0 (s) = Fm (s) , ⎪ 0 ⎨  Fm (s) − m−1 (34) i=m0 Fm (s), Gi (s) Gi (s) , Gm (s) = ⎪ m−1 ⎪ ⎪ Fm (s) − i=m0 Fm (s), Gi (s) Gi (s) ⎪ ⎩ for m > m0 . 3.5. Examples When the differentiation order  = 0.4 and the eigenvalue = 1.5, condition (29) imposes starting from m0 = 2. Hence, the first three orthonormal functions are G2 (s) = G3 (s) = G4 (s) =

3.18 (s 0.4 + 1.5) 4.62

(s 0.4 + 1.5)2 5.87 (s 0.4

+ 1.5)

2

+ +

1.73 , s + 1.5

G1 (s) = G2 (s) = G3 (s) =

1.49 , + 1.5

s 1.5

−2.38 1.06 + , 1.5 + 1.5 (s + 1.5)2

s 1.5

−2.77 1.06 3.56 + + . s 1.5 + 1.5 (s 1.5 + 1.5)2 (s 1.5 + 1.5)3

(37)

Fig. 2 shows impulse responses of the first three orthogonal functions for = 1.5 and various values of . One can notice that fractional orthonormal functions are more oscillatory when the differentiation order increases beyond 1. Also the impulse response at t = 0 can be finite or infinite depending on . Proof is given by the initial value theorem:  ∞ if 0 <  < 1, √ lim gm (t) = lim sGm (s) = (38) 2 if  = 1, s→∞ t→0 0 if1 <  < 2. 4. System identification using fractional Laguerre basis Fractional Laguerre basis can be used in output error identification with fixed denominator models (Wahlberg, 1991). Prior knowledge can be used to fix both tuning parameters  and . These two parameters are then plugged in (26), and the orthonormal basis is synthesized using (25). Then, Fourier coefficients are computed using a least-squares estimate. Assume u(t), y(t), t ∈ [0, T ] input and output data generated using a finite energy linear fractional model H. So H (s) can be approximated by a linear combination of orthonormal functions with: M 

gm Gm (s) = gT G(s),

(39)

m=m0

−18.48 (s 0.4 + 1.5)

where G is defined by (20) and , 3

−53.05 (s 0.4

+ 1.5)

3

+

g = [gm0 gm0 +1 . . . gM ]T . 95.50 (s 0.4

+ 1.5)

4

.

(35)

When the differentiation order is integer,  = 1, and the eigenvalue = 1.5, the first three orthonormal functions are G1 (s) =

This corresponds exactly to the classical definition of Laguerre functions (17), which are thus considered as a special case of the new fractional basis when  = 1. When the differentiation order  = 1.5 and the eigenvalue = 1.5, the first three orthonormal functions are

H (s) ≈

, 2

(36)

The truncation order M is fixed to obtain a satisfactory approximation. Akaike and Young information criteria can be used. The identification procedure consists of computing optimal coefficient vector g = [gm0 , gm0 +1 , . . . , gM ]T which minimizes the least-squares error:  1 T J= (ε(t))2 dt = ε, ε , (40) T 0

M. Aoun et al. / Automatica 43 (2007) 1640 – 1648

1645

1 2 G3 G2 G1

Input

1.5

0.5 0 −0.5

1

−1

0.5

0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

2

0

1 Output

−0.5 −1 0

5

10

15

0 −1 −2

Time(s)

Fig. 3. Input and output identification data.

2 G3 G2 G1

1.5

Setting

1

uG (t) = [uGm0 (t) uGm0 +1 (t) · · · uGM (t)], the optimum estimation of Fourier coefficients gˆ is given by the least-squares formula: −1  T  T uG (t)T uG (t) dt uG (t)T y(t) dt. (42) gˆ =

0.5 0 −0.5

0

−1 0

5

10

15

Time(s)

0

Or, after a numerical discretization, by defining Y as a column vector of system’s output and  as a regression matrix where columns are filter outputs, gˆ can be approximated by gˆ = (T )−1 T Y.

1

All properties of least-squares estimates (persistant excitation, variance on estimates) as stated in Ljung (1999) apply in this context.

G3 G2 G1

0.8 0.6

(43)

0.4 0.2

4.1. Example

0

To illustrate the use of fractional Laguerre bases in system identification, the following academic system is simulated and then identified in a noisy context:

−0.2 −0.4 −0.6 −0.8 0

5

10

15

Time(s) Fig. 2. Impulse responses of the first three orthogonal functions for = 1.5 and for different values of .

where ε(t) = y(t) −

M 

gm uGm (t),

(41)

m=m0

where y(t) and uGm (t) are, respectively, the system and orthogonal network outputs: uGm (t) = Gm (t) ⊗ u(t).

H (s) =

1 1 1 + + . s 0.7 + 2 s 0.8 + 2 s 0.9 + 1

(44)

The pseudorandom binary sequence used as input signal and the output signal are plotted in Fig. 3. Output signal is corrupted by a stationary zero mean Gaussian white noise with a signalenergy to-noise ratio (SNR) arbitrarily set to 10 log10 ( Signal Noise energy ) = 13 dB. H is approximated using two fractional Laguerre functions. Least-squares error, J, is computed in terms of differentiation order  varying from 0.55 to 1.9 and eigenvalue varying from 0.55 to 3. For each pair of (, ), generating functions (26) are orthogonalized, and then Fourier coefficients, g, computed using the least-squares formula (42). Isocontours of J are plotted in Fig. 4 for  ∈ [0.55, 1] and ∈ [0.5, 3]. Then, tuning parameters are chosen around the optimal values (, )=(0.75, 1.75).

1646

M. Aoun et al. / Automatica 43 (2007) 1640 – 1648

1

3

True model Fractional Laguerre Classical Laguerre Classical GOB

2 −12.75

−12

2.5 −12 .75

.75

0

2.

−12.93

5

−11

93

−12 .7

5

−9

−12

−12

2.7

−10

−1

−1

1.5

0.5

3 −12

Eigenvalue

−12 .9

2

−10

−11

−1

1 0 −1

−12

1

0.5 0.55

−12 .75 −12.93

−1

18 1

−1 12

−

0.6

0.65

0.7

−0.5

0.75 0.8 0.85 Fractional order

0 −1 −9

0.9

−8

−7

0.95

1

18.5

19

19.5

20

20.5 21 Time(s)

21.5

22

22.5

23

Fig. 5. Validation output data—the difference between the true system, the fractional Laguerre model, the Laguerre model and the generalized orthogonal basis model is almost indistinguishable.

Fig. 4. Isocontour of the least-squares error in dB.

The optimal differentiation order is far from an integer value, as expected, because the initial system (44) is fractional. The optimal model is Hˆ (s) = 1.90G1 (s) − 0.29G2 (s),

(45)

Magnitude (dB)

Bode

G1 (s) = G2 (s) =

1.51 , + 1.75

s 0.75

1.34s 0.75 − 4.68 . s 1.50 + 3.50s 0.75 + 3.06

The normalized residual,  T 2 0 (ε(t)) dt ≈ −13 dB, JdB = 10 log  T 2 0 (y(t)) dt

(46) (47)

(48)

is mainly due to the injected noise (as SNR is 13 dB). The modelling error is hence very small, which is confirmed when noise-free validation data are applied to the true system and to the obtained model. For further comparison, the system is identified using classical orthogonal bases. Six Laguerre functions or two GOB functions with optimal poles provide similar approximation errors: JdB = −12.5 dB. Optimal Kautz poles converge to the optimal Laguerre pole (the imaginary parts equal zero). Fig. 5 shows that time-domain outputs of the true system (44) and the identified models are almost identical. A better distinction is obtained in the frequency domain and Fig. 6 shows that the fractional Laguerre model clearly exhibits a better approximation. Fractional models can indeed have any asymptotic slope in the gain diagram and any asymptotic phase lock in the phase diagram (Oustaloup, 1995). 5. Conclusion An orthonormal fractional Laguerre basis has been synthesized for system approximations. It has two tuning parame-

Phase (deg)

where G1 (s) and G2 (s) are the orthonormal functions:

10 0 −10 −20 −30 10−1

0 −20 −40 −60 −80 −100 10−1

100 101 Frequency rad/s

102

100 101 Frequency rad/s

102

Fig. 6. Bode diagrams of the true system, the fractional Laguerre model, the Laguerre model, and the GOB model.

ters: an eigenvalue and a commensurable differentiation order  ∈]0, 2[. Classical Laguerre functions are shown as a special case of the new basis corresponding to  = 1. This is the first fractional basis ever developed. An academic system identification example shows the advantage of using a fractional basis as compared to a rational one. Appendix A. Computing scalar product Fh , Fm  Fh and Fm are the orthogonalized functions defined by (26)   1 1 Fh , Fm = , (A.1) , (s  + )h (s  + )m where ∈ R+∗ ,  ∈ R+∗ , h m0 , m m0 and m0 =1/2+1:    +∞ 1 1 1 Fh , Fm = d (A.2) 2 −∞ ((j ) + )h ((j ) + )m

M. Aoun et al. / Automatica 43 (2007) 1640 – 1648

then

1647

The integral (A.6) is

Fh , Fm =

1 2



∞

0

1 + 2



d

((j ) ∞

0

+ )h (j 

+ ) d

m

((−j ) + )h ((−j)  + )m

. (A.3)

The following change of variable is applied:  = x, d = (1/)x 1/−1 dx. Also, by defining:   ∞  1 x 1/−1 dx I , , , h, m = . (A.4)  ( x + 1)h (x + 1)m 0 Eq. (A.3) can be rewritten as    1 −1 j 2  −1 −j 2  1 Fh , Fm = I e , e , , h, m  2 h+m    1 +I ( −1 e−j 2  , −1 ej 2  , , h, m) . (A.5)  Let q = 1/ and  be the non-integer part of 1/. Depending on the value of , integral I is computed differently. Hence the following two decompositions:

  1 I , , , h, m    h−1   −1 =  Al (−1)h−l−1 h−l−1

sin() l=0   m−1   −1 +  (A.10) . Bl (−1)m−l−1 m−l−1  sin() l=0

• If 0 <  < 2 and 1/ ∈ N   1 I , , , h, m   ∞# h−2   Ah−1 + Bm−1 Al = + ( x + 1)(x + 1) ( x + 1)h−l 0 l=0 $ m−2  Bl dx, (A.11) + (x + 1)m−l l=0 where ⎛

• If 0 <  < 2 and 1/ ∈ /N   ∞  h−1 1 Al x −1 I , , , h, m =  ( x + 1)h−l 0 l=0  m−1  Bl x −1 dx, + (x + 1)m−l l=0 

Al

=

(m − 1)! ⎛

l

(m − 1)! l

k=Max(0,l−q+1)

k=Max(0,l−q)

⎞

Bl =

(q − 1)! (h − 1)!l

⎛

l 

×

k=Max(0,l−q+1)

⎜ (h + k − 1)! ⎜ ⎜ ⎝ k!(l − k)!(q + k − l − 1)! ⎞

and q! Bl = (h − 1)!l ⎛

k=Max(0,l−q)

⎞

By applying formula [3.194, 4, p. 285] of Gradshteyn and Ryshik (1980), ∞

x −1 dx

(−1)l−1 = 

sin() ( x + 1) l

(− )k (−)l−q−k+1 ⎟ ⎟ ×  h+k ⎟ . ⎠

− +1 

l 

⎜ (h + k − 1)! (− )k (−)l−q−k ⎟ ⎟ ⎜ ×⎜  h+k ⎟ (. A.8) ⎠ ⎝ k!(l − k)!(q + k − l)!

− +1 

0

(A.12)

and

⎜ (m + k − 1)! (−)k (− )l−q−k ⎟ ⎟ ⎜ ×⎜  m+k ⎟ (A.7) ⎠ ⎝ k!(l − k)!(q + k − l)!  − +1



⎜ (m+k−1)! ⎜ ⎜ ⎝ k!(l−k)!(q+k−l−1)!

(−)k (− )l−q−k+1 ⎟ ⎟ ×  m+k ⎟ ⎠  − +1

(A.6)

l 

q!

l 

⎞

where Al =

(q − 1)!



−1 l−1

 if l > .

(A.13)

Integral I can be computed using: 

∞ 0

and 

(A.9) 0

dx ( x + 1)l

∞

=

1

(l − 1)

if l > 1

dx ln( ) − ln() = . ( x + 1)(x + 1)

−

(A.14)

(A.15)

1648

M. Aoun et al. / Automatica 43 (2007) 1640 – 1648

Integral (A.11) is now:   1 I , , , h, m   (Ah−1 + Bm−1 ) = (ln( ) − ln())

− h−2 m−2   Al Bl + + .

(h − l − 1) (m − l − 1) l=0

(A.16)

l=0

The scalar product (A.1) can be deduced from (A.5) and (A.10) when 1/ ∈ / N, and from (A.5) and (A.16) when 1/ ∈ N. References Abbott, P. C. (2000). Generalized Laguerre polynomials and quantum mechanics. Journal of Physics A: Mathematical and General, 33(42), 7659–7660. Akçay, H., & Ninness, B. (1999). Orthonormal basis functions for modelling continuous-time systems. Signal Processing, 77, 261–274. Aoun, M. (2005). Systèmes linéaires non entiers et identification par bases orthogonales non entières. Ph.D. Thesis, Université Bordeaux 1, France. Aoun, M., Malti, R., Levron, F., & Oustaloup, A. (2004). Numerical simulations of fractional systems: An overview of existing methods and improvements. International Journal of Nonlinear Dynamics and Chaos in Engineering Systems, Special Issue on ‘Fractional Derivatives and Their Applications’, 38(1–4), 117–131. Battaglia, J.-L., Le Lay, L., Batsale, J.-C., Oustaloup, A., & Cois, O. (2000). Heat flux estimation through inverted non integer identification models. International Journal of Thermal Sciences, 39(3), 374–389. Darling, R., & Newman, J. (1997). On the short behavior of porous intercalation electrodes. Journal of the Electrochemical Society, 144(9), 3057–3063. El Sayed, A.-M. (1999). On the generalized Laguerre polynomials of arbitrary (fractional) orders and quantum mechanics. Journal of Physics A: Mathematical and General, 32, 8647–8654. Garnett, J. B. (2007). Bounded analytic functions. (Revised 1st ed.), Springer: Berlin. Gradshteyn, I.-S., & Ryshik, I.-M. (1980). Table of integrals, series, and products. New York: Academic Press. Ljung, L. (1999). System identification—theory for the user. (2nd ed), Englewood Cliffs, NJ: Prentice-Hall. Malti, R., Aoun, M., Levron, F., & Oustaloup, A. (2003). H2 norm of fractional differential systems. In Proceedings of design engineering technical conferences and computer and information in engineering conferences, Chicago, Illinois, USA, September 2003. New York: ASME. Malti, R., Ekongolo, S.-B., & Ragot, J. (1998). Dynamic SISO and MIMO system approximations based on optimal Laguerre models. IEEE Transactions on Automatic Control, 43(9), Matignon, D. (1998). Stability properties for generalized fractional differential systems. ESAIM Proceedings—Systèmes Différentiels Fractionnaires—Modèles, Méthodes et Applications, 5, 145–158. Ninness, B., & Gustafsson, F. (1997). A unifying construction of orthonormal bases for system identification. IEEE Transactions on Automatic Control, 42(4), 515–521. Oldham, K.-B., & Spanier, J. (1974). The fractional calculus: Theory and applications of differentiation and integration to arbitrary order. New York, London: Academic Press. Ortigueira, M.D., & Tenreiro Machado, J.A. (Eds.) (2006). Signal processing—special issue: Fractional calculus applications in signals and systems, (Vol. 86). Amsterdam: Elsevier, October 2006.

Oustaloup, A. (1995). La dérivation non entière: théorie, synthèse et applications. Paris: Hermès. Poinot, T., & Trigeassou, J.-C. (2003). A method for modelling and simulation of fractional systems. Signal Processing, 83, 2319–2333. Rudin, W. (1987). Real and complex analysis. (3rd ed), New York: McGrawHill. The American Heritage䉷(Ed.) (2000). Dictionary of the English language. (4th ed.) Houghton Mifflin Company. Van den Hof, P.-M.-J., Heuberger, P.-S.-C., & Bokor, J. (1995). System identification with generalized orthonormal basis functions. Automatica, 31(12), 1821–1834. Wahlberg, B. (1991). System identification using Laguerre models. IEEE Transactions on Automatic Control, 36, 551–562. Mohamed Aoun was born in Tunisia in 1975. He received his Ph.D in automatic control in 2005 from the University of Bordeaux, France. He is currently an Associate Professor in automatic control, electrical engineering, and computer engineering at ENIG, Tunisia and a member of its MACS research laboratory. His research interests include automatic control, system identification, and fractional differentiation. Rachid Malti was born in Serbia in 1972. He received his electrical engineering diploma from INELEC, Boumerdès, Algeria in 1994, his M.Sc. and Ph.D. in 1996 and 1999 in automatic control from INPL, Nancy, France. He held the position of an Associate Professor in automatic control, electrical engineering, and computer engineering at the University of Paris from 1999 to 2004. In 2004, he moved to the University of Bordeaux where he joined the IMS-CRONE team to work in the field of fractional differentiation and its applications in automatic control and system identification. François Levron was born in France in 1944. He received his Ph.D. in Mathematics in 1969 from the University of Bordeaux. He is currently an Associate Professor in the Department of Mathematics of the University of Bordeaux. His research interests include numerical analysis, fractional calculus, and mathematics applied to physics and engineering.

Alain Oustaloup was born in France in 1950. He received his engineering diploma from ENSEIRB in 1973, and his Doctor of Engineering and Doctor of Science in 1975 and 1981 from the University of Bordeaux. He is currently a Professor and the head of the Automatic Control department at ENSEIRB. He set up and is the head of the IMS-CRONE team at the University of Bordeaux. He is currently a Research Evaluator at the French Ministry of Education and Research. His research interests include fractional differentiation, its synthesis, and its applications in engineering sciences, particularly in automatic control. He conceived CRONE control (French acronym for “Commande Robuste d’Ordre Non Entier”). He is the author and co-author of five books including Oustaloup (1995). He was awarded the Afcet Trophy in 1995 and the CNRS Silver Medal in 1997.