# Receding Horizon Least Squares Estimator with Application to Estimation of Process and Measurement Noise Covariances

Authors: Vladimir Shin; Rebbecca T. Y. Thien; Yoonsoo Kim
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper presents a noise covariance estimation method for dynamical models with rectangular noise gain matrices. A novel receding horizon least squares criterion to achieve high estimation accuracy and stability under environmental uncertainties and experimental errors is proposed. The solution to the optimization problem for the proposed criterion gives equations for a novel covariance estimator. The estimator uses a set of recent information with appropriately chosen horizon conditions. Of special interest is a constant rectangular noise gain matrices for which the key theoretical results are obtained. They include derivation of a recursive form for the receding horizon covariance estimator and iteration procedure for selection of the best horizon length. Efficiency of the covariance estimator is demonstrated through its implementation and performance on dynamical systems with an arbitrary number of process and measurement noises.


---

## Body

## 1. Introduction

Filtering and system identification are powerful techniques for building models of complex control systems in time and frequency domains [1–5]. The Kalman filtering and its variations are signal processing techniques in widely used applications, such as navigation, target tracking, vehicle state estimation, communications engineering, air traffic control, and many others [6–9]. As the well-known standard Kalman filter (KF) equations do not directly depend on the original process and measurement noise covariances Q and R, respectively. They depend on the transformed covariances,(1)Q~=GQGT,andR~=DRDT,where G=Gk and D=Dk are the rectangular gain matrices associated with process and measurement noises, respectively [1, 2, 6, 7]. In this paper, the original Q,R and transformed Q~,R~ noise covariances shall be denoted by “ONC” and “TNC”, respectively.The optimality of KF depends on the accuracy of TNC, which directly depend on the prior assumptions on the ONC. Inadequacy of the prior values for the TNC can lead to unexpected results and the KF divergence [1, 2, 7]. The adaptive filtering is one of the approaches to prevent the divergence problem when precise knowledge on the noise covariances is not available.Adaptive Kalman filtering (AKF) is not the subject of this paper, so we present only its brief survey. The AKF represents a promising strategy for dynamical adjustment of the TNCQ~,R~ for online estimation. Two popular types of the AKF algorithms are the innovation-based estimation approach [10–13] and the adaptive fading filtering approach [14, 15]. The latter is a type of covariance scaling method, into which suboptimal fading factors are incorporated. The innovation-based estimation approach calculates Q~ and/or R~, assuming that the innovation sequence of the KF represents a white noise. In this approach, one of the covariance matching, the correlation, and the maximum likelihood methods is used. In the covariance matching method at first the sample error covariance of a state is computed based on a moving window and then Q~ is estimated using the calculated sample covariance. In the same way, R~ is computed, but using the sample covariance of the KF innovation sequence. In the correlation method, the noise covariance matrices are estimated based on the sample autocorrelations between the innovations by exploiting the relations between the estimation error covariance and the innovation covariance [11, 16, 17]. The drawback of the covariance matching and the correlation methods is that they do not guarantee the positive definiteness of the calculated covariances and require a large window of data. Lastly, implementation of the maximum likelihood method is very complicated which makes it impractical.As noted the AKF procedures estimate not the ONC, but the TNC. Therefore the estimation of the ONC remains open. Knowledge of the ONCs is necessary to study the quality of a system model. The resulting ONCs help to understand how much environmental disturbance affects the behavior of dynamics and to evaluate the accuracy of sensor measurements. However, if the processG or measurement D (“noise gain matrix”) is rectangular, we then have both theoretical and practical challenges in calculating the ONC. As will be shown shortly, popular approaches including the usage of the standard least squares approach and current data may lead to unsatisfactory results. To the best of the authors’ knowledge, there are no existing results for the noise covariance estimation within the Kalman filtering framework utilizing not only current data, but finite-memory information from the most recent time interval (receding horizon).The primary aim of this paper is to expand the previous results concerning the standard least squares covariance estimation to a new robust receding horizon covariance estimation. We propose a novel general estimation method for calculation of the ONCs with rectangular noise gain matricesG and D. The method combines the least squares (LS) optimization and receding horizon strategy. The proposed estimator can be integrated into the AKF for simultaneous estimation of both transformed Q~,R~ and original Q,R noise covariances.The main contributions of the paper are listed in the following:(1)
A novel receding horizon least squares optimization criterion is proposed. Using this criterion, an optimal covariance estimator for an arbitrary structure of noise gain matrices is derived.(2)
Constant noise gain matrices are comprehensively investigated, including the derivation of a recursive procedure for the optimal covariance estimator with a horizon interval and a stopping condition during the iteration process.(3)
Performance of the proposed estimator is illustrated on theoretical and practical examples for proving its correctness and efficiency.This paper is organized as follows. In Section2, two motivation examples within the Kalman filtering framework are demonstrated. Section 3 presents the standard LS noise covariance estimator and describes consequences of processing current (single) data only. In Section 4, a novel receding horizon least squares estimator (RHLSE) is proposed. It is based on the generalized LS approach and the novel receding horizon LS criterion. Here, the cardinal equations for determining the RHLSE and its equivalent modifications are derived and discussed. In Section 5, the algorithm of searching for the best horizon length improving the accuracy of the RHLSE is proposed. Section 6 demonstrates the practical usage of the RHLSE in adaptive Kalman filtering. In Section 7, performance and effectiveness of the proposed estimator are studied via a theoretical example with rectangular noise gain matrices and a motion with random velocity. Finally, we conclude the paper in Section 8. The list of main notations is given in Table 1.Table 1
List of main notations.
R n Set ofn-dimensional real column-vectors R n × m Set ofn×m real matrices O n × m n × m zero matrix A T Transpose of matrix I n Identity matrix of sizen×n A - 1 Inverse ofn×n matrix A A + Moore–Penrose pseudoinverse ofn×m matrix A t r ( A ) Trace ofn×n matrix A A , B Row block matrix consisting ofA and B A > 0   A ≥ 0 Positive definite (semidefinite) matrix N m , C Normal distribution with meanm and covariance matrix C w k ~ N 0 , Q Normal process with zero-mean and covariance matrixQ w k ~ 0 , Q Random process with zero-mean and covariance matrixQ E ∗ Expectation operator c o v x k Covariance (covariance matrix) of random vectorxk c o v x k , y k Cross-covariance ofxk and yk A F Frobenius norm (F-norm) of matrix,AF=trAAT A 1 1-norm of matrix,A1=maxj⁡∑i=1naij,A∈Rn×m A ∞ Infinity norm (I-norm) of matrix,A∞=maxi⁡∑j=1naij,A∈Rn×m A m a x Maximum norm (M-norm),Amax=maxi,j⁡aij A 2 Spectral norm (S-norm),A2=maxλATA


## 2. Motivation Examples with Rectangular Noise Gain Matrices

