On the estimation of measurement errors in linear dynamical systems Dina Khadanovich Vladimir Shiryaev dinakhadanovich@mail.ru shiriaevvi@susu.ru South Ural State University (Chelyabinsk, Russia) 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 filtering problem under guaranteed statement the problem of measurement noise identification occurs. The approach to adaptive measurement noise estimation is proposed. It is based on sta- tistical processing of the innovation sequence in the Kalman filter. The special case of linear dynamical system with one-dimensional output the state of which is observed only on a short-time interval is consid- ered and it is shown that statistical characteristics of the innovation sequence can be used for measurement noise identification and also for adjustment of the guaranteed estimates. Simulation results are given to confirm the usefulness of the approach. Keywords: guaranteed estimation; Kalman filter; innovation se- quence; short sample of observations. 1 Introduction Control problems appear in various fields of engineering: aircraft control systems, inertial navigation systems, automatic process control systems, tracking and target acquisition systems [1, 2]. When designing a dynam- ical 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 efficient 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 influence 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 filter (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 c by the paper’s authors. Copying permitted for private and academic purposes. Copyright ° In: G.A. Timofeeva, A.V. Martynenko (eds.): Proceedings of 3rd Russian Conference ”Mathematical Modeling and Information Technologies” (MMIT 2016), Yekaterinburg, Russia, 16-Nov-2016, published at http://ceur-ws.org 35 obtained as a result of its implementation can be overestimated. This is explained by the fact that a solution of the mininimax filtering 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 specified difficulties of the KF and the guaranteed algorithm. In the present article the filtering 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 filter 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 difference equations: xk+1 = Axk + Γwk , (1) yk+1 = Gxk+1 + vk+1 , k = 0, 1, . . . , N − 1. (2) where k is a discrete-time variable that takes values on a short interval k = 1, N , N < 30; xk ∈ Rn is the state vector; wk ∈ Rn is the process noise vector; yk ∈ R1 is the measurement output; vk ∈ 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 specified by sets of their possible values: x0 ∈ X0 , wk ∈ W, vk ∈ V. (3) 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): xk+1 ∈ X̄k+1 = Xk+1/k ∩ X[yk+1 ], (4) where Xk+1/k is the predicted state set Xk+1/k = AX̄k + ΓW, (5) and X[yk+1 ] is the measurement consistent set X[yk+1 ] = {x ∈ Rn |Gx = yk+1 − v, ∀v ∈ V }. (6) 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 ∈ Ṽ ⊂ V . Therefore, the problem of measurement noise identification in observations (2) occurs. The use of the set Ṽ instead of a priori given set in Eq. (6) allows to enhance a solution accuracy of the minimax filtering problem. 3 An adaptive Kalman filter 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) 36 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, (8) where Kk+1 = Pk+1/k GT [GPk+1/k GT + R]−1 , (9) Pk+1/k = APk AT + ΓQΓT , (10) Pk+1 = (I − Kk+1 G)Pk+1/k . (11) Here, x̂k+1 is an estimate of the state vector xk+1 at the moment of time k + 1, Kk+1 is the filter gain, Pk+1/k and Pk+1 are the covariance matrices of predicted and updated filtering error respectively, I is n × n identity matrix. The sequence of estimates x̂N (•) = {x̂1 , . . . , x̂N } obtained by the filter 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 = {x ∈ Rn |(xk − x̂k )T Pk−1 (xk − x̂k ) ≤ l2 }, (12) which is an ellipsoid centered at the point x̂k (the probability that xk ∈ 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 filtering approach can be used to solve the estimation problem [8, 9, 11, 12, 18, 20]. The adaptive filtering algorithm is based on statistical analysis of the innovation sequence [8, 9, 11, 12] Λk+1 = yk+1 − Gx̂k+1/k , (13) which is a zero mean Gaussian white noise process with correlation properties ( 0 GP GT + R, k=0 Ck = 0 (14) G[A(I − KG)]k−1 A[P GT − KC0 ], k>0 0 where P is the covariance matrix of predicted filtering error of time invariant filter, 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 (•) = {Λ1 , . . . , ΛN }. The estimates of Ck , denoted as Ĉ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]: N 1 X Ĉk = Λk ΛT k−l . (15) N k=l The estimates Ĉk can be represented by rewriting (15) explicitly[8]: 0 C1 = GAP GT − GAKC0 ; 0 C2 = GA2 P GT − GAKC1 − GA2 KC0 ; (16) ... 0 Cn = GAn P GT − GAKCn−1 − . . . − GAn KC0 . 37 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 Ĉk (15). In general, the test is carried out at the interval of available observations N À 100 [9]. For an optimal KF, the estimates Ĉ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 filter, when the process noise covariance matrix Q and the measurement noise covariance matrix R used in the filter Eqs. (8)-(11) correspond to the real noises in the system (1), (2), a residual sum of squares (RSS) approaches zero [11, 12] N X −1 ΛT k+1 Λk+1 → 0. (17) 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 ΛT T k+1 Λk+1 > GPk+1/k G + 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 (•) = {Λ1 , . . . , ΛN } . 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 filter has reached steady-state conditions. In order to obtain characteristics of the actual estimation, a posteriori value of innovation sequence can be used ∆k+1 = yk+1 − Gx̂k+1 , (18) with covariance var{∆k+1 } = Σk+1 = GPk+1 GT + R. (19) 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 , (20) 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 defined vk+1 = ∆k+1 − Gek+1 . (21) The covariance matrix Pk is the matrix of an ellipsoid for possible values of an estimation error vector ek : (xk − x̂k )T Pk−1 (xk − x̂k ) ≤ l2 or eT −1 2 k Pk ek ≤ l . (22) 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 N 2 1 X ¯ 2, σ̂∆ = (∆k − ∆) (23) N k=1 ¯ = 1 PN ∆k is the innovation sequence mean value. where ∆ N k=1 2 The obtained estimate R̂ = σ̂∆ (23) can be used for construction of set Ṽ ⊂ V of measurement errors vk in the minimax filtering algorithm (6). The set Ṽ is defined as follows p p vk ∈ Ṽ = [v; v] = [−l R̂; +l R̂], (24) 38 where v, v are the lower and higher values of the measurement errors range. The values of implemented mea- surement errors vk with probability 0.997 (v ∈ 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 filtering process: • if the implemented KF is not optimal, to obtain an estimate of R one can be used the known adaptive filtering algorithm proposed by R.K. Mehra [8]; • if the implemented KF is optimal, the estimate R̂ can be obtained according to the (23). 2. Define the bounded set Ṽ of measurement errors in the minimax filtering 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 = . The process noise wk and the measurement noise vk are normally distributed 0.5 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.4 0.1 0.2 0 0 −0.1 −0.2 −0.2 −0.3 −0.4 0 5 10 15 20 25 30 0 5 10 15 20 25 30 k k Figure 1: Process noise wk . Figure 2: Measurement noise vk . System matrices are µ ¶ µ ¶ 0.8 −0.7 0 ¡ ¢ A= , Γ= , G= 1 0 . (25) 0.4 0.8 1 The set X0 is specified by square X0 = {x ∈ R2 | − 1 ≤ x(1) ≤ 1, −1 ≤ x(2) ≤ 1}. (26) The sets Wk and Vk are specified by intervals W = {w ∈ R1 | − 0.3 ≤ w ≤ 0.3}, V = {v ∈ R1 | − 0.6 ≤ v ≤ 0.6}. (27) The covariance matrix P0 is defined in a way that the value of state vector x0 at a three sigma level is in the set X0 . Then the initial conditions for the KF are defined as follows: · ¸ 0 R2 2 x̂0 = , P0 = I, q = σw , r = σv2 , (28) 0 9 39 1.5 0.5 1 0.4 0.5 0.3 0 0.2 −0.5 0.1 −1 0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 k k Figure 3: Noise-corrupted observations yk . Figure 4: RSS processing result. √ ¡ ¢ where R = 2 is a radius of circle circumscribed around the set X0 , P0 = diag 0.2222, 0.2222 . Fig. (4) shows computational results of RSS. The results of the guaranteed algorithm are reported in Fig. (6): the sets X[yk ] and Xk+1/k are obtained as 0 a result of set operations (4)–(6), when as a set of measurement errors vk (3) the set V is used; the sets X [yk ] 0 and Xk+1/k are obtained as a result of set operations (4)–(6), when as a set of measurement errors vk (3) the set Ṽ is used. 0 0 As it is seen from Fig. (6), the information sets obtained as intersection Xk+1/k ∩ X [yk+1 ] are smaller than the information sets obtained by Xk+1/k ∩ X[yk+1 ]. This confirms that the use of the set Ṽ instead of a priori given set in Eq. (6) allows to enhance a solution accuracy of the minimax filtering problem. 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 0 5 10 15 20 25 30 k Figure 5: Measurement noise vk and innovation sequence ∆k (the solid lines denote vk and boundaries of the set V , the dashed lines denote ∆k and boundaries of the set Ṽ ). 40 2 2 X[y5] X[y15] X’[y5] 1 1 X’[y15] 0 y5 0 y15 x15 x X5/4 5 X’ 15/14 −1 −1 X’ X15/14 5/4 −2 −2 −2 −1 0 1 2 −2 −1 0 1 2 0 0 0 0 a) The sets X[yk+1 ], Xk+1/k , X [yk+1 ], Xk+1/k b) The sets X[yk+1 ], Xk+1/k , X [yk+1 ], Xk+1/k at the time step k = 5 at the time step k = 15 2 2 X’[y20] X[y30] X’[y30] 1 1 X X 20/19 30/29 x30 X’30/29 x20 y 30 0 y 0 20 X’20/19 −1 −1 X[y20] −2 −2 −2 −1 0 1 2 −2 −1 0 1 2 0 0 0 0 c) The sets X[yk+1 ], Xk+1/k , X [yk+1 ], Xk+1/k d) The sets X[yk+1 ], Xk+1/k , X [yk+1 ], Xk+1/k at the time step k = 20 at the time step k = 30 Figure 6: The adaptive minimax filtering results. The solid line denotes the sets X[yk+1 ] and Xk+1/k . The 0 0 dashed line denotes the sets X [yk+1 ] and Xk+1/k . 41 6 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 filter. 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 filter and the guaranteed algorithm are applied together. The Kalman filter 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 filter: a review. IEEE Transactioncs Industrial Electronics, 60(12):5458– 5470, 2013. [3] R.E. Kalman. A new approach to linear filtering 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 filtering and prediction theory. Transactions of the ASME – Journal of Basic Engineering, 83D:95–108, 1961. [5] R.E. Kalman. Identification 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 identification of variances and adaptive Kalman filtering. IEEE Transactions on Automatic Control, AC-15(2):175–184, 1970. [9] R.K. Mehra. Approaches to adaptive filtering. IEEE Transactions on Automatic Control, 17(5):693–698, October 1972. [10] I.Ya. Kats, A.B. Kurzhanskii. Minimax multi-step filtering 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 es- timation of integrated navigation system. IEEE International Siberian Conference on Control and Communications (SIBCON-2013) Proceedings, 3:212–214, September 2013. [15] A.B. Kurzhanski. Identification – 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. 42 [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 filtration 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 filtration 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. 43