- Email: [email protected]
Diffusive Realisations of Fractional Integrodifferential Operators: Structural Analysis Under Approximation☆
Copyright 0 IFAC System Structure and Control, Nantes, France, 1998
DIFFUSIVE REALISATIONS OF FRACTIONAL INTEGRODIFFERENTIAL OPERATORS: STRUCTURAL ANALYSIS UNDER APPROXIMATION * David HELESCHEWITZ * Denis MATIGNON **
* ENST, Signal and Image Processing Dept. 46 Rue Barrault F-75634 Paris Cedex 13, FRANCE. e-mail: [email protected] ** ENST, Signal and Image Processing Dept. 8 CNRS, URA 820 e-mail: [email protected]
Abstract: Fractional integrals and derivatives are usually presented as non-standard convolutions, but they are also linked to classical diffusive operators on a Hilbert state space of infinite dimension. A variety of equivalent realisations are proposed, which can naturally be approximated by truncation in finite dimension. This process gives rise to convergent approximation algorithms from which some numerical simulations are derived. Copyright © 1998 IFAC Resume : La derivation et l'integration fractionnaires, habituellement presentees comme convolutions non standards, sont intrinsequement liees a des operateurs diffusifs classiques dans un espace d'etat hilbertien de dimension infinie. Diverses realisations equivalentes sont proposees, elles sont ensuite naturellement approchees par troncature en dimension finie , donnant lieu a des algorithmes d'approximation convergents a partir desquels des simulations numeriques sont effectuees. Keywords: Fractional Differential Equations, Realisation theory, Rational approximations, Numerical simulation
1. INTRODUCTION
provide a powerful equivalent description which happens to be classical in an infinite-dimensional framework and allows for convergent approximation.
Fractional differential systems have become more and more popular in control theory both for modelling (see e.g. [MAM98a]) and controller design purposes (see [MbM95,MAM97,MAM98b]). Such systems correspond to non-standard convolutions which are difficult to handle from a theoretical point of view (see e.g. [Mat96])j moreover numerical simulation is not straightforward. The diffusive realisations of fractional integrals and derivatives, originally introduced in [MAMb93j ,
The paper is organized as follows: some equivalent diffusive realisations are derived in section 2 and finite-dimensional approximations are presented in section 3.
2. DIFFUSIVE REALISATIONS AND CONSEQUENCES
*
The authors would like to tha.nk G. Montseny (LAASjCNRS - Toulouse) a.nd J. Audounet (MIP jCNRS - Toulouse) for fruitful discussions a.nd useful advices.
The goal of this section is to present diffusive realisations of fractional integrals. Analogous re227
with J..1.o ~ (sin(7l'0)j7l') €-o . Hence, using Fubini's theorem, fractional integration of u E L2(0, T; JR) now reads:
alisations also exist for fractional derivatives (see [MAM98b]) but will not be presented in this paper for simplicity sake. In subsection 2.1, the very infinite-dimensional nature of fractional integrals is revealed by rewriting the input-output relation in a more explicit way; which leads to a straightforward realisation on a convenient Hilbert space. For finite-dimensional systems, realisation of a given transfer function is highly not unique, and even more so for infinite-dimensional systems; therefore in subsection 2.2, a wide variety of equivalent realisations is being presented and analysed, which naturally encompasses the so-called original diffusive realisations introduced by Montseny and Audounet (see [MAMb93,MbM95]) .
+00
y(t)
OtX(€, t) {
(3)
= -€ X(€, t) + J..1.o(€) u(t)
y(t)=JX(€, t)~
(4)
o where € E JR+ 'aqd the state X(., t) belongs to an appropriate Hilbert space Vo: , with initial condition X(€ , t = 0) :::: O. Remark 1. The state space is necessarily of infinite dimension, an intrinsic property which is somehow hidden in the usual input-output convolution relation (1) . It must be noted that (4) is expressed with first order derivatives only, which recasts fractional integral operators in a more tractable form and in a commonly used framework (such as semigroup theory in functional analysis) ; moreover it allows for numerical approximation schemes that are also highly standard techniques in numerical analysis (see section 3).
t
(t - r)O-l u (r) dr
* u)(t) ~
+00
do)
J
(ee
which can happily be realised in a straightforward way, such as:
Using classical notations in control theory (u for input and y for output), recall first that fractional integration of order 0 , y = F"U, generalises cascaded integration to an arbitrary real order. Therefore, it takes the form of a convolution with a causal power law kernel Yo ~ t~-1 where r denotes the Euler gamma function , namely:
= r(~)
JJ..1.o(~) o
2.1 Infinite-dimensional realisation of fractional integrals
y(t)
=
(1)
o
=(Yo *u) (t)
2.2 Equivalent realisations
(2)
The Laplace transform of the kernel Yo is defined in the right-half complex plane as: .c[Yo ]
= s-o
Re(s)
System (4) is a natural realisation of (3) and can be expressed in a more general form:
>0
OtX = AX { y=CX
It must be noted that the impulse response Yo: does not belong to the class of polynomialexponentials that arises in classical Ordinary Differential Equations, when linear and timeinvariant. Nevertheless, in the case 0 E ]0,1[, Yo: can be shown to be a continuous superposition of purely damped causal exponentials edt) = e-~tlt~o with weight J..1.o: . Indeed, since the real axis JR+* belongs to the convergence domain of C [Yo:] (s) , it is possible to rewrite Yo:(t) ex t-(I-o) in the following form:
A : AX = -€X being fixed, it follows that the measure J..1.o is dispatched between 5 and C, with only constraint b(€)c(~) = J..1.o:(€). For instance: o
- r(o)r(1 - 0)
JCO: -~t~
+00
=
o
=
• b J..1.o: and c 1 is realisation (4). • b = 1 and c = J..1.o: is the dual realisation of (4) . Note that the state X(~, t) does not depend upon 0, which means that fractional integrals of different orders of a same input u can be obtained from only one infinitedimensional system with several outputs, i.e. '101,02,'" E]o,I[:
+00
1
(5)
where A is a dissipative infinite-dimensional operator (dissipativity is a remarkable property of fractional integrodifferential operators that is preserved by diffusive realisations), 5, the control operator, is a scalar multiplier b(€) and C, the observation operator, is a continuous linear funcoo tional CX = fo+ c(€)X(€,.)~. Therefore from (3), other realisations may be proposed:
=
-
+ 5u(t)
e
J
J..1.o:(Oe - et ~
o 228
= -~X(~,t) +u(t)
8tX(~,t) +00
Yl =
f f
Jlo 1
(~) X(~,.) cl{ =
carefully presented in subsection 3.5, with special attention paid to the zeros of H.
IOl u
(6)
o
+00
Y2 =
Jlo 2 (~) X(~,.) cl{ = I
0
2
3.1 State-space realisation
u
Finite-dimensional approximations of the diffusive realisation (4) are obtained over a finite network of n points, say N n = {~dk=l,n C JR+, namely:
o
= c = ~ is the balanced realisation which enables to prove the positivity property of the operator. Such a property is crucial, for it implies dissipativity (see e.g. [MAM98a]).
• b
-~kXk(t) + Jlo(~k)U(t)
!Xk(t) = { Xk(t) ~ X (~k' t)
e
= -47r2 X where ~ E IR leads to the original diffusive realisation introduced by Montseny and Audounet in [MAMb93,MbM95], which can be easily interpreted in a physical framework. Indeed, it corresponds to the Fouriertransform (with respect to a space variable r) of a so-called diffusion Partial Differential Equation, namely: o A : AX
8tZ(r, t) = 8;Z(r, t) { y(t) = Z(O, t)
+ zio(r) u(t)
(8)
for k = 1, ... n, with initial condition Xk(O) = o. The approximate output y(t) is defined in two steps: • A reconstruction of the state by means of interpolating functions A k , i.e. n
X(~, t)
=L
Xk(t)Ak(~)
(9)
with {Ad chosen such that X(~k' t) • Then,
= Xk(t).
(7)
k=l
where r E lR, Z( ., t) is the state with initial condition Z (r, t = 0) = 0, \:Ir E IR and zio(r) is the inverse Fourier transform of the measure l 20 lIo (O = 2 sin(cm)(27rIW .
y(t)
f
= X(~, t) ~
(10)
R+
Thus, the following n-dimensional state-space realisation arises:
o Conversely, it can be chosen to fix B and C first and derive A afterwards. The balanced measurefree realisation, which corresponds to b = c = 1, leads to A : AX = _~{3 X with f3 = (1 _ a)-l .
!X(t) = AX(t) { y(t) = CX(t)
This realisation is very convenient from a theoretical point of view (especially when proving convergence of the resulting approximation, see subsection 3.2). Note that the measure does not appear anymore in this realisation, which is therefore called measure-free.
+ Bu(t)
(11)
with • • • •
=
=
X(t) (Xl(t),·· · ,xn(t)f, X(O) 0, A = diag (-6, ... , -~n), B = (bb···, bn)T with bk = Jlcr(~k)' C = (Cl,··· ,Cn ) with Ck = fR+ Ak(~) ~
The corresponding transfer function is: 3. FINITE-DIMENSIONAL APPROXIMATION
H(s)
=
t
k=l
This section is concerned with finite-dimensional approximations of the diffusive realisations seen previously. Choice has been made to consider realisation (4) only and subsection 3.1 shows how this realisation can give rise to a rational approximation. However analogous approximations can also be derived for all equivalent realisations.
S
bkCk , + ~k
&(s) > 0
(12)
Remark 2. There is no a-priori constraint on the choice of the poles {~kh=l,n . This important property of the approximation (namely freedom of choice) will be discussed in subsection 3.3.
3.2 Convergence of the approximation
In subsection 3.2, the property of convergence of the approximation is presented and the main
Convergence of the rational approximation presented above is quite tricky to prove from an analytic point of view. The interested reader is therefore referred to [MAM98b] for further details. However, some elements of the proof need to be outlined in order to justify the approximation process presented later on in subsection 3.3.
steps of the proof are being described. Then in subsection 3.3, the approximation process is detailed by a step by step procedure, whereas subsection 3.4 deals with a constructive choice of the poles of the transfer function of approximation H. To end with, some numerical simulations are 229
Theorem 3. By a convenient choice of the interpolating functions Ak (see subsection 3.3), for a given U E L2(0, T; lR),
IIX -
-70
(13)
yllL2(O,T;R) -70
(14)
XIIL2(o,T;V",)
Ily -
step 1: range of approximation
nm = [~~in' ~~a:t]' step 2: choice the poles Any choices of the poles seem possible in nm and in fact they are, as long as they meet the convergence condition (18) . Poles can be chosen as follows:
where Vc> is an appropriate Hilbert space for X andX .
• an arithmetic sequence: ~k+1 = ~k + ra(n) • a polynomial sequence: ~k = krp(n) 6 • a geometric sequence: ~k+1 = rg(n) ~k
·
In order to define the finite-dimensional approximation in subsection 3.1, a network of points N n has been introduced. Convergence of the approximation scheme is based upon a natural filling property of a network which should cover the whole frequency domain as the order of approximation m increases. Let
N;:'
= {O < ~:in = ~;:n < ~~ < ~k < ~~I <~: = ~:a:t < +oo}
...
The geometric choice has been widely used in the literature since it produces a linear sequence of pulsations on a Bode plot. In practise, it realises a good compromise, between the width of the band of interest and the 'Gomplexity of the resulting approximation. Such a choice has been made, for instance, by Alain Oustaloup and co-authors (see e.g. [Ous83]) when fitting phase and magnitude of fractional integrals or derivatives over a bounded range of frequencies. But, it happens that this choice is not unique and many others are to be investigated. Examples and simulations are given in subsection 3.5.
(15)
Then, the proof of theorem 3 can be sketched as follows: • decompose IIX - Xliv", into three parts in the frequency domain, • choose m such as to bound IIX - Xliv", over the intervals ]0, ~:in] and [~~a:t' +oo[ by using a Lebesgue-type argument, • choose n such as to bound IIX - Xliv", over the fixed compact interval nm = [~~in'~~a:t] by using the convergence property of the interpolation process, • finally deduce convergence of the output y towards y in L2(0, T ; JR).
step 3: choice of the interpolating functions Finite elements interpolation corresponds to a standard method in numerical analysis (see [Sain96, chapter 3]); the better the quality of the approximation (convergence, speed, regularity) , the greater the order of the finite element. For simplicity sake, finite elements of order 1 have been chosen:
Hence, natural convergence conditions write out: 'om.n
cm .
m-.+oo -7
0
(16)
cm
m-.+oo -7
+00
(17)
'omax
Vm EN sup I~~I - ~kl n~oo 0
Note that the interpolating function Ak is null outside the interval [~k-I' ~k+l] and Ak(~k) = 1.
(18)
I~k~n
step 4: continuation of interpolation Clearly, the quality of the approximation depends on the quality of the reconstruction of the state X by means of the interpolating functions Ak. Since this reconstruction has to be performed over the whole frequency domain, it follows that the interpolating functions Al and An must have a particular shape, namely:
3.3 Approximation Process
Suppose that the band of interest is of the form [Wmin, wmax ] with Wmin > O. The question is then to build the network of poles N;: to realise a convergent approximation over this bounded range of frequencies. The construction of N;: is entirely based upon the convergence result and can be seen as a recursive process which stops as soon as the resulting approximation meets the desired accuracy. step 0: initialization m = O. Fix ~in = Wmin and ~a:t
Indeed, let X(~, t) be the state obtained for a given input u and denote X(~, 5) its Laplacetransform in 'Re(s) > O. Using (4), it comes:
= w max . 230
(21)
which implies:
X(~, t) (~o P,Q(~) Jot u(r) dr
!
X(~, t) (~~oo P,Qi~)
i
i t
(22)
~
r
iI
u(t)
l
Fig. 1. Bode plot of iI with a geometric choice of poles. Number of poles n = 20, [Wmin , Wmax ] = [10 2 ,10 5 ] . There is a good compromise between the size n of the system and the bandwidth.
step 5: accuracy criterion Either the approximation provides the desired accuracy or not. In the latter case, one should first increase n to coarsen the grid in om . Go back to step 2. IT accuracy is not improved then change ~:in and ~:ax according to conditions (16) and (17). For instance ~mtn cm+l = 110 'm cm.In and ":.maxcm+ 1 10 ~:ax' Then go back to step 1. Convergence of the approximation ensures that one will get the desired accuracy after a finite number of trials.
3.5 Numerical simulations
This subsection is devoted to numerical simulations of fractional integration of order Q = 0.2 over a bounded range of frequencies [Wmi n, w max ]. Two different choices of poles location have been examined, geometric and polynomial. For each cho~ce , there is one figure showing the Bode plot of H and another one describing its zeros location.
3.4 Constructive choice of the poles
Bode plots Figures (1) and (2) show Bode plots of iI for a geometric and a polynomial choice of poles respectively. They both fit the behaviour of the fractional integrator over a bounded range of frequencies. For the polynomial choice, there are not enough poles in the range of frequencies (n too small) , which gives rise to oscillations on the phase plot. Conversely, for the geometric choice, n is great enough to avoid oscillations. <>
A constructive choice for the network of poles
N:;' is presented here in the case of a geometric sequence. For a fixed order of approximation m, condition (18) writes out:
11
sup I~fl n~oo 0
·20
- 60
is completely determined.
Ir(n) -
-60
~
Recalling definition (9) of X(~, t) yields the desired result. Note that the final expression of A1 and An is determined by continuity in 6 and ~n respectively. At this stage,
-
(23)
lSk~n
with r(n) = (~:ax/~:in)1/(n-1) . For a given m , (23) is straightforward. Now, as m increases, om+1 ::> om. Therefore, n(m + 1), the number of points in om+l , should be greater than n(m) in order to get at least the same accuracyj it yields the condition:
Zeros location Zeros of the transfer function are not chosen. They appear as a consequence of the location of the poles and of the choice of matrices B and C. From simulations, poles and zeros alternate on the real axis 1R- and are asymptotically distributed in a regular way depending exclusively on Q and not on the AkSj
<>
Hence, two choices are proposed: ~:in cc
m-It
~:ax cc ml2 n(m) cc m IZ + E
or
zf E [~k' ~k-rl]
~:in ccexp(-mrl)
~:ax cc m r2 n(m) cc mrl+r2+E
with
zf denoting the kth
k
= 1, .. . ,n -
1
(25)
zero. Location of zf in
[~k ' ~k+l] can be described by the following local
where Lt, r1, 12 , r2 denote the speed of convergence of ~:in and ~:ax towards 0 and +00 respectively and where to > 0 must be chosen as small as possible to keep n to a reasonable value. However, other choices can be investigated.
slopes on Bode plot:
Ri ght SIope
231
m _ log(zk) -log(~k) k log( ~k+l) - log( ~k )
(26)
0 .•
0.•
0.7
t.
i
~
i°.5 a:
o
0
0.'
0.3 Q.2
o ... .
....... ~ . . . . . ~ .-.
O"!-D--=----::-----:----:--~s--7----:~---:---'---..J,D
"'* '" poIK
FrtlQ\WlCy(l'IIdIMe)
Fig. 4. Right (0) and left (+) local slopes with a polynomial choi'ce of poles. Number of poles n = 10. '.
Fig. 2. Bode plot of if with a polynomial choice of poles. Number of poles n = 10, [Wmin, w max ] = [10 3 ,10 6 ]. The number of poles n is too small to describe the range of frequencies and a careful look at the graph shows little oscillations on the phase plot.
5. REFERENCES [MAM98b] G. Montseny, J. Audounet, and D. Matignon. Diffusive wave-absorbing feedback controls for the boundary stabilization of a flexible beam. submitted for publication. [MAMb93] G. Montseny, J . Audounet, and B. Mbodje. Optimal models of fractional integrators and application to systems with fading memory. In Conference IEEE Systems, Man and Cybernetics, Le Touquet, France, 1993. IEEE-SMC . [MAM97] G. Montseny, J . Audounet, and D. Matignon. Fractional integrodifferential boundary control of the Euler-Bernoulli beam. In Conference on Decision and Control, pages 4973-4978, San Diego, California, December 1997. IEEE-CSS , SIAM. [MAM98a] D. Matignon, J. Audounet, and G. Montseny. Energy decay for wave equations with damping of fractional order. In Fourth international conference on mathematical and numerical aspects of wave propagation phenomena, 3 pages, Golden, Colorado, June 1998. INRIA, SIAM. [Mat96] D. Matignon. Recent results in fractional differential systems theory. Technical Report 96 C 004, Ecole Nationale Superieure des Telecommunications, 1996. [MbM95] B. Mbodje and G. Montseny. Boundary fractional derivative control of the wave equation. IEEE Trans. Aut. Cont., 40(2):378-382, Feb. 1995. [Ous83] A. Oustaloup. Systemes asservis lineaires d'ordre /ractionnaire. Automatique. Masson, 1983. [Sain96] L. Sainsaulieu. Calcul scientifique. Masson, 1996.
0.' 0
0.35
1 0.3
I
~D.2S
If
o 0.2
+
~- - . • ... • -.• --. - .- - . -. -. - .• .. ~ .- . - . - 9 o · +
. . 0
o 0.15
0.10~---7---:---;-----7---:'':-D--:'::,,---',.':---,:'::-.--'-".....:o~20 indexafpatn
Fig. 3. Right (0) and left (+) local slopes with a geometric choice of poles. Number of poles n = 20. Left SIope
m _
k
-
10g(ZM-l) -log(ek+l) log(zk+l) -log(zk)
(27)
Figures (3) and (4) show the value of the right and left slopes in the geometric and polynomial case respectively. They are both located around Cl = 0.2 in the range of approximation whatever the choice of the poles.
4. CONCLUSION A powerful method of approximation of fractional integrals based upon diffusive realisations has been thoroughly examined; it proves to be very flexible and its extension to fractional derivatives and other pseudo-differential operators is shortly in prospect. Hence it will be applied to the simulation of Fractional Differential Equations. 232
DIFFUSIVE REALISATIONS OF FRACTIONAL INTEGRODIFFERENTIAL OPERATORS: STRUCTURAL ANALYSIS UNDER APPROXIMATION * David HELESCHEWITZ * Denis MATIGNON **
* ENST, Signal and Image Processing Dept. 46 Rue Barrault F-75634 Paris Cedex 13, FRANCE. e-mail: [email protected] ** ENST, Signal and Image Processing Dept. 8 CNRS, URA 820 e-mail: [email protected]
Abstract: Fractional integrals and derivatives are usually presented as non-standard convolutions, but they are also linked to classical diffusive operators on a Hilbert state space of infinite dimension. A variety of equivalent realisations are proposed, which can naturally be approximated by truncation in finite dimension. This process gives rise to convergent approximation algorithms from which some numerical simulations are derived. Copyright © 1998 IFAC Resume : La derivation et l'integration fractionnaires, habituellement presentees comme convolutions non standards, sont intrinsequement liees a des operateurs diffusifs classiques dans un espace d'etat hilbertien de dimension infinie. Diverses realisations equivalentes sont proposees, elles sont ensuite naturellement approchees par troncature en dimension finie , donnant lieu a des algorithmes d'approximation convergents a partir desquels des simulations numeriques sont effectuees. Keywords: Fractional Differential Equations, Realisation theory, Rational approximations, Numerical simulation
1. INTRODUCTION
provide a powerful equivalent description which happens to be classical in an infinite-dimensional framework and allows for convergent approximation.
Fractional differential systems have become more and more popular in control theory both for modelling (see e.g. [MAM98a]) and controller design purposes (see [MbM95,MAM97,MAM98b]). Such systems correspond to non-standard convolutions which are difficult to handle from a theoretical point of view (see e.g. [Mat96])j moreover numerical simulation is not straightforward. The diffusive realisations of fractional integrals and derivatives, originally introduced in [MAMb93j ,
The paper is organized as follows: some equivalent diffusive realisations are derived in section 2 and finite-dimensional approximations are presented in section 3.
2. DIFFUSIVE REALISATIONS AND CONSEQUENCES
*
The authors would like to tha.nk G. Montseny (LAASjCNRS - Toulouse) a.nd J. Audounet (MIP jCNRS - Toulouse) for fruitful discussions a.nd useful advices.
The goal of this section is to present diffusive realisations of fractional integrals. Analogous re227
with J..1.o ~ (sin(7l'0)j7l') €-o . Hence, using Fubini's theorem, fractional integration of u E L2(0, T; JR) now reads:
alisations also exist for fractional derivatives (see [MAM98b]) but will not be presented in this paper for simplicity sake. In subsection 2.1, the very infinite-dimensional nature of fractional integrals is revealed by rewriting the input-output relation in a more explicit way; which leads to a straightforward realisation on a convenient Hilbert space. For finite-dimensional systems, realisation of a given transfer function is highly not unique, and even more so for infinite-dimensional systems; therefore in subsection 2.2, a wide variety of equivalent realisations is being presented and analysed, which naturally encompasses the so-called original diffusive realisations introduced by Montseny and Audounet (see [MAMb93,MbM95]) .
+00
y(t)
OtX(€, t) {
(3)
= -€ X(€, t) + J..1.o(€) u(t)
y(t)=JX(€, t)~
(4)
o where € E JR+ 'aqd the state X(., t) belongs to an appropriate Hilbert space Vo: , with initial condition X(€ , t = 0) :::: O. Remark 1. The state space is necessarily of infinite dimension, an intrinsic property which is somehow hidden in the usual input-output convolution relation (1) . It must be noted that (4) is expressed with first order derivatives only, which recasts fractional integral operators in a more tractable form and in a commonly used framework (such as semigroup theory in functional analysis) ; moreover it allows for numerical approximation schemes that are also highly standard techniques in numerical analysis (see section 3).
t
(t - r)O-l u (r) dr
* u)(t) ~
+00
do)
J
(ee
which can happily be realised in a straightforward way, such as:
Using classical notations in control theory (u for input and y for output), recall first that fractional integration of order 0 , y = F"U, generalises cascaded integration to an arbitrary real order. Therefore, it takes the form of a convolution with a causal power law kernel Yo ~ t~-1 where r denotes the Euler gamma function , namely:
= r(~)
JJ..1.o(~) o
2.1 Infinite-dimensional realisation of fractional integrals
y(t)
=
(1)
o
=(Yo *u) (t)
2.2 Equivalent realisations
(2)
The Laplace transform of the kernel Yo is defined in the right-half complex plane as: .c[Yo ]
= s-o
Re(s)
System (4) is a natural realisation of (3) and can be expressed in a more general form:
>0
OtX = AX { y=CX
It must be noted that the impulse response Yo: does not belong to the class of polynomialexponentials that arises in classical Ordinary Differential Equations, when linear and timeinvariant. Nevertheless, in the case 0 E ]0,1[, Yo: can be shown to be a continuous superposition of purely damped causal exponentials edt) = e-~tlt~o with weight J..1.o: . Indeed, since the real axis JR+* belongs to the convergence domain of C [Yo:] (s) , it is possible to rewrite Yo:(t) ex t-(I-o) in the following form:
A : AX = -€X being fixed, it follows that the measure J..1.o is dispatched between 5 and C, with only constraint b(€)c(~) = J..1.o:(€). For instance: o
- r(o)r(1 - 0)
JCO: -~t~
+00
=
o
=
• b J..1.o: and c 1 is realisation (4). • b = 1 and c = J..1.o: is the dual realisation of (4) . Note that the state X(~, t) does not depend upon 0, which means that fractional integrals of different orders of a same input u can be obtained from only one infinitedimensional system with several outputs, i.e. '101,02,'" E]o,I[:
+00
1
(5)
where A is a dissipative infinite-dimensional operator (dissipativity is a remarkable property of fractional integrodifferential operators that is preserved by diffusive realisations), 5, the control operator, is a scalar multiplier b(€) and C, the observation operator, is a continuous linear funcoo tional CX = fo+ c(€)X(€,.)~. Therefore from (3), other realisations may be proposed:
=
-
+ 5u(t)
e
J
J..1.o:(Oe - et ~
o 228
= -~X(~,t) +u(t)
8tX(~,t) +00
Yl =
f f
Jlo 1
(~) X(~,.) cl{ =
carefully presented in subsection 3.5, with special attention paid to the zeros of H.
IOl u
(6)
o
+00
Y2 =
Jlo 2 (~) X(~,.) cl{ = I
0
2
3.1 State-space realisation
u
Finite-dimensional approximations of the diffusive realisation (4) are obtained over a finite network of n points, say N n = {~dk=l,n C JR+, namely:
o
= c = ~ is the balanced realisation which enables to prove the positivity property of the operator. Such a property is crucial, for it implies dissipativity (see e.g. [MAM98a]).
• b
-~kXk(t) + Jlo(~k)U(t)
!Xk(t) = { Xk(t) ~ X (~k' t)
e
= -47r2 X where ~ E IR leads to the original diffusive realisation introduced by Montseny and Audounet in [MAMb93,MbM95], which can be easily interpreted in a physical framework. Indeed, it corresponds to the Fouriertransform (with respect to a space variable r) of a so-called diffusion Partial Differential Equation, namely: o A : AX
8tZ(r, t) = 8;Z(r, t) { y(t) = Z(O, t)
+ zio(r) u(t)
(8)
for k = 1, ... n, with initial condition Xk(O) = o. The approximate output y(t) is defined in two steps: • A reconstruction of the state by means of interpolating functions A k , i.e. n
X(~, t)
=L
Xk(t)Ak(~)
(9)
with {Ad chosen such that X(~k' t) • Then,
= Xk(t).
(7)
k=l
where r E lR, Z( ., t) is the state with initial condition Z (r, t = 0) = 0, \:Ir E IR and zio(r) is the inverse Fourier transform of the measure l 20 lIo (O = 2 sin(cm)(27rIW .
y(t)
f
= X(~, t) ~
(10)
R+
Thus, the following n-dimensional state-space realisation arises:
o Conversely, it can be chosen to fix B and C first and derive A afterwards. The balanced measurefree realisation, which corresponds to b = c = 1, leads to A : AX = _~{3 X with f3 = (1 _ a)-l .
!X(t) = AX(t) { y(t) = CX(t)
This realisation is very convenient from a theoretical point of view (especially when proving convergence of the resulting approximation, see subsection 3.2). Note that the measure does not appear anymore in this realisation, which is therefore called measure-free.
+ Bu(t)
(11)
with • • • •
=
=
X(t) (Xl(t),·· · ,xn(t)f, X(O) 0, A = diag (-6, ... , -~n), B = (bb···, bn)T with bk = Jlcr(~k)' C = (Cl,··· ,Cn ) with Ck = fR+ Ak(~) ~
The corresponding transfer function is: 3. FINITE-DIMENSIONAL APPROXIMATION
H(s)
=
t
k=l
This section is concerned with finite-dimensional approximations of the diffusive realisations seen previously. Choice has been made to consider realisation (4) only and subsection 3.1 shows how this realisation can give rise to a rational approximation. However analogous approximations can also be derived for all equivalent realisations.
S
bkCk , + ~k
&(s) > 0
(12)
Remark 2. There is no a-priori constraint on the choice of the poles {~kh=l,n . This important property of the approximation (namely freedom of choice) will be discussed in subsection 3.3.
3.2 Convergence of the approximation
In subsection 3.2, the property of convergence of the approximation is presented and the main
Convergence of the rational approximation presented above is quite tricky to prove from an analytic point of view. The interested reader is therefore referred to [MAM98b] for further details. However, some elements of the proof need to be outlined in order to justify the approximation process presented later on in subsection 3.3.
steps of the proof are being described. Then in subsection 3.3, the approximation process is detailed by a step by step procedure, whereas subsection 3.4 deals with a constructive choice of the poles of the transfer function of approximation H. To end with, some numerical simulations are 229
Theorem 3. By a convenient choice of the interpolating functions Ak (see subsection 3.3), for a given U E L2(0, T; lR),
IIX -
-70
(13)
yllL2(O,T;R) -70
(14)
XIIL2(o,T;V",)
Ily -
step 1: range of approximation
nm = [~~in' ~~a:t]' step 2: choice the poles Any choices of the poles seem possible in nm and in fact they are, as long as they meet the convergence condition (18) . Poles can be chosen as follows:
where Vc> is an appropriate Hilbert space for X andX .
• an arithmetic sequence: ~k+1 = ~k + ra(n) • a polynomial sequence: ~k = krp(n) 6 • a geometric sequence: ~k+1 = rg(n) ~k
·
In order to define the finite-dimensional approximation in subsection 3.1, a network of points N n has been introduced. Convergence of the approximation scheme is based upon a natural filling property of a network which should cover the whole frequency domain as the order of approximation m increases. Let
N;:'
= {O < ~:in = ~;:n < ~~ < ~k < ~~I <~: = ~:a:t < +oo}
...
The geometric choice has been widely used in the literature since it produces a linear sequence of pulsations on a Bode plot. In practise, it realises a good compromise, between the width of the band of interest and the 'Gomplexity of the resulting approximation. Such a choice has been made, for instance, by Alain Oustaloup and co-authors (see e.g. [Ous83]) when fitting phase and magnitude of fractional integrals or derivatives over a bounded range of frequencies. But, it happens that this choice is not unique and many others are to be investigated. Examples and simulations are given in subsection 3.5.
(15)
Then, the proof of theorem 3 can be sketched as follows: • decompose IIX - Xliv", into three parts in the frequency domain, • choose m such as to bound IIX - Xliv", over the intervals ]0, ~:in] and [~~a:t' +oo[ by using a Lebesgue-type argument, • choose n such as to bound IIX - Xliv", over the fixed compact interval nm = [~~in'~~a:t] by using the convergence property of the interpolation process, • finally deduce convergence of the output y towards y in L2(0, T ; JR).
step 3: choice of the interpolating functions Finite elements interpolation corresponds to a standard method in numerical analysis (see [Sain96, chapter 3]); the better the quality of the approximation (convergence, speed, regularity) , the greater the order of the finite element. For simplicity sake, finite elements of order 1 have been chosen:
Hence, natural convergence conditions write out: 'om.n
cm .
m-.+oo -7
0
(16)
cm
m-.+oo -7
+00
(17)
'omax
Vm EN sup I~~I - ~kl n~oo 0
Note that the interpolating function Ak is null outside the interval [~k-I' ~k+l] and Ak(~k) = 1.
(18)
I~k~n
step 4: continuation of interpolation Clearly, the quality of the approximation depends on the quality of the reconstruction of the state X by means of the interpolating functions Ak. Since this reconstruction has to be performed over the whole frequency domain, it follows that the interpolating functions Al and An must have a particular shape, namely:
3.3 Approximation Process
Suppose that the band of interest is of the form [Wmin, wmax ] with Wmin > O. The question is then to build the network of poles N;: to realise a convergent approximation over this bounded range of frequencies. The construction of N;: is entirely based upon the convergence result and can be seen as a recursive process which stops as soon as the resulting approximation meets the desired accuracy. step 0: initialization m = O. Fix ~in = Wmin and ~a:t
Indeed, let X(~, t) be the state obtained for a given input u and denote X(~, 5) its Laplacetransform in 'Re(s) > O. Using (4), it comes:
= w max . 230
(21)
which implies:
X(~, t) (~o P,Q(~) Jot u(r) dr
!
X(~, t) (~~oo P,Qi~)
i
i t
(22)
~
r
iI
u(t)
l
Fig. 1. Bode plot of iI with a geometric choice of poles. Number of poles n = 20, [Wmin , Wmax ] = [10 2 ,10 5 ] . There is a good compromise between the size n of the system and the bandwidth.
step 5: accuracy criterion Either the approximation provides the desired accuracy or not. In the latter case, one should first increase n to coarsen the grid in om . Go back to step 2. IT accuracy is not improved then change ~:in and ~:ax according to conditions (16) and (17). For instance ~mtn cm+l = 110 'm cm.In and ":.maxcm+ 1 10 ~:ax' Then go back to step 1. Convergence of the approximation ensures that one will get the desired accuracy after a finite number of trials.
3.5 Numerical simulations
This subsection is devoted to numerical simulations of fractional integration of order Q = 0.2 over a bounded range of frequencies [Wmi n, w max ]. Two different choices of poles location have been examined, geometric and polynomial. For each cho~ce , there is one figure showing the Bode plot of H and another one describing its zeros location.
3.4 Constructive choice of the poles
Bode plots Figures (1) and (2) show Bode plots of iI for a geometric and a polynomial choice of poles respectively. They both fit the behaviour of the fractional integrator over a bounded range of frequencies. For the polynomial choice, there are not enough poles in the range of frequencies (n too small) , which gives rise to oscillations on the phase plot. Conversely, for the geometric choice, n is great enough to avoid oscillations. <>
A constructive choice for the network of poles
N:;' is presented here in the case of a geometric sequence. For a fixed order of approximation m, condition (18) writes out:
11
sup I~fl n~oo 0
·20
- 60
is completely determined.
Ir(n) -
-60
~
Recalling definition (9) of X(~, t) yields the desired result. Note that the final expression of A1 and An is determined by continuity in 6 and ~n respectively. At this stage,
-
(23)
lSk~n
with r(n) = (~:ax/~:in)1/(n-1) . For a given m , (23) is straightforward. Now, as m increases, om+1 ::> om. Therefore, n(m + 1), the number of points in om+l , should be greater than n(m) in order to get at least the same accuracyj it yields the condition:
Zeros location Zeros of the transfer function are not chosen. They appear as a consequence of the location of the poles and of the choice of matrices B and C. From simulations, poles and zeros alternate on the real axis 1R- and are asymptotically distributed in a regular way depending exclusively on Q and not on the AkSj
<>
Hence, two choices are proposed: ~:in cc
m-It
~:ax cc ml2 n(m) cc m IZ + E
or
zf E [~k' ~k-rl]
~:in ccexp(-mrl)
~:ax cc m r2 n(m) cc mrl+r2+E
with
zf denoting the kth
k
= 1, .. . ,n -
1
(25)
zero. Location of zf in
[~k ' ~k+l] can be described by the following local
where Lt, r1, 12 , r2 denote the speed of convergence of ~:in and ~:ax towards 0 and +00 respectively and where to > 0 must be chosen as small as possible to keep n to a reasonable value. However, other choices can be investigated.
slopes on Bode plot:
Ri ght SIope
231
m _ log(zk) -log(~k) k log( ~k+l) - log( ~k )
(26)
0 .•
0.•
0.7
t.
i
~
i°.5 a:
o
0
0.'
0.3 Q.2
o ... .
....... ~ . . . . . ~ .-.
O"!-D--=----::-----:----:--~s--7----:~---:---'---..J,D
"'* '" poIK
FrtlQ\WlCy(l'IIdIMe)
Fig. 4. Right (0) and left (+) local slopes with a polynomial choi'ce of poles. Number of poles n = 10. '.
Fig. 2. Bode plot of if with a polynomial choice of poles. Number of poles n = 10, [Wmin, w max ] = [10 3 ,10 6 ]. The number of poles n is too small to describe the range of frequencies and a careful look at the graph shows little oscillations on the phase plot.
5. REFERENCES [MAM98b] G. Montseny, J. Audounet, and D. Matignon. Diffusive wave-absorbing feedback controls for the boundary stabilization of a flexible beam. submitted for publication. [MAMb93] G. Montseny, J . Audounet, and B. Mbodje. Optimal models of fractional integrators and application to systems with fading memory. In Conference IEEE Systems, Man and Cybernetics, Le Touquet, France, 1993. IEEE-SMC . [MAM97] G. Montseny, J . Audounet, and D. Matignon. Fractional integrodifferential boundary control of the Euler-Bernoulli beam. In Conference on Decision and Control, pages 4973-4978, San Diego, California, December 1997. IEEE-CSS , SIAM. [MAM98a] D. Matignon, J. Audounet, and G. Montseny. Energy decay for wave equations with damping of fractional order. In Fourth international conference on mathematical and numerical aspects of wave propagation phenomena, 3 pages, Golden, Colorado, June 1998. INRIA, SIAM. [Mat96] D. Matignon. Recent results in fractional differential systems theory. Technical Report 96 C 004, Ecole Nationale Superieure des Telecommunications, 1996. [MbM95] B. Mbodje and G. Montseny. Boundary fractional derivative control of the wave equation. IEEE Trans. Aut. Cont., 40(2):378-382, Feb. 1995. [Ous83] A. Oustaloup. Systemes asservis lineaires d'ordre /ractionnaire. Automatique. Masson, 1983. [Sain96] L. Sainsaulieu. Calcul scientifique. Masson, 1996.
0.' 0
0.35
1 0.3
I
~D.2S
If
o 0.2
+
~- - . • ... • -.• --. - .- - . -. -. - .• .. ~ .- . - . - 9 o · +
. . 0
o 0.15
0.10~---7---:---;-----7---:'':-D--:'::,,---',.':---,:'::-.--'-".....:o~20 indexafpatn
Fig. 3. Right (0) and left (+) local slopes with a geometric choice of poles. Number of poles n = 20. Left SIope
m _
k
-
10g(ZM-l) -log(ek+l) log(zk+l) -log(zk)
(27)
Figures (3) and (4) show the value of the right and left slopes in the geometric and polynomial case respectively. They are both located around Cl = 0.2 in the range of approximation whatever the choice of the poles.
4. CONCLUSION A powerful method of approximation of fractional integrals based upon diffusive realisations has been thoroughly examined; it proves to be very flexible and its extension to fractional derivatives and other pseudo-differential operators is shortly in prospect. Hence it will be applied to the simulation of Fractional Differential Equations. 232