To begin, two typical motivation examples of a dynamical system with rectangular noise gain matrices are presented.Let us describe the noise estimation problem of present interest via two motivation examples. Suppose we have a linear discrete-time system,(2)xk+1=Fkxk+Gkwk,wk~0,Q,yk=Hkxk+Dkvk,vk~0,R.Then the TNCs Q~k=GkQGkT and R~k=DkRDkT represent covariance matrices of the linearly transformed noises w~k=Gkwk and v~k=Dkvk with the rectangular gain matrices Gk and Dk, respectively, i.e., Q~k=covw~k and R~k=covv~k.Motivation 1 (“colored process noise”).
Consider the scalar linear system with colored process noise:(3)xk+1=θxk+ηk,yk=xk+vk,where vk~N0,r is a measurement white noise.
Suppose thatηk is a colored process noise generated by the autoregressive equation of order p, ARp,(4)ηk=a1ηk-1+a2ηk-2…+apηk-p+wk,where wk~N0,q is a white noise.
Equation (4) may be considered as some system shaping the colored noise ηk from a white noise wk. Then the following is a state-space representation of the ARp model (4):(5)εk+1=F~εk+G~wk,where(6)εk=ηk-1ηk-2⋮ηk-p∈Rp,F~=a1a2⋯ap-1ap10⋯0001⋯00⋮⋮⋱⋮⋮0⋯010∈Rp×p,G~=10⋮0∈Rp.We combine (3) and (6) to obtain a system model with p+1–dimensional augmented state Xk and a scalar process noise wk~N0,q,(7)Xk+1=FkXk+Gkwk,yk=hkTXk+vk,where(8)Xk=xkεk∈Rp+1,Fk=θG~T0p×1F~∈Rp+1×p+1,Gk=0G~∈Rp+1,hk=10p×1∈Rp+1,wk,vk∈R.Motivation 2 (“composite process white noise and colored measurement noise”).
Consider the following discrete state-space model,(9a)xk+1=Fxk+G1ξ1,k+G2ξ2,k,(9b)yk=Cxk+ηk,(9c)ηk=Bηk-1+ek-1,where xk∈Rn, yk∈Rm, and ξ1,k∈Rr1, ξ2,k∈Rr2 are two uncorrelated process white noises with covariances Q1 and Q2, respectively; and ek∈Rm,ek~N0,R.
The state equation (9a) can be written in standard form with one composite process white noise wk∈Rr, and rectangular gain matrix G∈Rn×r,(10)xk+1=Fxk+Gwk,where G=G1,G2, wkT=ξ1,kT,ξ2,kT, wk∈Rr, r=r1+r2.
Following [2] we define an auxiliary measurement zk=yk+1-Byk described by the new measurement equation with white noise vk,(11)zk=Hxk+Dvk,where H=CF-BC, D=CG,Im, vkT=wkT,ekT,vk∈Rr+m.
Combining (10) and (11) we obtain the standard structure of the system model (2) with rectangular noise gain matrices G∈Rn×r and D∈Rm×(r+m).


## 3. Least Squares (LS) Estimation of Noise Covariance Using Single Data

As noted above, the purpose of this paper is to calculate the unknown ONCsQ and R based on the available TNCs Q~k and R~k obtained by the AKF algorithm. In this section, we present a simple LS solution to this problem.For a clear presentation, let us first rewrite (1) in the following form:(12)Akn×pXkp×pAkTp×n=Bkn×n,where(13a)Xk=XkT∈Rp×pis  an  unknown  ONC;(13b)Bk=BkT∈Rn×nis  a  known  TNC  (single  data);(13c)Ak∈Rn×pis  a  noise  gain  matrix,Ak=GkorDk.In a special case where the gain matrix Ak∈Rn×n is square and nonsingular, the LS solution to (12) takes the form Xk=Ak-1BkAk-T, where the known matrix Bk represents a current TNC at time instant k. But in practice, the gain matrix Ak∈Rn×p is rectangular, as the two motivation examples illustrated. In this rectangular case, the LS solution to (12) can be obtained using a pseudoinverse matrix.The linear matrix equation (12) has been considered by many authors [18–21]. They have studied general solutions to (12) for special solution structures, e.g., symmetric, triangular, or diagonal X using matrix decomposition techniques such as the singular value decomposition (SVD), the generalized SVD, and the canonical correlation decomposition.Theorem 3 (see [18]).
The least squares solution to matrix equation (12) takes the form(14)XkLS≜XkLSBk=Ak+BkAk+T.Note that the Moore–Penrose pseudoinverseAk+ always exists and is unique.However, the known TNCBk in practice is usually obtained from an experiment. It is calculated by using the AKF. Since the LS solution XkLSBk uses only the single (approximately calculated) TNC Bk, it is likely to be unstable or inaccurate. Therefore to overcome this disadvantage we propose a more realistic idea to use the set of data (TNCs), Bk,Bk-1,…,Bk-l+1, for the calculation of unknown covariance X. This set of data is obtained over the most recent time interval (receding horizon), k-l+1,k, saving in a finite memory. Using these data(15)Bk-l+1k=Bk-l+1,Bk-l+2,…,Bk,we can improve the LS solution (14), and, as a result, one can expect that the obtained receding horizon LS solution XkLSBk-l+1k will be more robust against processing and measurement errors than the LS solution XkLSBk which utilizes only the current data Bk.


## 4. Receding Horizon Least Squares Estimator (RHLSE)

In this section, we propose the novelreceding horizon least squares estimatorreferred to asRHLSE.In literature, the estimation problem in dynamical systems having various uncertainties has been solved in different ways, including finite-memory (receding horizon or sliding window) estimations [7, 22]. As a result, the finite-memory estimators are known to be more robust against model uncertainties and numerical errors than the standard LS estimator or Kalman filter [23–25]. For this reason, this receding horizon strategy is chosen to improve the LS solution (14).According to the strategy proposed above, the unknown solutionXk minimizes a new receding horizon LS criterion representing the sum of squared residuals over the receding horizon t∈k-l+1,k, i.e.,(16)minXkJk,lXkwith(17)Jk,lXk≜AkXkAkT-BkF2+Ak-1XkAk-1T-Bk-1F2+⋯+Ak-l+1XkAk-l+1T-Bk-l+1F2=∑t=k-l+1kAtXkAtT-BtF2where l≥1 is the receding horizon length (HL).Theorem 4 (RHLSE for time-varying noise gain matrices).
The solutionXk≜Xk,lRHLSE of the receding horizon least squares optimization problem represented by (17) satisfies the equation(18)∑t=k-l+1kMtXk,lRHLSEMt=∑t=k-l+1kAtTBtAt,where Mt=AtTAt.The proof of Theorem4 is given in Appendix.Next consider the important case in which the noise gain matrixAk≡A is time-invariant. Then the receding horizon LS problem (17) is rewritten as(19)minXkJk,lXkwith(20)Jk,lXk=∑t=k-l+1kAXkAT-BtF2,and we get the following result.Theorem 5 (RHLSE for time-invariant noise gain matrices).
The optimal receding horizon least squares optimization problem for constant noise gain matrices represented by (20) has the explicit solution,(21)Xk,lRHLSEBk-l+1k=A+B-k,lA+T,where B-k,l is the horizon average covariance,(22)B-k,l=1l∑t=k-l+1kBt.The proof of Theorem5 is given in Appendix.As we see from (21) and (22), the explicit solution Xk,lRHLSE depends only on the single data representing the horizon average covariance B-k,l. Therefore we get the following.Corollary 6.
The RHLSE (21) represents the LS estimator (14) calculated for the horizon average covariance, B-k,l, i.e., (23)Xk,lRHLSEBk-l+1k=XkLSB-k,l=A+B-k,lA+T.


## 5. Selection of Best Horizon Length for Time-Invariant Noise Gain Matrices

