On the estimation of measurement errors in linear dynamical systems

Abstract

The article is concerned with the problem of linear dynamical system state estimation subject to noise-corrupted observations. In solution of the linear optimal ¯ltering problem under guaranteed statement the problem of measurement noise identi¯cation occurs. The approach to adaptive measurement noise estimation is proposed. It is based on statistical processing of the innovation sequence in the Kalman ¯lter. The special case of linear dynamical system with one-dimensional output the state of which is observed only on a short-time interval is considered and it is shown that statistical characteristics of the innovation sequence can be used for measurement noise identi¯cation and also for adjustment of the guaranteed estimates. Simulation results are given to con¯rm the usefulness of the approach.

Keywords

• guaranteed estimation
• Kalman ¯lter
• innovation sequence
• short sample of observations

Introduction

Control problems appear in various ¯elds of engineering: aircraft control systems, inertial navigation systems, automatic process control systems, tracking and target acquisition systems [1, 2]. When designing a dynamical object control system, it is necessary to estimate a state of object operating under conditions of a priori uncertainty and incomplete measurement data.

One of the basic requirements to dynamical object control system design is e±cient algorithms development for object current state estimation. In its turn, the state estimation is impossible without taking into account the combination of random disturbances that in°uence an object and measurement errors. Depending on assumptions of the nature of uncontrolled factors, the estimation problem is solved either under stochastic statement [3, 4], assuming a priori knowledge of statistical information, or under guaranteed statement [6, 7, 10, 13, 15, 16, 17, 19, 21], when only possible ranges of uncontrolled factors are known.

A common problem is the state estimation problem of dynamical systems in the presence of disturbances that can be both probabilistic and deterministic [10]. On the one hand, implementation of the Kalman ¯lter (KF) as a probabilistic estimation technique requires complete a priori knowledge of noise statistics. For the only measurement realization noise statistics cannot be obtained, are unknown or only their rough estimates are known. On the other hand, the main problem with guaranteed approach is that the state vector estimate obtained as a result of its implementation can be overestimated. This is explained by the fact that a solution of the mininimax ¯ltering problem is chosen taking into account the worst combination of uncertain factors (although those is very unlikely event). Probabilistically guaranteed approach can be used to overcome speci¯ed di±culties of the KF and the guaranteed algorithm.

In the present article the ¯ltering problem is considered under probabilistically guaranteed statement when the KF and the guaranteed algorithm are applied together. The estimation problem is solved using the guaranteed algorithm. The Kalman ¯lter implementation is performed for measurement data preprocessing, particularly, for measurement errors estimation. The special case of linear dynamical system with Gaussian random inputs and one-dimensional output is considered. 2

Statement of the problem

Consider the state estimation problem of a discrete linear dynamical system with the state and measurement di®erence equations:

xk+1 = Axk + ¡wk; yk+1 = Gxk+1 + vk+1; k = 0; 1; : : : ; N ¡ 1: where k is a discrete-time variable that takes values on a short interval k = 1; N , N < 30; xk 2 Rn is the state vector; wk 2 Rn is the process noise vector; yk 2 R1 is the measurement output; vk 2 R1 is the measurement noise; state transition matrix A, constant input matrix ¡ and constant output matrix G are known.

A priori information about the initial state x0, input disturbances wk and measurement errors vk is speci¯ed by sets of their possible values: ( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) ( 6 )

The guaranteed (or set-membership, minimax) estimation of the state vector xk involving the linear model ( 1 ), observations ( 2 ) and boundary conditions ( 3 ) assumes recursive construction of the bounded sets X¹k+1 (information sets, feasible sets), k = 0; 1; : : : ; N ¡ 1, i.e. the sets of state vector possible values. The set X¹k+1 contains the true value of the state xk+1 and its size depends largely on given sets ( 3 ): where Xk+1=k is the predicted state set and X[yk+1] is the measurement consistent set xk+1 2 X¹k+1 = Xk+1=k \ X[yk+1];

Xk+1=k = AX¹k + ¡W;

