Discrete fractional Kalman filter

Discrete fractional Kalman filter Moufida ABDELHAMID ∗ , Mohamed AOUN ∗∗ , Slaheddine NAJAR, Mohamed Naceur ABDELKRIM National Engineers School of Gabes, Road of Medenine 6029 Gabes, Tunisia. (e-mail: [email protected]). ∗∗ National Engineers School of Gabes, Road of Medenine 6029 Gabes, Tunisia. (e-mail: [email protected]) ∗

Abstract: This paper presents a generalization of the classical discrete time Kalman filter algorithm to the case of the fractional systems. Motivations for the use of this filter are given and the algorithm is detailed. The document also shows a simple numerical example of linear state estimation. Keywords: Fractional differentiation; Fractional Kalman filter; state estimation. 3. DISCRETE FRACTIONAL KALMAN FILTER

1. INTRODUCTION The concept of differentiation to an arbitrary order (fractional differentiation) was defined in the 19th century by Riemann and Liouville. Their main concern was to extend differentiation by using not only integer but also noninteger orders. In the field of system identification studies on real systems, such as electrochemical (Ichise and al. 1971; Darling and Newman 1997) or thermal (Battaglia and al. 2001; Cois and al. 2003), reveal inherent fractional differentiation behavior. The use of classical models (based on integer order differentiation) is thus inappropriate in representing these fractional systems. The filtering problem in the case of linear systems is solved by using the classic Kalman filter; in the case of nonlinear systems the extended Kalman filter is used. The two already quoted cases require an integer state representation. However, as it was specified previously, several studies on physical systems show that they are modelled by a fractional state space representation, thus the need to define fractional Kalman filter. Compared to (Sierocuik and Dzielinski, 2006), this paper brings a new representation of the fractional Kalman model. In section 2 a definition of fractional calculus is outlined. The generalization of the discrete fractional Kalman filter algorithm is developed in section 3. Section 4 presents numerical example for the state estimation. Finally a whiteness and validation test are shown. 2. FRACTIONAL DIFFERENTIATION The Gr¨ unwald fractional derivative of order ν of f (t) is defined as:   ∞ 1 X ν j ν f (t − jh) (1) (−1) D f (t) = lim ν j h→0 h j=0

where h is a sampling time and    1 ν = ν(ν−1)...(ν−j+1) j j!

f or j = 0 f or j > 0

(2)

3.1 Fractional Kalman model A linear fractional order stochastic system is define by the following state space representation:  D(ν) x(t) = Ax(t) + Bu(t) + Gw(t) (3) y(t) = Cx(t) + v(t) Where: ν is the non integer order, w(t) and v(t) are state and output noises. A discretisation, of the above equations, allows to write, by applying the Gr¨ unwald definition:   ∞ 1 X ν ν j [D x(t)]t=Kh = lim ν (−1) x((K − j)h) (4) j h→0 h j=0

For h sufficiently small, the discrete state equation is written, according to:   ∞ X ν j (−1) xk+1−j = hν [Axk + Buk + Gwk ] (5) j j=0

For k < 0, the state is supposed null, the infinite sum in (5) is truncated to k + 1. By isolating the xk+1 term: xk+1 = hν [Axk +Bu k + Gwk ] + νxk k+1 P ν − (−1)j xk+1−j j j=2

(6)

Equation (6) could be rewritten as follows:   k+1 X ν j (−1) xk+1−j (7) xk+1 = Aν xk +Bν uk +Gν wk − j j=2

where: Aν = (hν A + νI), Bν = hν B and Gν = hν G. Thus the fractional Kalman model could be obtained by:

 xk+1 = Aν xk + B + Gν wk ν uk     k+1 P ν j xk+1−j − (−1) j  j=2   yk = Cxk + vk

In fact for the classic filter (ν = 1), the predicted value: (8)