The optimal noise covariance (21) is sequentially calculated for each HL (iteration), l=1,2,…,l∗. Here we propose a stopping criterion for selecting the best HL, l∗. This stopping criterion is based on the errors between iterations. Let us first represent the solutions (21) and (22) in recursive form.The recursive formula for the horizon average covariance (22) takes the form(24)B-k,l=l-1lB-k,l-1+1lBk-l+1,⍈l=1,2,…,lmax.Using (21) and (24) the explicit solution Xk,lRHLSE is also rewritten in recursive form,(25)Xk,lRHLSE=A+B-k,lA+T=A+l-1lB-k,l-1+1lBk-l+1A+T=l-1lXk,l-1RHLSE+1lXk-l+1LS,⍈l=1,2,…,lmax,where(26)Xk-l+1LS=A+Bk-l+1A+T.For the stopping criterion based on the norm of iteration error(27)ΔXl,l-1≜Xk,lRHLSE-Xk,l-1RHLSE,⍈l≥2,we have the following result.Theorem 7.
LetXk,lRHLSE be the RHLSE (21). Then the Frobenius norm of the iteration error ΔXl,l-1 is calculated as(28)ΔXl,l-1F=1ltrM~2ΔBk,l-12,⍈l≥2,M~=A+TA+,ΔBk,l-1=B-k,l-1-Bk-l+1.The proof of Theorem7 is given in Appendix.Using (28), the best HL l∗ is defined as(29)l∗=min⁡l:ΔXl,l-1F<ε,⍈l≥2,where ε is a given tolerance.Remark 8.
As noted in Section3, the initial covariance Xk,1RHLSE=XkLSBk=A+BkA+T, calculated at l=1, may be inaccurate. Numerical examples given below confirm this fact. Therefore to improve the covariance Xk,1RHLSE, we need to continue the iterative process Xk,lRHLSE,l≥2, using the stopping condition (29).Remark 9.
In general, the stopping condition (29) does not guarantee the convergence of the iterative process (25); that is, the inequality X-Xk,lRHLSEF<ρX-Xk,l-1RHLSEF,ρ∈0,1 may not hold. Besides, the number of iterations is limited, l≤lmax. However, (29) seems to be still useful for stopping the iterative process, as will be demonstrated by examples.Remark 10.
The stopping condition (29) can be carried out in real time. In fact, the term trM~2ΔBk,l-12 depends on the TNCs (B-k,l-1 and Bk-l+1) for a given horizon. Therefore, the complexity of the real-time implementation of the stopping condition is not an issue.


## 6. Application of RHLSE

In this section, an application of the RHLSE to the estimation of process and measurement noise covariances is considered.


### 6.1. Kalman Filtering

The Kalman filtering framework involves estimation of the state of a discrete-time linear dynamical system with additive white Gaussian noise,(30)xk+1=Fkxk+Gkwk,yk=Hkxk+Dkvk,where state xk∈Rnx, measurement yk∈Rny, process noise wk∈Rnw, wk~N0,Q, and measurement noise vk∈Rnv, vk~N0,R. All matrices Fk,Gk,Hk,Dk,Q≥0 and R>0 are with appropriate dimensions. We assume that the initial state x0~Nx-0,C0, and process and measurement noises wk,vk are mutually uncorrelated with constant covariances (ONCs), Q∈Rnw×nw and R∈Rnv×nv, respectively.As mentioned in Introduction the KF equations depend on the known TNCsQ~k=GkQGkT and R~k=DkRDkT, rather than the ONCs Q,R. We have(31)Prediction  step:x^k-=Fk-1x^k-1,x^0=x-0,Pk-=Fk-1Pk-1Fk-1T+Q~k-1,P0=C0,Correction  step:Kk=Pk-HkTHkPk-HkT+R~k-1,x^k=x^k-+Kkyk-Hkx^k-,Pk=Inx-KkHkPk-.


### 6.2. Adaptive Kalman Filtering

In a real problem the ONCs and, as a consequence, the TNCs are unknown. In this case the AKF simultaneously estimates both the statex^k and TNCs Q~,R~. The general scheme for adaptive estimation of the TNCs can be presented in Table 2.Table 2
Adaptive estimation of transformed noise covariances.
R estimation Q estimation Sequence Innovation sequence Residual sequence y ~ k = y k - H k x ^ k - x ~ k = x ^ k - x ^ k -  Sample covariance C y ~ k = 1 k 0  ∑ t = k - k 0 + 1 k y ~ t y ~ t T C x ~ k = 1 k 0  ∑ t = k - k 0 + 1 k x ~ t x ~ t T  Estimated noise covariance (TNC) R ~ k = C y ~ k - H k P k - H k T Q ~ k = C x ~ kThen the AKF equations for the system model (30) are given by(32)Initial  values:x^0=x-0,P0=C0,Q0>0,R0>0.Q~0=G0Q0G0T,R~0=D0R0D0T.x^k-=Fk-1x^k-1,Pk-=Fk-1Pk-1Fk-1T+Q~k-1,y~k=yk-Hkx^k-,Cy~k=1k0∑t=k-k0+1ky~ty~tT,R~k=Cy~k-HkPk-HkT,Kk=Pk-HkTHkPk-HkT+R~k-1,x^k=x^k-+Kkyk-Hkx^k-,Pk=Inx-KkHkPk-,x~k=x^k-x^k-,Q~k=1n0∑t=k-n0+1kx~tx~tT.To estimate the unknown ONCs Q,R using the available TNCs Q~k,R~k we apply the RHLSE. Equation (21) for the corresponding estimates R^k,l≜Rk,lRHLSE and Q^k,l≜Qk,lRHLSE are(33)R^k,l=D+B-k,lrD+T,B-k,lr=1l∑t=k-l+1kR~t,Q^k,l=G+B-k,lqG+T,B-k,lq=1l∑t=k-l+1kQ~t.Here Q~t and R~t are calculated in (32), and the best HLs l∗=lr for R^k,l and l∗=lq for Q^k,l are determined by (29), i.e.,(34)lr=min⁡l:ΔRl,l-1F<ε,l≥2,ΔRl,l-1F=1ltrM~D2ΔR~k,l-12,M~D=D+TD+,ΔR~k,l-1=B-k,l-1r-R~k-l+1,and(35)lq=min⁡l:ΔQl,l-1F<ε,l≥2,ΔQl,l-1F=1ltrM~G2ΔQ~k,l-12,M~G=G+TG+,ΔQ~k,l-1=B-k,l-1q-Q~k-l+1,respectively.


## 6.1. Kalman Filtering

The Kalman filtering framework involves estimation of the state of a discrete-time linear dynamical system with additive white Gaussian noise,(30)xk+1=Fkxk+Gkwk,yk=Hkxk+Dkvk,where state xk∈Rnx, measurement yk∈Rny, process noise wk∈Rnw, wk~N0,Q, and measurement noise vk∈Rnv, vk~N0,R. All matrices Fk,Gk,Hk,Dk,Q≥0 and R>0 are with appropriate dimensions. We assume that the initial state x0~Nx-0,C0, and process and measurement noises wk,vk are mutually uncorrelated with constant covariances (ONCs), Q∈Rnw×nw and R∈Rnv×nv, respectively.As mentioned in Introduction the KF equations depend on the known TNCsQ~k=GkQGkT and R~k=DkRDkT, rather than the ONCs Q,R. We have(31)Prediction  step:x^k-=Fk-1x^k-1,x^0=x-0,Pk-=Fk-1Pk-1Fk-1T+Q~k-1,P0=C0,Correction  step:Kk=Pk-HkTHkPk-HkT+R~k-1,x^k=x^k-+Kkyk-Hkx^k-,Pk=Inx-KkHkPk-.


## 6.2. Adaptive Kalman Filtering