X[yk+1] = fx 2 RnjGx = yk+1 ¡ v; 8v 2 V g: All operations ( 4 ){( 6 ) in the guaranteed algorithm are performed on sets: set intersection, linear mapping of sets, Minkowsky sum, which in turn requires more computational power for implementation of the algorithm in real-time mode [13].

When a process is implemented, a priori given set of measurement errors V can be exceeding. For instance, measurement errors can be actually realized from a subset vk 2 V~ ½ V . Therefore, the problem of measurement noise identi¯cation in observations ( 2 ) occurs. The use of the set V~ instead of a priori given set in Eq. ( 6 ) allows to enhance a solution accuracy of the minimax ¯ltering problem. 3

An adaptive Kalman ¯lter

The stochastic approach to the state vector estimation problem in the system ( 1 ), ( 2 ) is to assume that the initial state x0 is a random variable with known mean value E[x0] = x^0 and known covariance matrix P0, the process noise wk and the measurement noise vk are uncorrelated white zero mean noise processes with covariance matrices Q and R respectively: x0 » N (0; P0); wk » N (0; Q); vk » N (0; R):

( 7 ) Under these circumstances the estimation problem involving the linear model ( 1 ), observations ( 2 ) and initial conditions ( 7 ) is solved by the KF [3, 4]: x^k+1 = x^k+1=k + Kk+1(yk+1 ¡ Gx^k+1=k); x^k+1=k = Ax^k; k = 0; 1; : : : ; N ¡ 1; where

Kk+1 = Pk+1=kGT[GPk+1=kGT + R]¡1;

Pk+1=k = APkAT + ¡Q¡T;

Pk+1 = (I ¡ Kk+1G)Pk+1=k: Here, x^k+1 is an estimate of the state vector xk+1 at the moment of time k + 1, Kk+1 is the ¯lter gain, Pk+1=k and Pk+1 are the covariance matrices of predicted and updated ¯ltering error respectively, I is n £ n identity matrix.

The sequence of estimates x^N (²) = fx^1; : : : ; x^N g obtained by the ¯lter Eqs. ( 8 ){( 11 ) is optimal in terms of minimum mean square error (MSE) only for a great number of measurement realizations [5, 14]. It is known that the true value of the state vector xk with required probability will belong to the set

Ek = fx 2 Rnj(xk ¡ x^k)TPk¡1(xk ¡ x^k) · l2g; which is an ellipsoid centered at the point x^k (the probability that xk 2 R2 can be found in ellipse with l = 3 is 0.989).

Implementation of the KF requires complete a priori knowledge of noise statistics ( 7 ) in Eqs. ( 1 ), ( 2 ) [5, 14]. In certain real systems data processing is performed for the only measurement realization. In this case, statistics of the process and measurement noise cannot be obtained, the process noise covariance matrix Q and the measurement noise covariance matrix R are unknown or only their rough estimates are known. In this case an adaptive Kalman ¯ltering approach can be used to solve the estimation problem [8, 9, 11, 12, 18, 20].

The adaptive ¯ltering algorithm is based on statistical analysis of the innovation sequence [8, 9, 11, 12] which is a zero mean Gaussian white noise process with correlation properties

¤k+1 = yk+1 ¡ Gx^k+1=k; Ck = ( GP 0 GT + R;

G[A(I ¡ KG)]k¡1A[P 0 GT ¡ KC0]; k = 0 k > 0 where P 0 is the covariance matrix of predicted ¯ltering error of time invariant ¯lter, K is the Kalman gain.

To verify the linear model ( 1 ) and observations ( 2 ) (i.e. to check whether the KF constructed using a priori given covariance matrices Q and R is close to optimal or not) it is possible to use correlation analysis methods for processing the sequence ¤N (²) = f¤1; : : : ; ¤N g. The estimates of Ck, denoted as C^k, can be obtained by using the ergodic property of a stationary innovation sequence for lag l = 1; 2; : : : ; n, n is the state vector dimension [9]:

The estimates C^k can be represented by rewriting ( 15 ) explicitly[8]:

C^k = 1 N

X ¤k¤kT¡l:

N k=l C1 = GAP 0 GT ¡ GAKC0; C2 = GA2P 0 GT ¡ GAKC1 ¡ GA2KC0; : : : Cn = GAnP 0 GT ¡ GAKCn¡1 ¡ : : : ¡ GAnKC0: ( 8 ) ( 9 ) ( 10 ) ( 11 ) ( 12 ) ( 13 ) ( 14 ) ( 15 ) ( 16 )

Using ( 14 ){( 16 ) the correlation of terms of the sequence ¤k is checked. In other words, according to the Eqs. ( 14 ){( 16 ) an optimality of implemented KF is checked. However, this test requires considerable computing time for the accumulation of innovation sequence statistic, i.e. to obtain C^k ( 15 ). In general, the test is carried out at the interval of available observations N À 100 [9]. For an optimal KF, the estimates C^k are unbiased and consistent, Ck vanishes for all k > 1. In the case of short measurements realization, e.g. when N < 30 measurements are taken, the estimates of Ck are inconsistent. If N is small, other tests can be used. 4

An adaptive algorithm of measurement noise estimation For an optimal ¯lter, when the process noise covariance matrix Q and the measurement noise covariance matrix R used in the ¯lter Eqs. ( 8 )-( 11 ) correspond to the real noises in the system ( 1 ), ( 2 ), a residual sum of squares (RSS) approaches zero [11, 12]

N¡1 X ¤kT+1¤k+1 ! 0: k=0

The case of biased innovation sequence E[¤k] 6= 0, or if the actual covariance of innovation sequence ¤k is substantially greater than its expected value ¤kT+1¤k+1 > GPk+1=kGT + R, can indicate the suboptimality of the KF [11, 12]. Detection of an innovation sequence bias or deviation of its actual covariance from expected covariance is carried out by averaging over some interval of the innovation sequence ¤N (²) = f¤1; : : : ; ¤N g . It is assumed that a number of sample points N may be low (N = 5 : : : 10). The required number of sample points mainly depends on time interval length on which it can be assumed that the ¯lter has reached steady-state conditions.

In order to obtain characteristics of the actual estimation, a posteriori value of innovation sequence can be used where ¢¹ = N1 PN

k=1 ¢k is the innovation sequence mean value.

The obtained estimate R^ = ¾^¢2 (23) can be used for construction of set V~ ½ V of measurement errors vk in the minimax ¯ltering algorithm ( 6 ). The set V~ is de¯ned as follows

N ¾^¢2 = 1 X(¢k ¡ ¢¹ )2;

N

k=1 vk 2 V~ = [v; v] = [¡lpR^;

+lpR^]; with covariance

¢k+1 = yk+1 ¡ Gx^k+1; varf¢k+1g = §k+1 = GPk+1GT + R: ( 17 ) ( 18 ) ( 19 ) ( 20 ) ( 21 ) (22) (23) (24) Substituting the observations ( 2 ) in ( 18 ) the expression for the innovation sequence can be rewritten ¢k+1 = Gxk+1 + vk+1 ¡ Gx^k+1 = Gek+1 + vk+1; where ek+1 = xk+1 ¡ x^k+1 denotes a vector of estimation errors in the KF. Then the expression for measurement errors vk can be de¯ned

vk+1 = ¢k+1 ¡ Gek+1: The covariance matrix Pk is the matrix of an ellipsoid for possible values of an estimation error vector ek: (xk ¡ x^k)TPk¡1(xk ¡ x^k) · l2 or ekTPk¡1ek · l2:

Checking an optimality of the KF ( 17 ), i.e. computing a posteriori value of innovation sequence ¢k for each time step k, we can obtain an estimate of the measurement noise covariance R^. For linear dynamical system ( 1 ) with one-dimensional measurement output ( 2 ), as an estimate R^ it is possible to use a variance of the innovation sequence where v, v are the lower and higher values of the measurement errors range. The values of implemented measurement errors vk with probability 0.997 (v 2 R1 l = 3) are in the set (24).

Thus, the adaptive algorithm of one-dimensional measurement noise estimation can be implemented. It is required to 1. Compute RSS ( 17 ) of the innovation sequence ¢k to verify an optimality of the ¯ltering process: ² if the implemented KF is not optimal, to obtain an estimate of R one can be used the known adaptive ¯ltering algorithm proposed by R.K. Mehra [8]; ² if the implemented KF is optimal, the estimate R^ can be obtained according to the (23). 2. De¯ne the bounded set V~ of measurement errors in the minimax ¯ltering algorithm according to the (24). 5

Numerical example

Consider implementation of proposed algorithm on an example of linear dynamical system ( 1 ), ( 2 ). It is assumed ·0:5¸ that the initial state x0 = 0:5 . The process noise wk and the measurement noise vk are normally distributed random numbers with zero means and standard deviations ¾w = 0:1 and ¾v = 0:2 respectively (Fig. 1, 2). The available N = 30 observations yk are shown in (Fig. 3).

0.2 0.1

0 −0.1 −0.2 −0.3 0 15 k 5 10 20 25 30 5 10 20 25

30 15 k 0.2

0 −0.2 −0.4 0 (25) (26) (27) (28) 20 25 30 5 10 20 25

30 15 k the information sets obtained by Xk+1=k \ X[yk+1]. This con¯rms that the use of the set V~ instead of a priori given set in Eq. ( 6 ) allows to enhance a solution accuracy of the minimax ¯ltering problem.

0 −0.2 −0.4 −0.6 15 k 0.5 0.4 X’[y ]

20 20/19

X’ 20/19 X[y ]

20 −1 0 1 0 2 1 0 −1

Summary and conclusions

The adaptive algorithm of one-dimensional measurement noise estimation is proposed. It is not necessary to assume the model of measurement errors. The algorithm is based on statistical analysis of the innovation sequence values in the Kalman ¯lter. The state vector estimation problem is considered on an example of linear dynamical system the state of which is observed only on a short-time interval. The estimation problem is solved under probabilistically guaranteed statement. It is assumed that the Kalman ¯lter and the guaranteed algorithm are applied together. The Kalman ¯lter implementation is performed for measurement data preprocessing. A numerical example is given to illustrate the results of proposed algorithm.

The work was supported by Act 211 Government of the Russian Federation, contract 02.A03.21.0011

References

[1] Y. Bar-Shalom , X.R. Li , T. Kirubarajan . Estimation with applications to tracking and navigation: theory, algorithms and software . John Wiley and Sons, 2001 .

[2] F. Auger , M. Hilairet , J.M. Guerrero , E. Monmasson, T. Orlowska-Kowalska , S. Katsura . Industrial applications of the Kalman ¯lter: a review . IEEE Transactioncs Industrial Electronics , 60 ( 12 ): 5458 { 5470 , 2013 .

[3] R.E. Kalman . A new approach to linear ¯ltering and prediction problems . Transactions of the ASME { Journal of Basic Engineering , 82D : 35 { 45 , 1960 .

[4] R.E. Kalman , R.S. Bucy . New results in linear ¯ltering and prediction theory . Transactions of the ASME { Journal of Basic Engineering , 83D : 95 { 108 , 1961 .

[5] R.E. Kalman . Identi¯cation of noisy systems . Russian Mathematical Surveys , 40 ( 4 ): 25 { 42 , 1985 .

[6] F.C. Schweppe . Recursive state estimation: unknown but bounded errors and system inputs . IEEE Transactions on Automatic Control , 13 ( 1 ): 22 { 28 , 1968 .

[7] F.L. Chernousko . Minimax control for a class of linear systems subject to disturbances . Journal of Optimization Theory and Applications , 127 ( 3 ): 535 { 548 , 2005 .

[8] R.K. Mehra . On the identi¯cation of variances and adaptive Kalman ¯ltering . IEEE Transactions on Automatic Control , AC- 15 ( 2 ): 175 { 184 , 1970 .

[9] R.K. Mehra . Approaches to adaptive ¯ltering . IEEE Transactions on Automatic Control , 17 ( 5 ): 693 { 698 , October 1972 .

[10] I. Ya. Kats , A.B. Kurzhanskii . Minimax multi-step ¯ltering in statistically uncertain situations . Avtomat. i Telemekh , 11 : 79 { 87 , 1978 .

[11] A.H. Jazwinski . Adaptive Filtering. Automatica , 5 ( 4 ): 475 { 485 , July 1969 .

[12] Jazwinski , A.H. Stochastic Processes and Filtering Theory / A.H. Jazwinski . New York: Academic. 1970 .

[13] V.I. Shiryaev. Synthesis of control of linear systems in incomplete information . Izv . Ross. Akad. Nauk, Tekhn. Kibern., 3 : 229 { 237 , 1994 .

[14] L.A. Fokin , V.I. Shiryaev. Preliminary comparison of Kalman and minimax approaches to error estimation of integrated navigation system . IEEE International Siberian Conference on Control and Communications (SIBCON-2013) Proceedings , 3 : 212 { 214 , September 2013 .

[15] A.B. Kurzhanski . Identi¯cation { a theory of guaranteed estimates, in From data to model (ed . J.C. Willems). Springer-Verlag Berlin Heidelberg, 1989 .

[16] A.B. Kurzhanski , K. Sugimoto , I. Valyi. Guaranteed state estimation for dynamical systems: ellipsoidal techniques . International Journal of Adaptive Control and Signal Processing , 8 : 85 { 101 , 1994 .

[17] C. Combastel . Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence . Automatica , 55 : 265 { 273 , 2015 .

[18] M.M. Kogan . Optimal estimation and ¯ltration under unknown covariances of random factors . Autom. Remote Control , 75 ( 11 ): 1964 { 1981 , 2014 .

[19] V.T.H. Le , C. Stoica , T. Alamo , E.F. Camacho , D. Dumur . Zonotopic guaranteed state estimation for uncertain systems . Automatica , 49 ( 1 ): 3418 { 3424 , 2013 .

[20] A.E. Barabanov , Yu.A. Lukomskij , A.N. Miroshnikov . Adaptive ¯ltration under unknown intensities of disturbances and measurement noises . Autom. Remote Control , 53 ( 11 ): 1731 { 1738 , 1992 .

[21] T. Alamo , J.M. Bravo , E.F. Camacho . Guaranteed state estimation by zonotopes . Automatica , 41 ( 6 ): 1035 { 1043 , 2005 .