The orders νk in (8) are all identical and equal to the commensurate order. For more generalization, orders of differentiation of the elementary differential equations of the fractional system can be considered distinct real numbers, ν = [ν1 . . . νN ]. N is the number of the elementary differential equations of the system. The Kalman model is then given by:  xk+1 = Aν xk + Bν uk + Gν wk    k+1 P j − (−1) Υj xk+1−j  j=2   yk = Cxk + vk where:



  hν1 0 0 ν1    Aν =  0 . . . 0 A +  0 0 0 hνN 0 ν Bν = h B Gν = hν G

(a)

(9)

(b)  0 0  .. . 0  0 νN

Υk = diag

ν1 k

...



νN k

−The prediction error covariance matrix The prediction error covariance matrix of the fractional Kalman filter is deduced, as for integer case by the following equation: h i T P˜k = E (˜ xk − xk ) (˜ xk − xk ) (15)

The prediction error (˜ xk − xk ) is written by using the equations (9.a) and (14): x ˜k − xk = Aν xˆk−1 + Bν uk−1 − 

 νj , k

j = 1...N . 

− does not depend on smoothing terms, thus we can simplify x ˜k+1 and without a great loss of the filter quality, it is legitimate to approximate the smoothed states by those estimated according to: ∗   E [xk+1−j |zk∗ ] ∼ =x ˆk+1−j (13) = E xk+1−j zk+1−j Using this assumption the following equation is obtained: k+1 X j x ˜k+1 = Aν xˆk + Bν uk − (−1) Υj x ˆk+1−j (14) j=2

and Υk is a diagonal matrix formed by the factor 

− depends only on the previous estimated value, the sum in (12) is null (for ν = 1, Υj = 0 ∀j > 2).



− Aν xk−1 − Bν uk−1 − Gν wk−1 +

j

(−1) Υj xk−j (16)

thus:

2

(10)

x ˜k − xk = Aν (ˆ xk−1 − xk−1 ) − Gν wk−1 k P − (−1)j Υj (ˆ xk−j − xk−j )

Step of prediction −State vector prediction The state vector prediction is written, in function of the previous state vector estimation (Brown and Hwang, 1997), (Haykin, 2001). x ˜k+1 = E [xk+1 |zk∗ ]
= E[ Aν xk + Bν uk + Gν wk (11) k+1 P − (−1)j Υj xk+1−j | zk∗ ] j=2

(12)

The covariance matrix P˜k is given by: P˜k = E[(˜ xk − xk )(˜ xk − xk )T ] = E[(Aν (ˆ xk−1 − xk−1 ) − Gν wk−1 k P − (−1)j Υj (ˆ xk−j − xk−j ))

(18)

(Aν (ˆ xk−1 − xk−1 ) − Gν wk−1 k P − [(−1)j Υj (ˆ xk−j − xk−j )])T ] j=2

Thus, we can write: h i T P˜k = Aν E (ˆ xk−1 − xk−1 ) (ˆ xk−1 − xk−1 ) Aν T   T T −Aν E (ˆ xk−1 − xk−1 ) wk−1 Gν h i k P j T −Aν (−1) Υj E (ˆ xk−1 − xk−1 ) (ˆ xk−j − xk−j ) j=2 h i   T T T −Gν E wk−1 (ˆ xk−1 − xk−1 ) Aν T + Gν E wk−1 wk−1 Gν h i k P j T +Gν (−1) Υj E wk−1 (ˆ xk−j − xk−j ) j=2

− +

j=2

From equation (12), note that the state vector prediction depends on the whole past of the state. It is a characteristic of the fractional Kalman filter.

(17)

j=2

j=2

As for the classic case the fractional Kalman filter is obtained in two steps: a step of prediction and a step of correction.

x ˜k+1 = Aν E [xk |zk∗ ] + Bν uk k+1 P j − (−1) Υj E [xk+1−j | zk∗ ]

k X j=2

Kalman filter is a stochastic, recursive estimator, which estimates the state xk of a system based on the knowledge of the series of measures zk∗ = {y0 , y1 ..., yk }, using the minimum conditional variance criterion: i h

Thus:

j

(−1) Υj x ˆk−j

j=2

3.2 The fractional Kalman equations