In a real problem the ONCs and, as a consequence, the TNCs are unknown. In this case the AKF simultaneously estimates both the statex^k and TNCs Q~,R~. The general scheme for adaptive estimation of the TNCs can be presented in Table 2.Table 2
Adaptive estimation of transformed noise covariances.
R estimation Q estimation Sequence Innovation sequence Residual sequence y ~ k = y k - H k x ^ k - x ~ k = x ^ k - x ^ k -  Sample covariance C y ~ k = 1 k 0  ∑ t = k - k 0 + 1 k y ~ t y ~ t T C x ~ k = 1 k 0  ∑ t = k - k 0 + 1 k x ~ t x ~ t T  Estimated noise covariance (TNC) R ~ k = C y ~ k - H k P k - H k T Q ~ k = C x ~ kThen the AKF equations for the system model (30) are given by(32)Initial  values:x^0=x-0,P0=C0,Q0>0,R0>0.Q~0=G0Q0G0T,R~0=D0R0D0T.x^k-=Fk-1x^k-1,Pk-=Fk-1Pk-1Fk-1T+Q~k-1,y~k=yk-Hkx^k-,Cy~k=1k0∑t=k-k0+1ky~ty~tT,R~k=Cy~k-HkPk-HkT,Kk=Pk-HkTHkPk-HkT+R~k-1,x^k=x^k-+Kkyk-Hkx^k-,Pk=Inx-KkHkPk-,x~k=x^k-x^k-,Q~k=1n0∑t=k-n0+1kx~tx~tT.To estimate the unknown ONCs Q,R using the available TNCs Q~k,R~k we apply the RHLSE. Equation (21) for the corresponding estimates R^k,l≜Rk,lRHLSE and Q^k,l≜Qk,lRHLSE are(33)R^k,l=D+B-k,lrD+T,B-k,lr=1l∑t=k-l+1kR~t,Q^k,l=G+B-k,lqG+T,B-k,lq=1l∑t=k-l+1kQ~t.Here Q~t and R~t are calculated in (32), and the best HLs l∗=lr for R^k,l and l∗=lq for Q^k,l are determined by (29), i.e.,(34)lr=min⁡l:ΔRl,l-1F<ε,l≥2,ΔRl,l-1F=1ltrM~D2ΔR~k,l-12,M~D=D+TD+,ΔR~k,l-1=B-k,l-1r-R~k-l+1,and(35)lq=min⁡l:ΔQl,l-1F<ε,l≥2,ΔQl,l-1F=1ltrM~G2ΔQ~k,l-12,M~G=G+TG+,ΔQ~k,l-1=B-k,l-1q-Q~k-l+1,respectively.


## 7. Simulation Examples

In this section, we first test the efficiency of the proposed RHLSE and the stopping criterion on a theoretical example for which the true ONC is known (Section7.1). And then, performance of the RHLSEs (33)-(35) shall be demonstrated with the example of a moving object (Section 7.2).


### 7.1. Theoretical Example

