Least-Mean-Square Receding Horizon Estimation

We propose a least-mean-square LMS receding horizon RH estimator for state estimation. The proposed LMS RH estimator is obtained from the conditional expectation of the estimated state given a finite number of inputs and outputs over the recent finite horizon. Any a priori state information is not required, and existing artificial constraints for easy derivation are not imposed. For a general stochastic discrete-time state space model with both system and measurement noise, the LMS RH estimator is explicitly represented in a closed form. For numerical reliability, the iterative form is presented with forward and backward computations. It is shown through a numerical example that the proposed LMS RH estimator has better robust performance than conventional Kalman estimators when uncertainties exist.


Introduction
Several criteria have been often employed for the design of optimal estimators.In particular, the mean-square-error criterion is the most popular and has many applications since it offers a simple closed-form solution as well as important geometric and physical interpretations.It is well known that the optimal estimators based on the mean-square-error criterion is obtained from the conditional expectation of the estimated variable given the known measurements.
For state estimation, many trials have been conducted to obtain a receding horizon RH estimator based on the mean-square-error criterion.At the current time k, the leastmean-square LMS RH estimator is to estimate the state x k−h at the time k −h from the inputs u k−• and the outputs y k−• over the recent finite horizon k − N, k − 1 , which can be written as 1.1

System
Receding horizon estimator Estimated state The structure of the receding horizon estimator D is a unit delay component .
where h and N, that is, the lag size and the memory size are design parameters to be determined, respectively.Available inputs and outputs are regarded as given conditions and the corresponding conditional expectation of the state at time k − h is obtained as its optimal estimate.Practically, inputs and outputs are known variables and hence can be considered as conditions.The structure of the LMS RH estimator 1.1 is depicted in Figure 1.As mentioned before, the LMS RH estimator 1.1 minimizes the mean-square-error criterion, E For h ≥ 2, h 1, and h ≤ 0, the estimators 1.1 are often called the smoothers, the filters, and the predictors, respectively."Receding horizon" is traced from the fact that the finite-time horizon, where inputs and outputs necessary for estimating unknown states are available, recedes with time.
In designing the controls, receding horizon schemes have already been popular 1, 2 .If N approaches ∞, the estimator 1.1 reduces to the well-known stationary Kalman estimator with infinite-memory.So, the LMS RH estimator 1.1 can be also called a finite memory estimator for the finite and fixed memory size N.It has been illustrated through numerical simulation and analysis that the LMS RH estimators with finite memory have better robust performance than conventional Kalman estimators with infinite memory 3, 4 .Also in input/output models arising in signal processing area, it is acknowledged that the finitememory or finite-impulse-response FIR filters have been preferable for practical reasons 5 .
In spite of the good performance and the usefulness of the LMS RH estimators, no one has proposed a general result on the conditional expectation 1.1 .Since it was difficult to obtain the conditional expectation 1.1 for a general state space model, some assumptions have been made to simplify the problem and then obtain a solution easily.At first, system or measurement noise were set to zero, which offers a closed-form solution easily 6-10 .In 11, 12 , a priori information on the initial state on the horizon was assumed to be known for obtaining a solution easily.Instead of directly obtaining a closed-form solution, the duality of a control and the complicated scattering theory from a physical phenomenon were used to show the feasibility of implementation for a general state space model 13, 14 .For easy derivation, the Kalman filter was also employed with infinite covariance 4 .However, this approach is so heuristic that the optimality is not guaranteed in the sense of the meansquare-error criterion.Besides, the system matrix is required to be nonsingular.As in other conventional estimators, there were also some trials that unbiased and linear constraints are imposed to obtain the optimal estimator 3, 15-17 .However, external control inputs are not considered 3 and the system matrix is required to be nonsingular 15 .Furthermore, it is not guaranteed that even though such constraints are removed, the optimality is still preserved.Though computing the conditional expectation 1.1 looks like a very simple problem, there is no result on a closed-form solution corresponding to 1.1 for a general state space model without any artificial assumptions and requirements.
In this paper, existing artificial assumptions for obtaining a solution easily are not made and any conditions are not required.Unlike the existing results, the system matrix is not required to be nonsingular.Both system and measurement noise are considered together with external control inputs.The LMS RH estimator will be directly obtained from the conditional expectation 1.1 , which automatically guarantees its optimality.It turns out that the proposed LMS RH estimator has the deadbeat property and the linear structure with inputs and outputs over the recent finite horizon.
The rest of this paper is organized as follows.In Section 2, the LMS RH estimator is obtained from the conditional expectation and its iterative computation is introduced in Section 3. A numerical simulation is carried out in Section 4 to illustrate the performance of the proposed LMS RH estimator.Finally, conclusions are presented in Section 5.