E kxk − xˆk / {y0 , y1 , ..., yk } 6 h i 2 E kxk − zk / {y0 , y1 , ..., yk }

k X

+

k P

j=2 k P

h i j T (−1) Υj E (ˆ xk−j − xk−j ) (ˆ xk−1 − xk−1 ) Aν T

j=2 k P j=2

  T j T (−1) Υj E (ˆ xk−j − xk−j ) wk−1 Gν

h i T Υj E (ˆ xk−j − xk−j ) (ˆ xk−j − xk−j ) ΥTj

(19)

Using the same assumptions of the classic Kalman filter, we can write:   T E (ˆ xk−1 − xk−1 ) wk−1 =0 h i T E wk−1 (ˆ xk−1 − xk−1 ) =0

h i T ∀j = 2...k, E wk−1 (ˆ xk−j − xk−j ) =0

(20)

  T ∀j = 2...k, E (ˆ xk−j − xk−j ) wk−1 =0 The equation (19) is reduced to: h i T P˜k = Aν E (ˆ xk−1 − xk−1 ) (ˆ xk−1 − xk−1 ) Aν T h i k P −Aν (−1)j Υj E (ˆ xk−1 − xk−1 ) (ˆ xk−j − xk−j )T j=2   T T Gν +Gν E wk−1 wk−1 h i k P j T − (−1) Υj E (ˆ xk−j − xk−j ) (ˆ xk−1 − xk−1 ) Aν T +

j=2 k P

j=2

h i T Υj E (ˆ xk−j − xk−j ) (ˆ xk−j − xk−j ) ΥTj

(21) We note that the sums, presenting the whole past, in the expression of the prediction error covariance matrix. Here too it is a characteristic of fractional Kalman filter. Sums are null for ν = 1, because (Υj = 0, ∀j > 2). The evaluation of intercorrelation terms of the estimation error between two distinct moments, namely E [(ˆ xl − xl )(ˆ xm − xm )] is very difficult. Without a great loss of the filter quality these intercorrelations are supposed null (estimation error at two distinct moments are independent), thus: h i E (ˆ xk−1 − xk−1 ) (ˆ xk−j − xk−j )T = 0 (22)

and

h i E (ˆ xk−j − xk−j ) (ˆ xk−1 − xk−1 )T = 0

(23)

As a result the prediction error covariance matrix is expressed by: P˜k = Aν Pk−1 Aν T + Gν Qk−1 Gν T +

k X

Υj Pk−j ΥTj (24)

j=2

that the minimization of previous criterion is equivalent to minimizing the following: x ˆk = arg min[(˜ xk − x)P˜k−1 (˜ xk − x)T x (27) +(yk − Cx)Rk−1 (yk − Cx)T ] By deriving the criterion (27), according to x, it is shown that the estimated state is a linear combination of the predicted state at the moment k and the innovation which is the difference between real measure and the predicted mesure at the same moment. The state vector estimation is written according to: (Bertein and Ceschi, 2005) x ˆk = x ˜ k + K k ( yk − C x ˜ ) (28) | {z k} Innovation

where Kk is the Kalman gain.