Let us estimate a random constantθ~Nmθ,σθ2 given two single sensor measurements y1,k and y2,k of θ corrupted by their own uncorrelated Gaussian white noises (measurement errors) v1,k and v2,k, respectively. We assume that the two sensors contain a common environmental Gaussian white noise ηk, which is uncorrelated with v1,k and v2,k. Then the dynamic equations describing this situation are(36)State:xk+1=xk,k=0,1,…;x0=θ~Nmθ,σθ2,Sensor  1:y1,k=xk+v1,k+ηk,v1,k~N0,r1,ηk~N0,ση2,Sensor  2:y2,k=xk+v2,k+ηk,v1,k~N0,r2.Here we calculate the estimation accuracy of the proposed RHLSE in terms of the receding HL, l≥1, for the two measurement models.Model 1 (“measurement model without environmental noiseηk”).
(37) y 1 , k = x k + v 1 , k , y 2 , k = x k + v 2 , k ,or in matrix form(38)Yk=Hxk+D1Vk1,Yk=y1,ky2,k,H=11,D1=1001,Vk1=v1,kv2,k~N0,R1.Here the ONC and TNC take the form(39)R1=r100r2,R~1=D1R1D1T=r100r2,respectively.
For the model (36)-(39) the measurement noise gain matrix D≜D1=I2, and D+≜D1+=I2 and M~D=D+TD+=I2. Then using (33) and (34) the RHLSE R^k,l≜Rk,lRHLSE for the measurement noise covariance X≜R(1) and the best HL l∗≜lr take the following form,(40)R^k,l=B-k,lr=1l∑t=k-l+1kR~t1,l=2,3,…,l∗,l∗=min⁡l:ΔR~l,l-1F<ε,where(41)ΔR~l,l-1F=1ltrΔR~k,l-12,l≥2,ΔR~k,l-1=B-k,l-1r-R~k-l+11,R~t1=Cy~t-HPt-HT,Cy~t=1t-t0∑i=t0+1tY~iY~iT,t=k-l+1,,…,k,Y~i=Yi-Hx^i-1,i=t0+1,…,t,and t0 is chosen empirically to give some statistical smoothing.The simulation results are illustrated in Table3 and Figure 1 for k=50, r1=0.1, r2=0.8, and lmax=30. Our point of interest are the values of different norms of the true error Eltrue≜X-Xk,lRHLSE=R(1)-R^k,l, and the best HL l∗. In Table 3, a sign ∗ represents one of the five considered matrix norms, ∗F, ∗M, ∗∞, ∗1, and ∗2 (Note: ∗∞=∗1 for symmetric matrices). The columns II and III show the absolute norm (A-norm) and the relative norm (R-norm) of the iteration error ΔRl,l-1=R^k,l-R^k,l-1, respectively.Table 3
Norms of true and iteration errors for Model1.
HLl Absolute and relative norms for iteration errorΔRl,l-1 Norm of true error∗=Eltrue A-normΔRl,l-1F R-normΔRl,l-1FR^k,lF F-norm∗F M-norm∗max I-norm∗∞=∗1 S-norm∗2 I II III IV V VI VII 1 - - - - - - 0.3832 0.3365 0.3572 0.3392 2 0.0276 0.0417 0.3674 0.3308 0.3398 0.3312 3 0.0207 0.0316 0.3717 0.3330 0.3381 0.3331 4 0.0122 0.0185 0.3712 0.3341 0.3472 0.3351 5 0.0114 0.0171 0.3693 0.3325 0.3534 0.3350 6 0.0085 0.0129 0.3680 0.3279 0.3520 0.3314 7 0.0220 0.0330 0.3716 0.3330 0.3713 0.3410 8 0.0194 0.0287 0.3793 0.3434 0.3912 0.3544 9 0.0180 0.0263 0.3852 0.3515 0.4080 0.3654 10 0.0222 0.0318 0.3887 0.3565 0.4245 0.3748 11 0.0114 0.0162 0.3933 0.3637 0.4331 0.3816 12 0.0084 0.0118 0.3974 0.3700 0.4391 0.3871 13 0.0076 0.0106 0.4008 0.3740 0.4459 0.3919 14 0.0137 0.0188 0.4017 0.3778 0.4511 0.3953 15 0.0110 0.0149 0.4033 0.3817 0.4553 0.3986 16 0.0090 0.0121 0.4063 0.3868 0.4590 0.4025 17 0.0046 0.0062 0.4103 0.3912 0.4625 0.4064 18 0.0103 0.0138 0.4120 0.3951 0.4647 0.4091 19 0.0062 0.0083 0.4135 0.3982 0.4642 0.4107 20 0.0073 0.0097 0.4155 0.4017 0.4631 0.4124 21 0.0102 0.0138 0.4188 0.4045 0.4628 0.4144 22 0.0099 0.0135 0.4209 0.4064 0.4608 0.4151 23 0.0103 0.0142 0.4223 0.4080 0.4567 0.4152 24 0.0125 0.0175 0.4226 0.4086 0.4496 0.4138 25 0.0117 0.0165 0.4230 0.4090 0.4428 0.4127 26 0.0096 0.0136 0.4236 0.4095 0.4372 0.4120 27 0.0068 0.0097 0.4238 0.4097 0.4328 0.4114 28 0.0143 0.0207 0.4264 0.4092 0.4284 0.4105 29 0.0188 0.0280 0.4307 0.4082 0.4251 0.4093 30 0.0180 0.0275 0.4353 0.4071 0.4217 0.4079Figure 1
Norms of true and iteration errors for Model1.Comments on Table3 and Figure 1 are in order.(1a)
The values of norms in Table3 confirm the known inequalities for matrix norms A of the symmetric matrix A=Eltrue∈Rn×n, i.e.,(42)Amax≤A2≤AF≤nA2,A2≤A∞.(1b)
The best HL for each of the presented norms is the following,lA∗=lR∗=lF∗=lM∗=6; lI∗=3; and lS∗=2or6. So a small HL, 2≤l≤6, is suitable for Model 1 in which the sizes of the ONC and TNC are equal, R(1),R~(1)∈R2×2.(1c)
The values of all norms of the true errorEltrue confirm Remark 8 (see columns IV~VII). They show that the initial HL l=1 is not optimal, i.e., the single ONC R^k,1 is not good.(1d)
The curves of all norms of the true error have a similar shape. The curves slightly fluctuate and then increase around a small HL,l<8.(1e)
The curves of the A- and R-norms are similar. Therefore, only one of them is shown in Figure1.(1f)
As noted in Remark9, fulfillment of the condition ΔRl,l-1F<ε does not guarantee the convergence. In Table 3 we can observe that the small iteration errors ΔR19,18F<ΔR20,19F<ΔR6,5F=0.0085<ε=0.01, correspond to the large true errors, E20trueF>E19trueF>E6trueF. This fact is also illustrated in Figure 2, where F-norm EltrueF and scaled A-norm 10ΔXl,l-1A are plotted for a clear comparison.Figure 2
Selection of best HL for Model1: comparison of true error (F-norm) and absolute iteration error (A-norm).Next, we consider more complicated Model2, in which the noise gain matrix D(2) is rectangular, and the sizes of the original R(2) and transformed R~(2) measurement covariances are different. In this case (12) represents an undefined system of equations.Model 2 (“measurement model with environmental noiseηk”).
(43) y 1 , k = x k + v 1 , k + η k , y 2 , k = x k + v 2 , k + η k ,or in matrix form(44)Yk=Hxk+D2Vk2,Yk=y1,ky2,k,H=11,D2=101011,Vk2=v1,kv2,kηk~N0,R2.Here the ONC and TNC take the form(45)R2=r1000r2000ση2,R~2=D2R2D2T=r1+ση2ση2ση2r2+ση2,respectively.
For the model (43)-(45),(46)D≜D2=101011,D+≜D2+=132-1-1211,and the RHLSE (33) takes the form R^k,l=D+B-k,lrD+T.
Setting upk=50, r1=0.1, r2=0.8, and ση2=0.3, we calculate the norms of iteration error ΔRl,l-1 and true error Eltrue. The simulation results are illustrated in Table 4 and Figure 3.Table 4
Norms of true and iteration errors for Model2.
HLl Absolute and relative norms for iteration errorΔRl,l-1 Norm of true error∗=Eltrue A-normΔRl,l-1F R-normΔRl,l-1FR^k,lF F-norm∗F M-norm∗max I-norm∗∞=∗1 S-norm∗2 I II III V VI VII VIII 1 - - - - - - 0.1813 0.1151 0.2303 0.1724 2 0.0858 0.1636 0.1230 0.0651 0.1303 0.1197 3 0.0726 0.1248 0.0752 0.0478 0.0957 0.0727 4 0.0582 0.0926 0.0549 0.0317 0.0635 0.0541 5 0.0105 0.0165 0.0532 0.0277 0.0555 0.0501 6 0.0095 0.0147 0.0543 0.0272 0.0544 0.0476 7 0.0146 0.0222 0.0608 0.0280 0.0561 0.0457 8 0.0237 0.0359 0.0713 0.0317 0.0634 0.0520 9 0.0089 0.0136 0.0786 0.0364 0.0728 0.0582 10 0.0207 0.0319 0.0894 0.0475 0.0950 0.0726 11 0.0110 0.0171 0.0990 0.0544 0.1088 0.0825 12 0.0115 0.0181 0.1024 0.0601 0.1202 0.0908 13 0.0147 0.0236 0.1100 0.0678 0.1357 0.1020 14 0.0162 0.0267 0.1179 0.0753 0.1507 0.1130 15 0.0122 0.0203 0.1241 0.0806 0.1613 0.1210 16 0.0079 0.0133 0.1308 0.0852 0.1705 0.1280 17 0.0084 0.0143 0.1319 0.0868 0.1736 0.1305 18 0.0112 0.0195 0.1335 0.0884 0.1768 0.1332 19 0.0143 0.0255 0.1361 0.0896 0.1793 0.1361 20 0.0125 0.0227 0.1378 0.0897 0.1795 0.1376 21 0.0092 0.0167 0.1395 0.0900 0.1801 0.1395 22 0.0099 0.0176 0.1360 0.0878 0.1757 0.1358 23 0.0103 0.0181 0.1320 0.0852 0.1705 0.1312 24 0.0109 0.0188 0.1269 0.0817 0.1635 0.1250 25 0.0115 0.0194 0.1241 0.0794 0.1589 0.1206 26 0.0166 0.0274 0.1207 0.0757 0.1515 0.1141 27 0.0122 0.0197 0.1187 0.0726 0.1453 0.1091 28 0.0150 0.0237 0.1165 0.0684 0.1369 0.1027 29 0.0123 0.0191 0.1181 0.0662 0.1325 0.0996 30 0.0206 0.0313 0.1216 0.0625 0.1250 0.0952Figure 3
Norms of true and iteration errors for Model2.In general, comments (1a)~(1f) made for Model1 are the same for Table 4 and Figure 3, but there are some minor changes concerning (1b) and (1d).(2b)
The best HLs for the norms are little different, i.e.,lA∗=lR∗=9, lF∗=5, lM∗=lI∗=6, lS∗=7. So a small HL, l≤9, is also suitable for Model 2, when the noise gain matrix D is rectangular.(2d)
The curves of all norms of the true error have a similar shape. The curves increase and then slightly decrease around a medium HL,9<l<20.Here is a summary of the simulation results presented for Models1 and 2.(i)
We propose to use the RHLSEXk,lRHLSE=A+B-k,lA+T with a small HL, 2≤l≤9, for both square and rectangular noise matrices A=GorD. Usage of short l=1 and long l≫1 horizon lengths may lead to unsatisfactory results.(ii)
In parallel with testing the efficiency of the stopping criterion we also study the behavior of the mean square errors (MSEs) for the covariance estimatesXk,lRHLSE∈R2×2 (Model 1) and Xk,lRHLSE∈R3×3 (Model 2),(47)MSEk,ij≜EXij-Xk,l;ijRHLSE2,i,j=1,…,p,whereXij is the true simulated value, Xk,l;ijRHLSE is the RHLSE estimate of the element Xij,k(=0,1,…) is the time instant, and l=l∗ is the best HL, X=Xij, Xk,lRHLSE=Xk,l;ijRHLSE. To calculate the MSEs the Monte-Carlo simulation with 1000 runs was used, i.e., the expectation operator in (47) is taken for the 1000 error data at each time instant k. We observe that the relative errors, ΔMSEk,ii(1), for Model 1,(48)ΔMSEk,ii1=MSEk,iiXii×100%,i=1,2,X11=r1,X22=r2,l∗=6,range from 10.8% to 13.7% within the time zonek∈0,50, and then they decrease. At k>500 the maximum values of the ΔMSEk,11(1) and ΔMSEk,22(1) reach the values of 7.4% and 8.8%, respectively. Similar results, but slightly worse due to the rectangularity of the noise gain matrix D(2), are obtained for Model 2. For this model the relative errors, (49)ΔMSEk,ii2=MSEk,iiXii×100%,i=1,2,3,X11=r1,X22=r2,X33=ση2,l∗=6,range from 12% to 16.5% within the time zonek∈0,50. At k>500 the maximum values of the ΔMSEk,11(2), ΔMSEk,22(2), and ΔMSEk,33(2) amount to 7.1%,8.4%, and 9.2%, respectively. Thus, this analysis of the MSEs shows that the RHLSE is suitable for practical applications.


### 7.2. Numerical Example: 1D Tracking with Random Velocity