LMS RH Estimator
Consider a linear discrete-time state space model with an external control input: where x i ∈ n , u i ∈ l , and y i ∈ q are the state, the input, and the output, respectively, and the system noise w i ∈ p and the measurement noise v i ∈ q are assumed to be zero-mean white Gaussian with where nonzero S often happens when I/O models are converted to state space models.It is also assumed that A, C of the system is observable.Through this paper, the current time is denoted by k.For mathematical tractability, the state space model 2.1 is equivalently changed to where A s A − GSR −1 C and w s,i w i − SR −1 v i .It can be easily shown that w s,i and v i in 2.3 are not correlated.In other words, we have where Q − SR −1 S T ≥ 0. It is noted that off-diagonal blocks are filled with zeros.

Mathematical Problems in Engineering
For simple representations, several variables are defined as where ⊕ denotes the direct sum of the matrices.Additionally, W k−1 , V k−1 , and U k−1 are defined by replacing y in 2.5 with w s , v, and u, respectively.G N and S N are also defined by replacing B in B N with G and GSR −1 , respectively.First, we consider the case of 0 ≤ h ≤ N − 1 in 1.1 .A prediction problem for h ≤ −1 will be discussed later on.
By using the defined variables, the state x k−h to be estimated is represented in terms of the initial state x k−N on the horizon, inputs, and system noise on the recent finite horizon k − N, k − 1 as where L g,N is given by and L b,N and L s,N are obtained by replacing G in L g,N with B and GSR −1 , respectively.From 2.9 , the conditional expectation 1 can be represented as For this purpose, the system 2.3 is represented in a batch form on the recent finite horizon k − N, k − 1 as follows: 12 from which we obtain We can see from 2.13 that W k−1 and V k−1 are linearly, more correctly affinely, dependent on W k−1 and x k−N .The joint probability density function of W k−1 and x k−N , that is, for an appropriate scaling factor k. How to choose k will be discussed later on.By using the probability density function of W k−1 and V k−1 given as where a constant k is chosen to satisfy the following normalization condition:

Mathematical Problems in Engineering
In the same way, is an exponential function with the following exponent: where * denotes symmetric parts and W 1,1 , W 1,2 , and W 2,2 are given by

2.18
By completing the square of the terms of the exponent and recalling the integration of Gaussian functions from −∞ to ∞, we can compute To begin with, we introduce the following completion of squares: for some vectors a and b, and some matrices α, β, and γ of appropriate dimensions.The relation 2.20 can be easily obtained in a similar way to the mean of normal distribution.
From the following correspondences: 2.21 we have where W 1,1 , W 1,2 , and W 2,2 are given by 2.18 .It is noted that if A, C of the system 2.1 is observable, W 1,1 is positive definite and W 2,2 − W T 1,2 W −1 1,1 W 1,2 is also positive definite, which implies that the block matrix including W 1,1 , W 1,2 , and W 2,2 is guaranteed to be positive definite and hence nonsingular.Substituting 2.22 into 2.11 , we can obtain Through a long and tedious algebraic calculation, we have the covariance matrix of the mean-square-error T as follows: where Ξ N is defined by What we have done until now can be summarized in the following theorem.
Theorem 2.1.Suppose that A, C of the system 2.1 is observable.Then, the LMS RH estimator 1.1 is given by

2.25
where , and W 2,2 are defined in 2.18 .The corresponding covariance matrix is given as 2.23 .
It is noted that the proposed LMS RH estimator 2.25 is designed without requirements of the removal of some noise and the nonsingular system matrix A. In previous work, those requirements are adopted to solve the problems more easily.If constraints or assumptions of existing results are applied to the proposed result, the latter reduces to the former.For example, Q and B N can be set to zero for comparison with the results on systems without system noise and inputs, respectively.As seen in 2.25 , the LMS RH estimator is linear with respect to inputs and outputs on the recent finite horizon k − N, k − 1 .So, "finite memory" may be called "finite impulse response" that is often used in linear signal processing systems.In addition, the deadbeat property is guaranteed in the following theorem.Theorem 2.2.If no noise is applied, the LMS RH estimator 2.25 is a deadbeat estimator, that is, x k−h|k x k−h at the fixed-lag (h ≥ 1) or the current (h 0) time.
Proof.The LMS RH estimator 2.25 can be rearranged in terms of Y k−1 and U k−1 , that is, HY k−1 LU k−1 for some gain matrices H and L. Substituting Y k−1 in 2.12 into HY k−1 LU k−1 , using 2.9 , and removing all noise, we have