−The estimation error covariance matrix The estimation error covariance matrix is deduced as for the prediction error covariance matrix. Pk = E[(ˆ xk − xk )(ˆ xk − xk )T (29) The difference between the estimated state and the real state at one moment k is: x ˆk − xk = x˜k + Kk (yk − C x ˜k ) − xk (30) = (I − Kk C) x ˜k + Kk yk − xk By substituting yk by its expression (Cxk + vk ), equation (30) could be rewritten as follows: x ˆk − xk = (I − Kk C) x ˜k + Kk Cxk + Kk vk − xk (31) consequently: x ˆk − xk = (I − Kk C) (˜ xk − xk ) + Kk vk (32) Pk , the estimation error covariance matrix is then: Pk = E[(ˆ xk − xk )(ˆ xk − xk )T ] = E[((I − Kk C)(˜ xk − xk ) + Kk vk ) ((I − Kk C)(˜ xk − xk ) + Kk vk )T ] = (I − Kk C)E[(˜ xk − xk )(˜ xk − xk )T ](I − Kk C)T T T +Kk E[vk vk ]Kk (33) Rk is a covariance matrix of an output noise vk in (9.b). Rk = E[vk vkT ] as a result we obtain: T Pk = (I − Kk C) P˜k (I − Kk C) + Kk Rk KkT

(34) (35)

−Kalman gain

Qk is a covariance matrix of a system noise wk in (9.a). Qk = E[wk wkT ] (25) Finally, according to (14) and (24), the equations of the step of prediction are: k+1 P x ˜k+1 = Aν x ˆk + Bν uk − (−1)j Υj x ˆk+1−j j=2

k P Υj Pk−j ΥTj P˜k = Aν Pk−1 Aν T + G−ν Qk−1 G−ν T + j=2

(26)

Step of correction −The state vector estimation The estimated state is obtained by minimizing the criterion (10). Indeed in (Schutter and al., 1999) one shows

Kalman gain Kk is obtained in order that, the estimation error is statistically orthogonal with the innovation (yk − C x ˜k ) (Auger, 1999):   E (ˆ xk − xk ) (yk − C x ˜k ) T = 0

By substituting xˆk and yk by their expressions (˜ xk + Kk (yk − C x ˜k )) and (Cxk + vk ): E [((I − Kk C)(˜ xk − xk ) + Kk vk ) (36) (Cxk + vk − C x ˜k ) T ] = 0 then: E [((I − Kk C)(˜ xk − xk ) + Kk vk ) (37) (−C(˜ xk − xk ) + vk ) T ] = 0 Taking account of the independence of the prediction error of the state (˜ xk − xk ) and the sequence of noise vk :   E (˜ xk − xk ) vk T = 0

the equation (37) becomes: h i T (I − Kk C) E (˜ xk − xk ) (˜ xk − xk ) C T (38)   −Kk E vk vkT = 0 also: − (I − Kk C) P˜k C T + Kk Rk = 0 (39) from where we extract the expression of the optimal Kalman gain:  −1 (40) Kk = P˜k C T C P˜k C T + Rk

Using (39) the estimation error covariance matrix(35) can be simplified. Indeed,

Kk Rk KkT = (I − Kk C) P˜k C T KkT (41) using equation (41), (35) is written as follows: T ˜ T T Pk = (I − Kk C) P˜k (I h − Kk C) + (I − Kk iC) Pk C Kk T T T ˜ = (I − Kk C) Pk (I − Kk C) + C K

A discretisation of (49) with the sampling time h = 1 leads to:      0 0.7 0.1995   u x + x =  k+1  0.063 k −0.0063 −1.1874 k  k+1 P (50) +wk − (−1)j Υj xk+1−j    j=2   yk = [ 0.1 0.3 ] xk + vk Let us suppose that: − wk and vk are two white noises characterized by:       0.003 0 and E vk vkT = 0.3. (51) E wk wkT = 0 0.003 − the initial state is a random variable independent of noises wk and vk characterized by: T

E [x0 ] = m0 =  [ 0 0 ]   0 0 E x0 wkT = 0 0   T T E x0 vk = [ 0 0 ]

k

(42) The simplified expression of the estimation error covariance matrix is then: (43) Pk = (I − Kk C) P˜k Finally, it helds from (28), (40) and (43) the step of correction equations: xˆk = x˜k + Kk (yk − C x ˜k ) Pk = (I − Kk C) P˜k (44)  −1 Kk = P˜k C T C P˜k C T + Rk

 h i  100 0 T P0init = E (x0 − m0 ) (x0 − m0 ) = 0 100

(53)

The system is excited by the following input signal. 1

Signal input

0.5

Discrete fractional Kalman filter:

0

The discrete fractional Kalman filter equations are written finally as follows:

(52)

0

50

100

150

200 Time(s)

250

300

350

400

Fig. 1. The input signal

− Prediction equations: x ˜k+1 = Aν x ˆk + Bν uk −

k+1 P

j

(−1) Υj x ˆk+1−j

j=2

k P P˜k = Aν Pk−1 Aν + G−ν Qk−1 G−ν T + Υj Pk−j ΥTj T

The noises wk and vk are generated respecting (51). The real and the noised output signals are traced on the figure (2).

j=2

(45)

2.5 Real output Noised output

2

(46)

1.5 1

(47)

The initial state and the initial covariance matrix of the system are: h i T x0 ∈ RN , P0init = E (x0 − m0 ) (x0 − m0 ) (48) 4. ACADEMIC EXAMPLE Let the fractional linear state space representation with two differentiation orders ν1 = 0.7 and ν2 = 1.2 (also treated in Sierocuik and Dzielinski, 2006).       ([0.7 0 1 0 D 1.2 ]) x(t) = x(t) + u(t) + w(t) −0.1 −0.2 1  y(t) = [0.1 0.3] x(t) + v(t) (49)

Outputs

− Correction equations: x ˆk = x ˜k + Kk (yk − C x ˜k ) Pk = (I − Kk C) P˜k where Kk is the optimal Kalman gain, given by:  −1 Kk = P˜k C T C P˜k C T + Rk

0.5 0 −0.5 −1 −1.5

0

50

100

150 200 Time(s)

250

300

350

Fig. 2. Evolution of the real and the noisy output signals The fractional Kalman filter previously described allows the estimation of the states variables x ˆ1 and x ˆ2 with high accuracy.

Real and estimated states are shown on the figure (3). This last shows that estimated states and real states are almost superimposed. 14 Real state x

1

12

Estimated state x

1

5. WHITENESS TEST By a whiteness test (computation of the ratio of normalized autocorrelation estimates) we have checked, that the intercorrelations leading to the assumptions (22) and (23) are white.

Real state x

2

10

Estimated state x

2

States

8 RN(i)

1

6 4

0.5

2 0

0

−0.5 −200

−2 −4

0

50

100

150 200 Time(s)

250

300

350

−150

−100

−50

0

50

100

150

200

Fig. 7. Computation of the ratio of normalized autocorrelation estimates

Fig. 3. Evolution of real and estimated states Figures (4) and (5), present respectively the estimation error and the prediction error on the output and the states. They are very weak and converge from a certain time tc ≃ 20s. 1.5

Estimation error output Prediction error output

Output errors

1 0.5 0 −0.5

0

50

100

150

200 Time(s)

250

300

350

400

Fig. 4. Estimation error and prediction error on the output 4

Estimation error on x

State errors

3

Estimation error on x

6. VALIDATION TEST The method consists in generating a group of measure yk for an input u, to which we associate each time two sequences of noises wk and vk . These last are used respectively to calculate the covariances matrix on the input and the output noises Q and R. So we can apply the fractional Kalman filter (N = 25) time to estimated the state x ˜k . The mean of the N state vectors estimation x ˜k is roughly the real state xk .

1 2

2

Prediction error on x

1

Prediction error on x2

2.5

1

Real outputs Estimated outputs Predicted outputs

2

0 1.5

−1 0

50

100

150

200 Time(s)

250

300

350

400

Fig. 5. Estimation error and prediction error on states x1 and x2

1 Outputs

−2

0.5 0 −0.5

The time of convergence (tc < 20s) can be checked by the analysis of the trace of the state error covariance matrix, shown in figure (6).

Error

200

−1.5

0

50

100

150 200 time(s)

250

300

350

Estimation error covariance matrix

Fig. 8. Evolution of different outputs

100

0

−1

0

50

100

150

200 Time(s)

250

300

350

Fig. 6. Evolution of the covariance matrix trace

400

Result of the N state estimation are shown in figure (9).

14 Differente estimated states 12 10 8

States

6 4 2 0 −2 −4 −6

0

50

100

150 200 Time(s)

250

300

350

Fig. 9. Evolution of different states As it could be seen in figure (10) the original output is estimated with high accuracy. 2 Real ouput Mean of the N prediction output 1.5

Outputs

1

0.5

0

−0.5

−1 −50

0

50

100

150 Time(s)

200

250

300

350

Fig. 10. Evolution of the real output and the mean of predicted output 7. CONCLUSION This paper establishes a generalization of the classical discrete time Kalman filter algorithm to the case of the fractional systems. The Kalman filter equations are modified to taking account of the fractional derivative order presented by the fractional system state equation. The detaily developped calculus computations are made and well justified and verified assumptions through the whiteness test and validation test with an academic fractional system example. REFERENCES Aoun, M. (2005). Syst`emes lin´eaires non entiers et identification par bases orthogonales non enti`eres. Th`ese de doctorat Universit´e Bordeaux 1 France. Arzelier, D. (2000). Introduction ` a la th´eorie d’estimation. ENSICA. Notes de cours, Version 1.2.

Auger, F. (1999). Introduction ` a la th´eorie du signal et de l’information: Cours et exercices volume 2. Collection sciences et technologies editions technip edition. Axtell, M. and Bise, E. (1990). Fractional calculus applications in control systems. In IEEE National Aerospace and Electronics Conference New York, USA. Bertein, J. and Ceschi, R. (1998). Processus stochastiques et filtrage de Kalman. Herm`es Paris. Bertein, J. and Ceschi, R. (2005). Processus stochastiques discret et filtrages optimaux. Number 271 pages. Herm`es science : Lavoisier. Brown, R. and Hwang, P. (1997). Introduction to Random Signals and Applied Kalman Filtering with Matlab Exercises and Solutions. New York: Wiley. Caputo, M. (1967). Linear models of dissipation whose q is almost frequency independent. Geophysical journal of the royal astronomical society 2(13):529–539. Cois, O. (2002). Syst`emes lin´eaires non entiers et identification par mod`ele non entier: application en thermique. Th`ese de doctorat Universit´e Bordeaux 1 France. Dugowson, S. (1994). Les Diff´erentielles m´etaphysiques : Histoire et philosophie de la g´en´eralisation de l’ordre de d´erivation. Th`ese de doctorat Universit´e Paris XIII France. Grewal, M. S. and Andrews, A. P. (2001). Kalman Filtering:Theory and Practice Using Matlab. John Wiley and Sons,Inc second edition. Haykin, S. (2001). Kalman Filtering and Neural Networks. New York: Wiley. Ichise, M., Y. Nagayanagi, and T. Kojima; (1971). An analog simulation of non integer order transfer functions for analysis ofelectrode processes. In J. Electroanal. Chem. Interfacial Electrochem pages 33–253. Kalman, R. (1960). A new approach to linear filtering and prediction problems. In Trans. ASME J. Basic Eng. volume 82 of D pages 35–45. Kalman, R. and Bucy, R. (1961). New results in linear filtering and prediction problems. In Trans. ASME J. Basic Eng. volume 83 of D pages 95–108. Miller, K.-S. and Ross, B. (1993). An introduction to the fractional calculus and fractional differential equations. A Wiley-Interscience Publication. Oldham, K.-B. and Spanier, J. (1974). The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order. Academic Press. Oustaloup, A. (1995). La d´erivation non enti`ere: th´eorie, synth`ese et applications. Herm`es Paris. Podlubny, I. (1999). Fractional Differential Equations. Academic press San Diego. Ross, B. (1974). A breif history and exposition of the fundamental theory of fractional calculus. Fractional calculus and its applications (457):1–37. Samko, S.-G., Kilbas, A.-A., and Marichev, O.-I. (1993). Fractional integrals and derivatives: theory and applications. Gordon and Breach Science. Sierocuik, D. and Dzielinski, A. (2006). Fractional kalman filter algorithm for the states,parameters and order of fractional system estimation. Int. J. Appl. Math. Comput. Sci. Schutter, J., Geeter, J., Lefebvre, T., and H., B. (1999). Kalman filters: A tutorial. http://citeseer.ist.psu.edu/443226.html.