In this section, the proposed RHLSE is used in adaptive filtering with unknown ONC (Q and R).Let the true target state at discrete timetk be defined as xk=pkυkT, where pk and υk denote the position and velocity in Cartesian coordinates, respectively. Assume that the constant velocity is subjected to a random disturbance wk. The measured position can be defined as yk=pk+ξk+ηk, where ξk is a sensory error and ηk is an environmental noise. The target state dynamics is then modeled as follows:(50)xk+1=Fxk+Gwk,xk∈R2,yk=Hxk+Dek,yk∈R,wherewk∈R,wk~N0,Q,ek∈R2,ek~N0,R,xk=x1,kx2,k=pkυk,F=1Δt01,G=01,ek=ξkηk,HT=10,DT=11,R=diagrξrη,x0~Nx-0,C0,Δt≪1.Simulation results with usage of the AKF (32) and RHLSE (12) are illustrated in Figures 4–8 for the following model parameters:(51)Δt=0.005,x-0=410T,C0=diag102.5,Q=0.25,lmax=100,k∈0;1000,rξ=0.64,rη=0.09.Figures 4 and 5 show the position and velocity estimates compared to the true values. The AKF appears to converge, in the sense that the estimation error tends to decrease, but slowly for velocity because the sensor measures only the position.Figure 4
True and estimated position,k∈[500,1000].Figure 5
True and estimated velocity,k∈[500,1000].Figure 6
Transformed measurement noise covariance (variance)R~k.Figure 7
Best HLlr=6 based on minimum of iteration error ΔR~l,l-1: comparison of the norms of true error and iteration error (scaled by 100 times).Figure 8
Best HLlq=5 based on minimum of iteration error ΔQ~l,l-1: comparison of the norms of true error and iteration error (scaled by 1000 times).The key value in adaptive filtering is the TNCR~k=Cy~k-HkPk-HkT. This covariance is plotted in Figure 6.As we observe, the covarianceR~k is not seriously changing for k≥500. For this reason we fixed the receding horizon interval k-l+1,k=500-l+1,500 for further estimation of the ONCs R and Q by formulas (33). We have(52)R^k,l=D+B-k,lrD+T,Q^k,l=G+B-k,lqG+T,D+=0.50.5,G+=01,k=500,l≥1.The best HLs lr and lq are calculated based on the iteration errors ΔR~l,l-1 and ΔQ~l,l-1, respectively,(53)lr=min⁡l:ΔR~l,l-1F<ε=10-5,ΔR~l,l-1F=1ltrM~D2ΔR~k,l-12,M~D=12,and(54)lq=min⁡l:ΔQ~l,l-1F<ε=2∙10-5,ΔQ~l,l-1F=1ltrM~G2ΔQ~k,l-12,M~G=0001.As indicated earlier in Remark 9, the best HLs lr and lq do not guarantee a minimum of the true errors ElR and ElQ, respectively. Therefore, in addition we calculate the norms of the true errors, ElRF and ElQF, which confirm that the proposed idea of choosing the best HL using the iterative error (53) or (54) leads to good results. Figure 7 presents the norms of the true error ElRF and the iteration error ΔR~l,l-1F for measurement covariance, respectively. Analogously, Figure 8 presents the norms of the true error ElQF and iteration error ΔQ~l,l-1F for process covariance, respectively. We observe that as well as in the theoretical example (Section 7) the obtained best HLs lr=6 and lq=5 are located within a small horizon interval. In Figure 8, both curves show one value for the best HL, lq=5. But in Figure 7 we observe five horizon lengths lr=6,46,52,5458, which satisfy condition ΔR~l,l-1F<10-5 within the horizon interval l∈1;100, and the last four lie in the middle of the horizon interval (see Table 5). Table 5 shows that the iteration errors, ΔR~l,l-1F, are different; however, the true errors ElRF remain almost the same. Hence, it can be concluded that the small length lr=6 is acceptable.Table 5
Best HLs for measurement noise covariance.
l r Δ R ~ l , l - 1 F E l R F 6 3.9237E-06 0.0070463 46 9.39816E-07 0.0070258 52 2.8749E-06 0.0069196 54 2.1627E-07 0.0070807 58 1.0062E-06 0.0070107To compare the relative MSEs,(55)ΔMSEkQ=MSEkQQ×100%,ΔMSEkξ=MSEkξrξ×100%,ΔMSEkη=MSEkηrη×100%,lq=5,lr=6,the Monte-Carlo simulation with 1000 runs was performed. The numerical values of all relative errors (55) within the time interval k∈0,300 range from 13.8% to 16.7%. And for large time instants, k>500, they do not change much. In this case, the maximum values of the errors ΔMSEk(Q), ΔMSEk(ξ), and ΔMSEk(η) amount to 14.5%, 12.7%, and 13%, respectively.


## 7.1. Theoretical Example