2.26
It can be easily seen that H C N A N−h s and terms associated with inputs U k−1 and outputs Y k−1 become zero.It follows then that we have HY k−1 LU k−1 x k−h .This completes the proof.
According to the deadbeat property, the LMS RH estimator tracks down the real state exactly if no noise is applied.
The prediction problem for h ≤ −1 can be easily solved from the previous results for h 0. Since x k−h is represented as 2.27 we have

Iterative Computation
The LMS FM estimator 2.25 is of the compact and simple form.However, this form requires the inverse of big matrices, which may lead to long computation time and large numerical errors.To overcome these weak points, in this section, we provide an effective iterative form of 2.25 .First, we represent 2.25 in another form for getting recursive equations.Recalling the following fact:

3.1
Mathematical Problems in Engineering 9 we obtain

3.2
Note that nonsingularity of C T N Π −1 N C N is guaranteed for the observability of A, C .From 3.2 , it can be easily seen that the LMS RH 2.25 can be represented as

3.3
Now, we consider two recursive equations for the batch form 3.3 .One is for obtaining The first one is computed in a backward time and the second one in a forward time.Next subsections deal with each recursive equation.

Recursive Equation for Backward Computation
Here, how to obtain 3 will be discussed.These values will be computed in a backward time.
will be denoted by β 0|k and P 0 , respectively, for consistency with the next section.The following theorem provides the main result.

be computed recursively as follows:
where P N and α N|k are obtained from , and P 1 and α 1|k are zero matrices with appropriate dimensions.
Proof.Before going into a main proof, we introduce some variables C j , N j , Π j , and P j as 3.9 for 2 ≤ j ≤ N. In terms of Π j in 3.10 and P j in 3.11 , β 0|k and P 0 can be represented as where I o , S o j , and B o j are obtained by removing the first row blocks of I, S j , and B j , respectively, and Y j , U j , and α j|k are given by for 2 ≤ j ≤ N. In order to obtain P 0 and β 0|k , we have only to know P N and α N|k .Now, in order to get P N and α N|k , we try to find recursive equations for P j and α j|k on 1 ≤ j ≤ N. By using recursions in 3.8 and 3.9 , we have the following equality: where Δ j is given by Pre-and postmultiplying 3.15 by respectively, we have 3.7 .Note that 3.6 and 3.7 hold for i ≥ 2. From 3.11 and 3.14 , P 2 and α 2|k can be written as

3.21
If P 1 and α 1|k are set to zero matrices with appropriate dimensions, P 2 in 3.20 and α 2|k in 3.21 can be calculated from 3.6 and 3.7 .Thus, we can say that 3.6 and 3.7 hold for i ≥ 1 and are initiated with P 1 0 and α 1|k 0. After obtaining P N from 3.6 , we can calculate P 0 from 3.4 .β 0|k in 3.5 can be obtained from α N|k that comes from 3.7 .
This completes the proof.

Recursive Equation for Forward Computation
Here, the recursive equation for 3.3 on the horizon is derived under the assumption that given.P 0 and β 0|k in Section 3.1 are computed in a backward time while variables introduced in this section are computed in a forward time.Before proceeding to a main result, we introduce some variables and the necessary lemma.
L G,i , M i , and N i are defined as 23 where 2 ≤ i ≤ N and G o i is the matrix obtained by removing the last zero column block from The following lemma shows how variables L G,i , N i , and M i are related to one another.

Lemma 3.2. The following relation is satisfied:
where L G,i , N i , and M i are defined in 3.22 -3.24 , respectively.
Proof.If Γ i is defined as

Mathematical Problems in Engineering 13
The four block elements in M i 1 can be expressed recursively as

3.29
Using 3.28 and 3.29 , we have This completes the proof.
Lemma 3.2 is useful for breaking up big matrices of 3.3 into small matrices.We now exploit the recursive equations for two following quantities: where L b,j and L s,j are defined in a form 2.24 , and Y j−1 , U j−1 , L B,j , L S,j , and Y m,j−1 are given by

3.32
Note that and x k−h|k γ N|k .However, β j|k is recursively computed from j 0 to j N and γ j|k from j N − h to j N with the initial value γ N−h|k β N−h|k .Even though β j|k does not look useful on k − h, k , it is still used to compute γ j|k on that horizon.Now, we try to find out recursive equations for β j|k and γ j|k in what follows.