Let us estimate a random constantθ~Nmθ,σθ2 given two single sensor measurements y1,k and y2,k of θ corrupted by their own uncorrelated Gaussian white noises (measurement errors) v1,k and v2,k, respectively. We assume that the two sensors contain a common environmental Gaussian white noise ηk, which is uncorrelated with v1,k and v2,k. Then the dynamic equations describing this situation are(36)State:xk+1=xk,k=0,1,…;x0=θ~Nmθ,σθ2,Sensor  1:y1,k=xk+v1,k+ηk,v1,k~N0,r1,ηk~N0,ση2,Sensor  2:y2,k=xk+v2,k+ηk,v1,k~N0,r2.Here we calculate the estimation accuracy of the proposed RHLSE in terms of the receding HL, l≥1, for the two measurement models.Model 1 (“measurement model without environmental noiseηk”).
(37) y 1 , k = x k + v 1 , k , y 2 , k = x k + v 2 , k ,or in matrix form(38)Yk=Hxk+D1Vk1,Yk=y1,ky2,k,H=11,D1=1001,Vk1=v1,kv2,k~N0,R1.Here the ONC and TNC take the form(39)R1=r100r2,R~1=D1R1D1T=r100r2,respectively.
For the model (36)-(39) the measurement noise gain matrix D≜D1=I2, and D+≜D1+=I2 and M~D=D+TD+=I2. Then using (33) and (34) the RHLSE R^k,l≜Rk,lRHLSE for the measurement noise covariance X≜R(1) and the best HL l∗≜lr take the following form,(40)R^k,l=B-k,lr=1l∑t=k-l+1kR~t1,l=2,3,…,l∗,l∗=min⁡l:ΔR~l,l-1F<ε,where(41)ΔR~l,l-1F=1ltrΔR~k,l-12,l≥2,ΔR~k,l-1=B-k,l-1r-R~k-l+11,R~t1=Cy~t-HPt-HT,Cy~t=1t-t0∑i=t0+1tY~iY~iT,t=k-l+1,,…,k,Y~i=Yi-Hx^i-1,i=t0+1,…,t,and t0 is chosen empirically to give some statistical smoothing.The simulation results are illustrated in Table3 and Figure 1 for k=50, r1=0.1, r2=0.8, and lmax=30. Our point of interest are the values of different norms of the true error Eltrue≜X-Xk,lRHLSE=R(1)-R^k,l, and the best HL l∗. In Table 3, a sign ∗ represents one of the five considered matrix norms, ∗F, ∗M, ∗∞, ∗1, and ∗2 (Note: ∗∞=∗1 for symmetric matrices). The columns II and III show the absolute norm (A-norm) and the relative norm (R-norm) of the iteration error ΔRl,l-1=R^k,l-R^k,l-1, respectively.Table 3
Norms of true and iteration errors for Model1.
HLl Absolute and relative norms for iteration errorΔRl,l-1 Norm of true error∗=Eltrue A-normΔRl,l-1F R-normΔRl,l-1FR^k,lF F-norm∗F M-norm∗max I-norm∗∞=∗1 S-norm∗2 I II III IV V VI VII 1 - - - - - - 0.3832 0.3365 0.3572 0.3392 2 0.0276 0.0417 0.3674 0.3308 0.3398 0.3312 3 0.0207 0.0316 0.3717 0.3330 0.3381 0.3331 4 0.0122 0.0185 0.3712 0.3341 0.3472 0.3351 5 0.0114 0.0171 0.3693 0.3325 0.3534 0.3350 6 0.0085 0.0129 0.3680 0.3279 0.3520 0.3314 7 0.0220 0.0330 0.3716 0.3330 0.3713 0.3410 8 0.0194 0.0287 0.3793 0.3434 0.3912 0.3544 9 0.0180 0.0263 0.3852 0.3515 0.4080 0.3654 10 0.0222 0.0318 0.3887 0.3565 0.4245 0.3748 11 0.0114 0.0162 0.3933 0.3637 0.4331 0.3816 12 0.0084 0.0118 0.3974 0.3700 0.4391 0.3871 13 0.0076 0.0106 0.4008 0.3740 0.4459 0.3919 14 0.0137 0.0188 0.4017 0.3778 0.4511 0.3953 15 0.0110 0.0149 0.4033 0.3817 0.4553 0.3986 16 0.0090 0.0121 0.4063 0.3868 0.4590 0.4025 17 0.0046 0.0062 0.4103 0.3912 0.4625 0.4064 18 0.0103 0.0138 0.4120 0.3951 0.4647 0.4091 19 0.0062 0.0083 0.4135 0.3982 0.4642 0.4107 20 0.0073 0.0097 0.4155 0.4017 0.4631 0.4124 21 0.0102 0.0138 0.4188 0.4045 0.4628 0.4144 22 0.0099 0.0135 0.4209 0.4064 0.4608 0.4151 23 0.0103 0.0142 0.4223 0.4080 0.4567 0.4152 24 0.0125 0.0175 0.4226 0.4086 0.4496 0.4138 25 0.0117 0.0165 0.4230 0.4090 0.4428 0.4127 26 0.0096 0.0136 0.4236 0.4095 0.4372 0.4120 27 0.0068 0.0097 0.4238 0.4097 0.4328 0.4114 28 0.0143 0.0207 0.4264 0.4092 0.4284 0.4105 29 0.0188 0.0280 0.4307 0.4082 0.4251 0.4093 30 0.0180 0.0275 0.4353 0.4071 0.4217 0.4079Figure 1
Norms of true and iteration errors for Model1.Comments on Table3 and Figure 1 are in order.(1a)
The values of norms in Table3 confirm the known inequalities for matrix norms A of the symmetric matrix A=Eltrue∈Rn×n, i.e.,(42)Amax≤A2≤AF≤nA2,A2≤A∞.(1b)
The best HL for each of the presented norms is the following,lA∗=lR∗=lF∗=lM∗=6; lI∗=3; and lS∗=2or6. So a small HL, 2≤l≤6, is suitable for Model 1 in which the sizes of the ONC and TNC are equal, R(1),R~(1)∈R2×2.(1c)
The values of all norms of the true errorEltrue confirm Remark 8 (see columns IV~VII). They show that the initial HL l=1 is not optimal, i.e., the single ONC R^k,1 is not good.(1d)
The curves of all norms of the true error have a similar shape. The curves slightly fluctuate and then increase around a small HL,l<8.(1e)
The curves of the A- and R-norms are similar. Therefore, only one of them is shown in Figure1.(1f)
As noted in Remark9, fulfillment of the condition ΔRl,l-1F<ε does not guarantee the convergence. In Table 3 we can observe that the small iteration errors ΔR19,18F<ΔR20,19F<ΔR6,5F=0.0085<ε=0.01, correspond to the large true errors, E20trueF>E19trueF>E6trueF. This fact is also illustrated in Figure 2, where F-norm EltrueF and scaled A-norm 10ΔXl,l-1A are plotted for a clear comparison.Figure 2
Selection of best HL for Model1: comparison of true error (F-norm) and absolute iteration error (A-norm).Next, we consider more complicated Model2, in which the noise gain matrix D(2) is rectangular, and the sizes of the original R(2) and transformed R~(2) measurement covariances are different. In this case (12) represents an undefined system of equations.Model 2 (“measurement model with environmental noiseηk”).
(43) y 1 , k = x k + v 1 , k + η k , y 2 , k = x k + v 2 , k + η k ,or in matrix form(44)Yk=Hxk+D2Vk2,Yk=y1,ky2,k,H=11,D2=101011,Vk2=v1,kv2,kηk~N0,R2.Here the ONC and TNC take the form(45)R2=r1000r2000ση2,R~2=D2R2D2T=r1+ση2ση2ση2r2+ση2,respectively.
For the model (43)-(45),(46)D≜D2=101011,D+≜D2+=132-1-1211,and the RHLSE (33) takes the form R^k,l=D+B-k,lrD+T.
Setting upk=50, r1=0.1, r2=0.8, and ση2=0.3, we calculate the norms of iteration error ΔRl,l-1 and true error Eltrue. The simulation results are illustrated in Table 4 and Figure 3.Table 4
Norms of true and iteration errors for Model2.
HLl Absolute and relative norms for iteration errorΔRl,l-1 Norm of true error∗=Eltrue A-normΔRl,l-1F R-normΔRl,l-1FR^k,lF F-norm∗F M-norm∗max I-norm∗∞=∗1 S-norm∗2 I II III V VI VII VIII 1 - - - - - - 0.1813 0.1151 0.2303 0.1724 2 0.0858 0.1636 0.1230 0.0651 0.1303 0.1197 3 0.0726 0.1248 0.0752 0.0478 0.0957 0.0727 4 0.0582 0.0926 0.0549 0.0317 0.0635 0.0541 5 0.0105 0.0165 0.0532 0.0277 0.0555 0.0501 6 0.0095 0.0147 0.0543 0.0272 0.0544 0.0476 7 0.0146 0.0222 0.0608 0.0280 0.0561 0.0457 8 0.0237 0.0359 0.0713 0.0317 0.0634 0.0520 9 0.0089 0.0136 0.0786 0.0364 0.0728 0.0582 10 0.0207 0.0319 0.0894 0.0475 0.0950 0.0726 11 0.0110 0.0171 0.0990 0.0544 0.1088 0.0825 12 0.0115 0.0181 0.1024 0.0601 0.1202 0.0908 13 0.0147 0.0236 0.1100 0.0678 0.1357 0.1020 14 0.0162 0.0267 0.1179 0.0753 0.1507 0.1130 15 0.0122 0.0203 0.1241 0.0806 0.1613 0.1210 16 0.0079 0.0133 0.1308 0.0852 0.1705 0.1280 17 0.0084 0.0143 0.1319 0.0868 0.1736 0.1305 18 0.0112 0.0195 0.1335 0.0884 0.1768 0.1332 19 0.0143 0.0255 0.1361 0.0896 0.1793 0.1361 20 0.0125 0.0227 0.1378 0.0897 0.1795 0.1376 21 0.0092 0.0167 0.1395 0.0900 0.1801 0.1395 22 0.0099 0.0176 0.1360 0.0878 0.1757 0.1358 23 0.0103 0.0181 0.1320 0.0852 0.1705 0.1312 24 0.0109 0.0188 0.1269 0.0817 0.1635 0.1250 25 0.0115 0.0194 0.1241 0.0794 0.1589 0.1206 26 0.0166 0.0274 0.1207 0.0757 0.1515 0.1141 27 0.0122 0.0197 0.1187 0.0726 0.1453 0.1091 28 0.0150 0.0237 0.1165 0.0684 0.1369 0.1027 29 0.0123 0.0191 0.1181 0.0662 0.1325 0.0996 30 0.0206 0.0313 0.1216 0.0625 0.1250 0.0952Figure 3
Norms of true and iteration errors for Model2.In general, comments (1a)~(1f) made for Model1 are the same for Table 4 and Figure 3, but there are some minor changes concerning (1b) and (1d).(2b)
The best HLs for the norms are little different, i.e.,lA∗=lR∗=9, lF∗=5, lM∗=lI∗=6, lS∗=7. So a small HL, l≤9, is also suitable for Model 2, when the noise gain matrix D is rectangular.(2d)
The curves of all norms of the true error have a similar shape. The curves increase and then slightly decrease around a medium HL,9<l<20.Here is a summary of the simulation results presented for Models1 and 2.(i)
We propose to use the RHLSEXk,lRHLSE=A+B-k,lA+T with a small HL, 2≤l≤9, for both square and rectangular noise matrices A=GorD. Usage of short l=1 and long l≫1 horizon lengths may lead to unsatisfactory results.(ii)
In parallel with testing the efficiency of the stopping criterion we also study the behavior of the mean square errors (MSEs) for the covariance estimatesXk,lRHLSE∈R2×2 (Model 1) and Xk,lRHLSE∈R3×3 (Model 2),(47)MSEk,ij≜EXij-Xk,l;ijRHLSE2,i,j=1,…,p,whereXij is the true simulated value, Xk,l;ijRHLSE is the RHLSE estimate of the element Xij,k(=0,1,…) is the time instant, and l=l∗ is the best HL, X=Xij, Xk,lRHLSE=Xk,l;ijRHLSE. To calculate the MSEs the Monte-Carlo simulation with 1000 runs was used, i.e., the expectation operator in (47) is taken for the 1000 error data at each time instant k. We observe that the relative errors, ΔMSEk,ii(1), for Model 1,(48)ΔMSEk,ii1=MSEk,iiXii×100%,i=1,2,X11=r1,X22=r2,l∗=6,range from 10.8% to 13.7% within the time zonek∈0,50, and then they decrease. At k>500 the maximum values of the ΔMSEk,11(1) and ΔMSEk,22(1) reach the values of 7.4% and 8.8%, respectively. Similar results, but slightly worse due to the rectangularity of the noise gain matrix D(2), are obtained for Model 2. For this model the relative errors, (49)ΔMSEk,ii2=MSEk,iiXii×100%,i=1,2,3,X11=r1,X22=r2,X33=ση2,l∗=6,range from 12% to 16.5% within the time zonek∈0,50. At k>500 the maximum values of the ΔMSEk,11(2), ΔMSEk,22(2), and ΔMSEk,33(2) amount to 7.1%,8.4%, and 9.2%, respectively. Thus, this analysis of the MSEs shows that the RHLSE is suitable for practical applications.


## 7.2. Numerical Example: 1D Tracking with Random Velocity

In this section, the proposed RHLSE is used in adaptive filtering with unknown ONC (Q and R).Let the true target state at discrete timetk be defined as xk=pkυkT, where pk and υk denote the position and velocity in Cartesian coordinates, respectively. Assume that the constant velocity is subjected to a random disturbance wk. The measured position can be defined as yk=pk+ξk+ηk, where ξk is a sensory error and ηk is an environmental noise. The target state dynamics is then modeled as follows:(50)xk+1=Fxk+Gwk,xk∈R2,yk=Hxk+Dek,yk∈R,wherewk∈R,wk~N0,Q,ek∈R2,ek~N0,R,xk=x1,kx2,k=pkυk,F=1Δt01,G=01,ek=ξkηk,HT=10,DT=11,R=diagrξrη,x0~Nx-0,C0,Δt≪1.Simulation results with usage of the AKF (32) and RHLSE (12) are illustrated in Figures 4–8 for the following model parameters:(51)Δt=0.005,x-0=410T,C0=diag102.5,Q=0.25,lmax=100,k∈0;1000,rξ=0.64,rη=0.09.Figures 4 and 5 show the position and velocity estimates compared to the true values. The AKF appears to converge, in the sense that the estimation error tends to decrease, but slowly for velocity because the sensor measures only the position.Figure 4
True and estimated position,k∈[500,1000].Figure 5
True and estimated velocity,k∈[500,1000].Figure 6
Transformed measurement noise covariance (variance)R~k.Figure 7
Best HLlr=6 based on minimum of iteration error ΔR~l,l-1: comparison of the norms of true error and iteration error (scaled by 100 times).Figure 8
Best HLlq=5 based on minimum of iteration error ΔQ~l,l-1: comparison of the norms of true error and iteration error (scaled by 1000 times).The key value in adaptive filtering is the TNCR~k=Cy~k-HkPk-HkT. This covariance is plotted in Figure 6.As we observe, the covarianceR~k is not seriously changing for k≥500. For this reason we fixed the receding horizon interval k-l+1,k=500-l+1,500 for further estimation of the ONCs R and Q by formulas (33). We have(52)R^k,l=D+B-k,lrD+T,Q^k,l=G+B-k,lqG+T,D+=0.50.5,G+=01,k=500,l≥1.The best HLs lr and lq are calculated based on the iteration errors ΔR~l,l-1 and ΔQ~l,l-1, respectively,(53)lr=min⁡l:ΔR~l,l-1F<ε=10-5,ΔR~l,l-1F=1ltrM~D2ΔR~k,l-12,M~D=12,and(54)lq=min⁡l:ΔQ~l,l-1F<ε=2∙10-5,ΔQ~l,l-1F=1ltrM~G2ΔQ~k,l-12,M~G=0001.As indicated earlier in Remark 9, the best HLs lr and lq do not guarantee a minimum of the true errors ElR and ElQ, respectively. Therefore, in addition we calculate the norms of the true errors, ElRF and ElQF, which confirm that the proposed idea of choosing the best HL using the iterative error (53) or (54) leads to good results. Figure 7 presents the norms of the true error ElRF and the iteration error ΔR~l,l-1F for measurement covariance, respectively. Analogously, Figure 8 presents the norms of the true error ElQF and iteration error ΔQ~l,l-1F for process covariance, respectively. We observe that as well as in the theoretical example (Section 7) the obtained best HLs lr=6 and lq=5 are located within a small horizon interval. In Figure 8, both curves show one value for the best HL, lq=5. But in Figure 7 we observe five horizon lengths lr=6,46,52,5458, which satisfy condition ΔR~l,l-1F<10-5 within the horizon interval l∈1;100, and the last four lie in the middle of the horizon interval (see Table 5). Table 5 shows that the iteration errors, ΔR~l,l-1F, are different; however, the true errors ElRF remain almost the same. Hence, it can be concluded that the small length lr=6 is acceptable.Table 5
Best HLs for measurement noise covariance.
l r Δ R ~ l , l - 1 F E l R F 6 3.9237E-06 0.0070463 46 9.39816E-07 0.0070258 52 2.8749E-06 0.0069196 54 2.1627E-07 0.0070807 58 1.0062E-06 0.0070107To compare the relative MSEs,(55)ΔMSEkQ=MSEkQQ×100%,ΔMSEkξ=MSEkξrξ×100%,ΔMSEkη=MSEkηrη×100%,lq=5,lr=6,the Monte-Carlo simulation with 1000 runs was performed. The numerical values of all relative errors (55) within the time interval k∈0,300 range from 13.8% to 16.7%. And for large time instants, k>500, they do not change much. In this case, the maximum values of the errors ΔMSEk(Q), ΔMSEk(ξ), and ΔMSEk(η) amount to 14.5%, 12.7%, and 13%, respectively.


## 8. Conclusion

In some application problems, knowledge of noise covariances is necessary to study quality of the system and measurement model. In order to estimate arbitrary process and measurement noise covariances, a receding horizon estimation criterion with a generalized least squares approach is proposed (Theorem4).Special attention is given for time-invariant noise gain matrices (Theorem5). In this case the proposed RHLSE is comprehensively investigated, including the derivation of compact matrix form for recursive RHLSE (25) and an effective stopping condition for iterative estimation process (Theorem 7). The computational procedure based on the Moore–Penrose pseudoinverse for evaluation of the best noise covariances is recommended.In view of the importance of noise covariances in practice, the proposed RHLSE is illustrated on theoretical and practical dynamical models with different types of noises. The examples show that the RHLSE with selection of the best HL within a small horizon interval yields reasonably good estimation accuracy.Simulation analysis and comparison of the proposed estimator by several norms of estimation error show that usage of the receding horizon strategy is a very effective approach for the achievement of robustness against any uncertainties and errors.



---

*Source: 5303694-2018-11-19.xml*