3.2.1.
Recursive Equation for β j|k on 0 ≤ j ≤ N Using Lemma 3.2, we will obtain a recursive equation for β j|k , which is introduced in the following theorem.Theorem 3.3.On 0 ≤ j ≤ N, β j|k in 3.30 can be computed as follows: where P j is given by and initial conditions are set to Proof.First, we will obtain a closed-form of P j in 3.34 and then, using the closed-form of P j , we show that β j|k in 3.30 can be computed recursively from 3.33 .
Using the defined variables 3.22 and 3.24 , we assume that P j in 3.34 is of the form:

3.35
By an induction method, we will prove 3.35 .For the first step, we check 3.35 for j 1.Given the initial value P 0 , we transform P 1 to the form 3.35 by 3.34 as follows: A s G T .

3.36
Equation 3.38 can be written in terms of L G,i and N i as . Now, we check P j 1 under the assumption that P j L G,j N −1 j L T G,j .From 3.34 , we have

3.37
Thus, we can see that P j in 3.34 can be represented as the form 3.35 in terms of L G,j and N j .Using this result, we are in a position to show that β j|k in 3.30 can be computed recursively from 3.33 .We can rewrite β j|k in 3.30 as where T j is given by

3.39
As in a derivation of P j in a batch form, we show by an induction method that β j|k in 3.38 is equivalent to the one in 3.33 .First, we show that β 1|k can be obtained from where, in terms of L G,i , N i , and T i , β 1|k can be written as Next, we show that β j 1|k in the form 3.38 can be obtained from β j|k of the batch form, that is, β j|k L G,j N −1 j T j L B,j U j−1 L S,j Y j−1 .First, note that we have the following relations:

3.42
where G and Q s in the first and second matrix blocks on the right-hand side of the first equality have no effect on the equation.This completes the proof.
It is observed that the recursive equations with 3.33 and 3.34 are the same as the Kalman filter with initial conditions P 0 β j|k on 1 ≤ j ≤ N provides initial values and inputs for a recursive equation of γ j|k 3.31 , which will be investigated in what follows.

Recursive Equation for γ j|k on N − h ≤ j ≤ N
Here, we discuss the recursive equation for γ j|k on N − h ≤ j ≤ N. As mentioned before, the recursive equation for γ j|k starts from j N − h with the initial value γ N−h|k β N−h|k .On N − h ≤ j ≤ N, γ j|k can be recursively computed with the help of β j|k .However, γ N|k x k−h|k is what we want to find out finally.The recursive equation for γ j|k is given in the following theorem.

Simulation
In this section, a numerical example is given to demonstrate the performance of the proposed LMS RH estimator.Suppose that we have a state space model represented as The system noise covariance Q and the measurement noise covariance R are set to 0.01 2 and 0.027 2 , respectively.The memory size and the fixed-lag size are taken as N 10 and h 3, respectively.A sinusoidal input is applied as an input.We carry out a simulation for the system 4.1 with temporary modeling uncertainties 4.2 .In Figure 3, we compare the estimation errors of the LMS RH estimator with the fixedlag Kalman estimator 18 .When uncertainties do not exist, the fixed-lag Kalman estimator has the smaller estimation error than the proposed LMS RH estimator.It can, however, be seen that the estimation error of the LMS RH estimator is considerably smaller than that of the fixed-lag Kalman estimator when modeling uncertainties exist.Actually, one of poles of the fixed-lag Kalman estimator is so close to a unit that even small uncertainties have a good chance of divergence.Additionally, we can see that the estimation error of the LMS RH estimator converges much more rapidly than that of the fixed-lag Kalman estimator after temporary modeling uncertainty disappears.The slow response is also related to the pole Mathematical Problems in Engineering near a unit.To be summarized, we can say that the proposed LMS RH estimator is more robust than other estimators with infinite memory when applied to systems with modeling uncertainties.

Conclusions
In this paper, we proposed a receding horizon RH estimator based on the mean-squareerror criterion for a discrete-time state space model, called a least-mean-square LMS RH estimator.An unknown state was estimated by making use of the finite number of inputs and outputs over the recent finite horizon without any arbitrary assumptions and any a priori state information.The proposed LMS RH estimator was obtained from the conditional expectations of the initial state and the system noise on the corresponding horizon.It was shown that the LMS RH estimator has a deadbeat property and has good robust performance through a numerical example.
To the best of authors' knowledge, the proposed LMS RH estimator would be the most general version among existing RH or finite-memory estimators in the mean-square-error sense.Furthermore, the LMS RH estimator could be extended to other stochastic systems with imperfect communications, uncertainties, and so on 19, 20 .

Figure 3 :
Figure 3: Estimation errors of LMS RH and Kalman estimators when temporary uncertainties exist.