# Least-Mean-Square Receding Horizon Estimation

Authors: Bokyu Kwon; Soohee Han
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

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. Anya 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.


---

## Body

## 1. 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 timek, the least-mean-square (LMS) RH estimator is to estimate the state xk-h at the time k-h from the inputs uk-· and the outputs yk-· over the recent finite horizon [k-N,k-1], which can be written as(1.1)x̂k-h∣k=E[xk-h∣uk-N⋯uk-1,yk-N⋯yk-1],
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[(x̂k-h∣k-xk-h)T(x̂k-h∣k-xk-h)∣uk-N⋯uk-1,yk-N⋯yk-1]. 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 finite-memory or finite-impulse-response (FIR) filters have been preferable for practical reasons [5].Figure 1
The structure of the receding horizon estimator (Dis a unit delay component).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 mean-square-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 Section2, 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.


## 2. LMS RH Estimator

Consider a linear discrete-time state space model with an external control input:(2.1)xi+1=Axi+Bui+Gwi,yi=Cxi+vi,
where xi∈ℜn, ui∈ℜl, and yi∈ℜq are the state, the input, and the output, respectively, and the system noise wi∈ℜp and the measurement noise vi∈ℜq are assumed to be zero-mean white Gaussian with(2.2)E[[wivi][wiTviT]]=[QSSTR],Q≥0,R>0,
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(2.3)xi+1=Asxi+Bui+Gws,i+GSR-1yi,yi=Cxi+vi,
where As≜A-GSR-1C and ws,i≜wi-SR-1vi. It can be easily shown that ws,i and vi in (2.3) are not correlated. In other words, we have(2.4)E[[ws,ivi][ws,iTviT]]=[Q-SR-1STOOR],
where Q-SR-1ST≥0. It is noted that off-diagonal blocks are filled with zeros.For simple representations, several variables are defined as(2.5)C̃N≜[CCAsCAs2⋮CAsN-1],Yk-1≜[yk-Nyk-N+1yk-N+2⋮yk-1],(2.6)B̃N≜[00⋯00CB0⋯00CAsBCB⋯00⋮⋮⋮⋮⋮CAsN-2BCAsN-3B⋯CB0],(2.7)QN≜Q-SR-1ST⊕⋯⊕Q-SR-1ST,︷N(2.8)RN≜R⊕R⊕⋯⊕R,︷N
where ⊕ denotes the direct sum of the matrices. Additionally, Wk-1,Vk-1,  and Uk-1 are defined by replacing y  in (2.5) with ws, 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 of0≤h≤N-1 in (1.1). A prediction problem for h≤-1 will be discussed later on.By using the defined variables, the statexk-h to be estimated is represented in terms of the initial state xk-N on the horizon, inputs, and system noise on the recent finite horizon [k-N,  k-1] as(2.9)xk-h=AsN-hxk-N+Lg,NWk-1+Lb,NUk-1+Ls,NYk-1,
where Lg,N is given by(2.10)Lg,N≜[AsN-h-1GAsN-h-2G⋯GO⋯O︷h],
and Lb,N and Ls,N are obtained by replacing G in Lg,N with B and GSR-1, respectively. From (2.9), the conditional expectation E[xk-h|Uk-1,Yk-1] in (1.1) can be represented as
(2.11)E[xk-h∣Uk-1,Yk-1]=AsN-hE[xk-N∣Uk-1,Yk-1]+Lg,NE[Wk-1∣Uk-1,Yk-1]+Lb,NUk-1+Ls,NYk-1.
Note that if E[xk-N∣Uk-1,Yk-1] and E[Wk-1∣Uk-1,Yk-1] are known, E[xk-h∣Uk-1,Yk-1] can be obtained. In order to compute E[xk-N∣Uk-1,Yk-1] and E[Wk-1∣Uk-1,Yk-1], we first try to find the relation among xk-N, Wk-1, Uk-1, and Yk-1. For this purpose, the system (2.3) is represented in a batch form on the recent finite horizon [k-N,k-1] as follows:(2.12)Yk-1=C̃Nxk-N+B̃NUk-1+G̃NWk-1+S̃NYk-1+Vk-1,
from which we obtain(2.13)[Wk-1Vk-1]=[I0-G̃N-C̃N][Wk-1xk-N]+[0Yk-1-S̃NYk-1-B̃NUk-1].
We can see from (2.13) that Wk-1 and Vk-1 are linearly, more correctly affinely, dependent on Wk-1 and xk-N. The joint probability density function of Wk-1 and xk-N, that is, fwx(Wk-1,xk-N), can be expressed as kfwv(Wk-1,Vk-1∣Uk-1,Yk-1)=kfwv(Wk-1,Yk-1-C̃Nxk-N-B̃NUk-1-G̃NWk-1-S̃NYk-1) for an appropriate scaling factor k. How to choose k will be discussed later on. By using the probability density function of Wk-1 and Vk-1 given as
(2.14)fwv(Wk-1,Vk-1)=1(2π)(p+q)Ndet⁡(QN)det⁡(RN)×exp⁡(-12[Wk-1Vk-1]T[QN00RN]-1[Wk-1Vk-1]),E[xk-N∣Uk-1,Yk-1] can be computed from(2.15)E[xk-N∣Uk-1,Yk-1]=∫xk-Nfwx(Wk-1,xk-N∣Uk-1,Yk-1)dWk-1dxk-N=∫xk-Nkfwv(Wk-1,Yk-1-C̃Nxk-N-B̃NUk-1-G̃NWk-1-S̃NYk-1)dWk-1dxk-N,
where a constant k is chosen to satisfy the following normalization condition:(2.16)∫fwx(Wk-1,xk-N∣Uk-1,Yk-1)dWk-1dxk-N=k∫fwv(Wk-1,Yk-1-C̃Nxk-N-B̃NUk-1-G̃NWk-1-S̃NYk-1)dWk-1dxk-N=1.
In the same way, E[Wk-1∣Uk-1,Yk-1] can be obtained. It is noted that fwv(Wk-1,Yk-1-C̃Nxk-N-B̃NUk-1-G̃NWk-1-S̃NYk-1) is an exponential function with the following exponent:(2.17)[Yk-1T(I-S̃N)TUk-1Txk-NTWk-1T][RN-1-RN-1B̃N-RN-1G̃N-RN-1C̃N*B̃NTRN-1B̃NB̃NTRN-1G̃NB̃NTRN-1C̃N**W1,1W1,2***W2,2][Yk-1T(I-S̃N)TUk-1Txk-NTWk-1T]T,
where * denotes symmetric parts and W1,1, W1,2, and W2,2 are given by(2.18)W1,1=C̃NTRN-1C̃N,W1,2=C̃NTRN-1G̃N,W2,2=G̃NTRN-1+QN-1.
By completing the square of the terms of the exponent and recalling the integration of Gaussian functions from -∞ to ∞, we can compute E[xk-N∣Uk-1,Yk-1] and E[Wk-1∣Uk-1,Yk-1]. To begin with, we introduce the following completion of squares:(2.19)[ab]T[αββTγ][ab]=(b+γ-1βTa)Tγ(b+γ-1βTa)+aT(α-βγ-1βT)a
to obtain(2.20)∫bexp⁡([ab]T[αββTγ][ab])db∫exp⁡([ab]T[αββTγ][ab])db=-γ-1βTa,
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)a⟵[Yk-1T(I-S̃N)TUk-1T],b⟵[xk-NTWk-1T],α⟵[RN-1-RN-1B̃N*B̃NTRN-1B̃N],β⟵[-RN-1G̃N-RN-1C̃NB̃NTRN-1G̃NB̃NTRN-1C̃N],γ⟵[W1,1W1,2*W2,2],
we have(2.22)[E[xk-N∣Uk-1,Yk-1]E[Wk-1∣Uk-1,Yk-1]]=[W1,1W1,2W1,2TW2,2]-1[C̃NG̃N]×RN-1((I-S̃N)Yk-1-B̃NUk-1),
where W1,1, W1,2, and W2,2 are given by (2.18). It is noted that if (A,C) of the system (2.1) is observable, W1,1 is positive definite and W2,2-W1,2TW1,1-1W1,2 is also positive definite, which implies that the block matrix including W1,1, W1,2, and W2,2 is guaranteed to be positive definite and hence nonsingular. Substituting (2.22) into (2.11), we can obtain E[xk-h∣Uk-1,Yk-1]. Through a long and tedious algebraic calculation, we have the covariance matrix of the mean-square-error E[(x̂k-h∣k-xk-h)(x̂k-h∣k-xk-h)T] as follows:(2.23)P=ΞN[W1,1W1,2W1,2TW2,2]-1ΞNT,
where ΞN is defined by(2.24)Ξj≜[AsN-hAsN-h-1GAsN-h-2G⋯GOO⋯O︷j-N+h],
for N-h≤j≤N.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)x̂k-h|k=ΞN[W1,1W1,2W1,2TW2,2]-1[C̃NG̃N]RN-1((I-S̃N)Yk-1-B̃NUk-1)+Lb,NUk-1+Ls,NYk-1,
where x̂k-h∣k denotes E[xk-h∣Uk-1,Yk-1] and W1,1, W1,2, and W2,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, xk-h|k=xk-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 Yk-1 and Uk-1, that is, HYk-1+LUk-1 for some gain matrices H and L. Substituting Yk-1 in (2.12) into HYk-1+LUk-1, using (2.9), and removing all noise, we have
(2.26)HYk-1+LUk-1=(HC̃N-AsN-h)xk-N+xk-h+(HB̃N-Lb,N+L)Uk-1+(HS̃N-Ls,N)Yk-1.
It can be easily seen that HC̃N=AsN-h and terms associated with inputs Uk-1 and outputs Yk-1 become zero. It follows then that we have HYk-1+LUk-1=xk-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 forh≤-1 can be easily solved from the previous results for h=0. Since xk-h is represented as(2.27)xk-h=A-hxk+∑i=0-h-1AiBuk+i+∑i=0-h-1AiGwk+i,
we have(2.28)E[xk-h∣Uk-1,Yk-1]=A-hE[xk∣Uk-1,Yk-1]+∑i=0-h-1AiBuk+i,
which means that E[xk-h∣Uk-1,Yk-1] can be obtained from E[xk∣Uk-1,Yk-1]. Setting h in (2.25) to zero, we can compute E[xk∣Uk-1,Yk-1] easily.


## 3. 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)(C̃NTΠN-1C̃N)-1=(W1,1-W1,2W2,2-1W1,2T)-1,C̃NTΠN-1=C̃NTRN-1-W1,2W2,2-1G̃NTRN-1,
we obtain(3.2)[C̃NTΠN-10]+[C̃NTG̃NT]RN-1=[2C̃NTRN-1-W1,2W2,2-1G̃NTRN-1G̃NTRN-1]=[2I-W1,2W2,2-10I][C̃NTG̃NT]RN-1,[W1,1+C̃NTΠN-1C̃NW2,2W1,2TW2,2]-1[2I-W1,2W2,2-10I]={[12I12W1,2W2,2-10I][W1,1+C̃NTΠN-1C̃NW1,2W1,2TW2,2]}-1=[W1,1W1,2W1,2TW2,2]-1.
Note that nonsingularity of C̃NTΠN-1C̃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)x̂k-h|k=ΞN[W1,1+C̃NTΠN-1C̃NW1,2W1,2TW2,2]-1×([C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1)0]+[C̃NTG̃NT]RN-1((I-S̃N)Yk-1-B̃NUk-1))+Ls,NYk-1+Lb,NUk-1.Now, we consider two recursive equations for the batch form (3.3). One is for obtaining C̃NTΠN-1C̃N and C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1). The other is for (3.3) given C̃NTΠN-1C̃N and C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1). The first one is computed in a backward time and the second one in a forward time. Next subsections deal with each recursive equation.


### 3.1. Recursive Equation for Backward Computation

Here, how to obtainC̃NTΠN-1C̃N and C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1) in (3.3) will be discussed. These values will be computed in a backward time. (C̃NTΠN-1C̃N)-1C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1) and (C̃NTΠN-1C̃N)-1 will be denoted by β0|k and P0, respectively, for consistency with the next section. The following theorem provides the main result.Theorem 3.1.
P0=(C̃NTΠN-1C̃N)-1 and β0|k=P0C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1) can be computed recursively as follows:
(3.4)P0=(CTR-1C+P̂N)-1,(3.5)β0|k=P0(CTR-1yk-N+αN|k),
where P̂N and αN∣k are obtained from
(3.6)P̂j+1=AsTCTR-1CAs+AsTP̂jAs-AsT(CTR-1C+P̂j)G×{Qs-1+GT(CTR-1C+P̂j)G}-1GT(CTR-1C+P̂j)As,(3.7)αj+1∣k={AsT-AsT(CTR-1C+P̂j)G{Qs-1+GT(CTR-1C+P̂j)G}-1GT}×{αj∣k+CTR-1yk-j-(P̂j+CTR-1C)(Buk-j-1+GSR-1yk-j-1)},
for 1≤j≤N-1, Qs=Q-SR-1ST,  and P̂1 and α1∣k are zero matrices with appropriate dimensions.Proof.
Before going into a main proof, we introduce some variablesĈj, N̂j, Π̂j, and P̂j as
(3.8)Ĉj≜[CAs⋮CAsj-3CAsj-2CAsj-1]=[CĈj-1]A,(3.9)N̂j≜[CG0⋯0CAsGCG⋯0⋮⋮⋮⋮CAsj-2GCAsj-3G⋯CG]=[CG0Ĉj-1GN̂j-1],(3.10)Π̂j≜N̂jQj-1N̂jT+Rj-1∈Rq(j-1)×q(j-1),(3.11)P̂j≜ĈjTΠ̂j-1Ĉj,
for 2≤j≤N. In terms of Π̂j in (3.10) and P̂j in (3.11), β0∣k and P0 can be represented as
(3.12)β0∣k=P0(CTR-1yk-N+ĈNTΠ̂N-1((Io-S̃No)ŶN-B̃NoÛN))=P0(CTR-1yk-N+αN∣k),P0=(CTR-1C+P̂N)-1,
where Io, S̃jo, and B̃jo 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
(3.13)Ŷj≜[yk-jT⋯yk-1T]T,Ûj≜[uk-jT⋯uk-1T]T,(3.14)αj∣k≜ĈjTΠ̂j-1((Io-S̃jo)Ŷj-B̃joÛj),
for 2≤j≤N. In order to obtain P0 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:
(3.15)Π̂j+1-1=Δ̂j-1-Δ̂j-1[CĈj]G{Qs-1+GT[CĈj]TΔ̂j-1[CĈj]G}-1×GT[CĈj]TΔ̂j-1,
where Δ̂j is given by
(3.16)Δ̂j≜[R00Π̂j].
Pre- and postmultiplying (3.15) by
(3.17)Ĉj+1T=AsT[CĈj]T,Ĉj+1=[CĈj]As,
respectively, we have (3.6). Pre- and postmultiplying (3.15) by
(3.18)Ĉj+1T=AsT[CĈj]T,(3.19)(Io-S̃j+1o)Ŷj+1-B̃j+1oÛj+1=[yk-j(Io-S̃jo)Ŷj-B̃joÛj]-[CĈj](Buk-j-1+GSR-1yk-j-1),
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.20)P̂2=Ĉ2TΠ̂2-1Ĉ2=AsTCT(CGQsGTCT+R)-1CAsT,(3.21)α2∣k=Ĉ2TΠ̂2-1((Io-S̃2o)Ŷ2-B̃2oÛ2)=AsTCT(CGQsGTCT+R)-1(yk-1-CGSR-1yk-2-CBuk-2).
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 P0 from (3.4). β0∣k in (3.5) can be obtained from αN|k that comes from (3.7).
This completes the proof.


### 3.2. Recursive Equation for Forward Computation

Here, the recursive equation for (3.3) on the horizon is derived under the assumption that P0=(C̃NTΠN-1C̃N)-1 and β0∣k=P0C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1) are given. P0 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.LG,i, Mi, and Ni are defined as(3.22)LG,i≜[AsiAsi-1GAsi-2G⋯AsGG],(3.23)Mi≜[C̃iTRi-1C̃i+P0-1C̃iTRi-1G̃ioG̃ioTRi-1C̃iG̃ioTRi-1G̃io+Qi-1-1],(3.24)Ni≜[Mi00Qs-1]=[C̃iTRi-1C̃i+P0-1C̃iTRi-1G̃iG̃iTRi-1C̃iG̃iTRi-1G̃i+Qi-1],
where 2≤i≤N and G̃io is the matrix obtained by removing the last zero column block from G̃i, that is, G̃i=[G̃io0]. In particular, LG,1, M1, and N1 are defined as [AsG], CTR-1C+P0-1 and M1⊕Qs-1. The following lemma shows how variables LG,i, Ni, and Mi are related to one another.Lemma 3.2.
The following relation is satisfied:(3.25)Ni+LG,iTCTR-1CLG,i=Mi+1,
where LG,i,  Ni,  and Mi are defined in (3.22)–(3.24), respectively.Proof.
IfΓi is defined as
(3.26)Γi≜[Asi-1GAsi-2G⋯AsGG],LG,i can be represented as
(3.27)LG,i=[AsiΓi],
from which we have
(3.28)LG,iTCTR-1CLG,i=[AsTiCTR-1CAsiAsTiCTR-1CΓiΓiTCTR-1CAsiΓiTCTR-1CΓi].
The four block elements inMi+1 can be expressed recursively as
(3.29)C̃i+1TRi+1-1C̃i+1=C̃iTRi-1C̃i+AsTiCTR-1CAsi,C̃i+1TRi+1-1G̃i+1o=[C̃iTRi-1G̃io0]+AsTiCTR-1CΓi,G̃i+1oTRi+1-1G̃i+1o=[G̃ioTRi-1G̃io000]+ΓiTCTR-1CΓi.
Using (3.28) and (3.29), we have Ni+LG,iTCTR-1CLG,i=Mi+1. This completes the proof.Lemma3.2 is useful for breaking up big matrices of (3.3) into small matrices. We now exploit the recursive equations for two following quantities:(3.30)βj∣k=LG,jNj-1([C̃NTΠN-1Ỹm,N-10]+[C̃jTG̃jT]Rj-1Ỹm,j-1)+LB,jŨj-1+LS,jỸj-1,
for 1≤j≤N(3.31)γj∣k=ΞjNj-1([C̃NTΠN-1Ỹm,N-10]+[C̃jTG̃jT]Rj-1Ỹm,j-1)+Lb,jŨj-1+Ls,jỸj-1,
for N-h≤j≤N, where Lb,j and Ls,j are defined in a form (2.24), and Ỹj-1, Ũj-1, LB,j, LS,j, and Ỹm,j-1 are given by(3.32)Ỹj-1≜[yk-NT⋯yk-N+j-1T]T,Ũj-1≜[uk-NT⋯uk-N+j-1T]T,LB,j≜[Asj-1BAsj-2B⋯B],LS,j≜[Asj-1GSR-1Asj-2GSR-1⋯GSR-1],Ỹm,j-1≜(I-S̃j)Ỹj-1-B̃jŨj-1.
Note that β0∣k=C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1), γN-h∣k=βN-h∣k, 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 Lemma3.2, we will obtain a recursive equation for βj∣k, which is introduced in the following theorem.Theorem 3.3.
On0≤j≤N, βj∣k in (3.30) can be computed as follows:
(3.33)βj+1∣k=Asβj∣k+Buk-N+j+GSR-1yk-N+j+AsPjCT(R+CPjCT)-1(yk-N+j-Cβj|k),
where Pj is given by
(3.34)Pj+1=AsPjAsT+GQsGT-AsPjCT(R+CPjCT)-1CPjAsT,
and initial conditions are set to P0=(C̃NTΠN-1C̃N)-1 and β0|k=P0C̃NTΠN-1Ỹm,N-1 computed in Section 3.1.Proof.
First, we will obtain a closed-form ofPj in (3.34) and then, using the closed-form of Pj, 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 Pj in (3.34) is of the form:
(3.35)Pj=LG,jNj-1LG,jT.
By an induction method, we will prove (3.35). For the first step, we check (3.35) for j=1. Given the initial value P0, we transform P1 to the form (3.35) by (3.34) as follows:
(3.36)P1=[AsG][P0-1+CTR-1C00Qs-1]-1[AsG]T.
Equation (3.38) can be written in terms of LG,i and Ni as P1=LG,1N1-1LG,1T. Now, we check Pj+1 under the assumption that Pj=LG,jNj-1LG,jT. From (3.34), we have
(3.37)Pj+1=AsLG,jMj+1-1LG,jTAsT+GQsGT=LG,j+1Nj+1-1LG,j+1T.
Thus, we can see that Pj in (3.34) can be represented as the form (3.35) in terms of LG,j and Nj. 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
(3.38)βj∣k=LG,jNj-1Tj+LB,jŨj-1+LS,jỸj-1,
where Tj is given by
(3.39)Tj≜[P0-10]β0|k+[C̃jTG̃jT]Rj-1Ỹm,j-1.
As in a derivation of Pj 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 P0=(C̃NTΠN-1C̃N)-1 and β0∣k=P0C̃NTΠN-1Ỹm,N-1:
(3.40)β1∣k=As(P0-1+CTR-1C)-1(P0-1β0∣k+CTR-1yk-N)+Buk-N+GSR-1yk-N,
where, in terms of LG,i, Ni, and Ti, β1|k can be written as β1|k=LG,1N1-1T1+LB,1Ũ0+LS,1Ỹ0. 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=LG,jNj-1Tj+LB,jŨj-1+LS,jỸj-1. First, note that we have the following relations:
(3.41){As-AsPjCT(R+CPjCT)-1C}LG,jNj-1Tj=AsLG,jMj+1-1Tj,AsPjCT(R+CPjCT)-1yk-N+j=AsLG,jMj+1-1LG,jTCTR-1yk-N+j.
Substituting (3.41) into (3.33) yields
(3.42)βj+1∣k=AsLG,jMj+1-1Tj+AsLG,jMj+1-1LG,jTCTR-1yk-N+j-AsLG,jMj+1-1LG,jTCTR-1C(LB,jŨj-1+LS,jỸj-1)+LB,j+1Ũj+LS,j+1Ỹj,=[AsLG,jG][Mj+1-100Qs]{[P0-10]β0∣k+[C̃j+1T[G̃j+1o0]T]Rj+1-1Ỹm,j}+LB,j+1Ũj+LS,j+1Ỹj=LG,j+1Nj+1-1Tj+1+LB,j+1Ũj+LS,j+1Ỹj,
where G and Qs 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 P0=(C̃NTΠN-1C̃N)-1 and β0∣k=P0C̃NTΠN-1Ỹm,N-1. β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.


#### 3.2.2. 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.Theorem 3.4.
OnN-h≤j≤N-1, γj|k in (3.31) can be computed as follows:
(3.43)γj+1∣k=γj∣k+KjCT(CPjCT+R)-1(yk-N+j-Cβj∣k),
where γN-h∣k=βN-h|k, βj∣k and Pj are obtained from (3.33) and (3.34), respectively, and 𝒦j is given by
(3.44)Kj=Kj-1(I-CT(CPj-1CT+R)-1CPj-1)AT,
with the initial condition 𝒦N-h=PN-h.Proof.
Using Lemma3.2, γj+1∣k in (3.31) for j≥N-h can be represented as
(3.45)γj+1∣k=γj∣k+ΞjNj-1LG,jTCT(CPjCT+R)-1(yk-N+j-Cβj∣k).
If we denote ΞjNj-1LG,jT by 𝒦j, we have only to prove (3.44). Since ΞN-h=LG,N-h, 𝒦N-h is equal to PN-h. Using Lemma 3.2, we can represent 𝒦j in a recursive form as follows:
(3.46)Kj=Kj-1(I-CT(CPjCT+R)-1CPj-1)AT.
This completes the proof.It is noted that the recursive equation (3.43) in Theorem 3.4 is a fixed-point smoother of the state xk-h, which runs with a recursive equation (3.33) in Theorem 3.3. Variables for recursive equations in Sections 3.1 and 3.2 are visualized in Figure 2. Starting from P̂1=0 and α1∣k=0, we compute P̂N and αN∣k recursively in a backward time. From P̂N and αN∣k, we compute P0 and β0∣k, from which we drive the forward recursive equation for βj∣k to get γN∣k=x̂k-h∣k.Figure 2
An iterative form of the LMS RH estimator.


## 3.1. Recursive Equation for Backward Computation

Here, how to obtainC̃NTΠN-1C̃N and C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1) in (3.3) will be discussed. These values will be computed in a backward time. (C̃NTΠN-1C̃N)-1C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1) and (C̃NTΠN-1C̃N)-1 will be denoted by β0|k and P0, respectively, for consistency with the next section. The following theorem provides the main result.Theorem 3.1.
P0=(C̃NTΠN-1C̃N)-1 and β0|k=P0C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1) can be computed recursively as follows:
(3.4)P0=(CTR-1C+P̂N)-1,(3.5)β0|k=P0(CTR-1yk-N+αN|k),
where P̂N and αN∣k are obtained from
(3.6)P̂j+1=AsTCTR-1CAs+AsTP̂jAs-AsT(CTR-1C+P̂j)G×{Qs-1+GT(CTR-1C+P̂j)G}-1GT(CTR-1C+P̂j)As,(3.7)αj+1∣k={AsT-AsT(CTR-1C+P̂j)G{Qs-1+GT(CTR-1C+P̂j)G}-1GT}×{αj∣k+CTR-1yk-j-(P̂j+CTR-1C)(Buk-j-1+GSR-1yk-j-1)},
for 1≤j≤N-1, Qs=Q-SR-1ST,  and P̂1 and α1∣k are zero matrices with appropriate dimensions.Proof.
Before going into a main proof, we introduce some variablesĈj, N̂j, Π̂j, and P̂j as
(3.8)Ĉj≜[CAs⋮CAsj-3CAsj-2CAsj-1]=[CĈj-1]A,(3.9)N̂j≜[CG0⋯0CAsGCG⋯0⋮⋮⋮⋮CAsj-2GCAsj-3G⋯CG]=[CG0Ĉj-1GN̂j-1],(3.10)Π̂j≜N̂jQj-1N̂jT+Rj-1∈Rq(j-1)×q(j-1),(3.11)P̂j≜ĈjTΠ̂j-1Ĉj,
for 2≤j≤N. In terms of Π̂j in (3.10) and P̂j in (3.11), β0∣k and P0 can be represented as
(3.12)β0∣k=P0(CTR-1yk-N+ĈNTΠ̂N-1((Io-S̃No)ŶN-B̃NoÛN))=P0(CTR-1yk-N+αN∣k),P0=(CTR-1C+P̂N)-1,
where Io, S̃jo, and B̃jo 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
(3.13)Ŷj≜[yk-jT⋯yk-1T]T,Ûj≜[uk-jT⋯uk-1T]T,(3.14)αj∣k≜ĈjTΠ̂j-1((Io-S̃jo)Ŷj-B̃joÛj),
for 2≤j≤N. In order to obtain P0 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:
(3.15)Π̂j+1-1=Δ̂j-1-Δ̂j-1[CĈj]G{Qs-1+GT[CĈj]TΔ̂j-1[CĈj]G}-1×GT[CĈj]TΔ̂j-1,
where Δ̂j is given by
(3.16)Δ̂j≜[R00Π̂j].
Pre- and postmultiplying (3.15) by
(3.17)Ĉj+1T=AsT[CĈj]T,Ĉj+1=[CĈj]As,
respectively, we have (3.6). Pre- and postmultiplying (3.15) by
(3.18)Ĉj+1T=AsT[CĈj]T,(3.19)(Io-S̃j+1o)Ŷj+1-B̃j+1oÛj+1=[yk-j(Io-S̃jo)Ŷj-B̃joÛj]-[CĈj](Buk-j-1+GSR-1yk-j-1),
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.20)P̂2=Ĉ2TΠ̂2-1Ĉ2=AsTCT(CGQsGTCT+R)-1CAsT,(3.21)α2∣k=Ĉ2TΠ̂2-1((Io-S̃2o)Ŷ2-B̃2oÛ2)=AsTCT(CGQsGTCT+R)-1(yk-1-CGSR-1yk-2-CBuk-2).
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 P0 from (3.4). β0∣k in (3.5) can be obtained from αN|k that comes from (3.7).
This completes the proof.


## 3.2. Recursive Equation for Forward Computation

Here, the recursive equation for (3.3) on the horizon is derived under the assumption that P0=(C̃NTΠN-1C̃N)-1 and β0∣k=P0C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1) are given. P0 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.LG,i, Mi, and Ni are defined as(3.22)LG,i≜[AsiAsi-1GAsi-2G⋯AsGG],(3.23)Mi≜[C̃iTRi-1C̃i+P0-1C̃iTRi-1G̃ioG̃ioTRi-1C̃iG̃ioTRi-1G̃io+Qi-1-1],(3.24)Ni≜[Mi00Qs-1]=[C̃iTRi-1C̃i+P0-1C̃iTRi-1G̃iG̃iTRi-1C̃iG̃iTRi-1G̃i+Qi-1],
where 2≤i≤N and G̃io is the matrix obtained by removing the last zero column block from G̃i, that is, G̃i=[G̃io0]. In particular, LG,1, M1, and N1 are defined as [AsG], CTR-1C+P0-1 and M1⊕Qs-1. The following lemma shows how variables LG,i, Ni, and Mi are related to one another.Lemma 3.2.
The following relation is satisfied:(3.25)Ni+LG,iTCTR-1CLG,i=Mi+1,
where LG,i,  Ni,  and Mi are defined in (3.22)–(3.24), respectively.Proof.
IfΓi is defined as
(3.26)Γi≜[Asi-1GAsi-2G⋯AsGG],LG,i can be represented as
(3.27)LG,i=[AsiΓi],
from which we have
(3.28)LG,iTCTR-1CLG,i=[AsTiCTR-1CAsiAsTiCTR-1CΓiΓiTCTR-1CAsiΓiTCTR-1CΓi].
The four block elements inMi+1 can be expressed recursively as
(3.29)C̃i+1TRi+1-1C̃i+1=C̃iTRi-1C̃i+AsTiCTR-1CAsi,C̃i+1TRi+1-1G̃i+1o=[C̃iTRi-1G̃io0]+AsTiCTR-1CΓi,G̃i+1oTRi+1-1G̃i+1o=[G̃ioTRi-1G̃io000]+ΓiTCTR-1CΓi.
Using (3.28) and (3.29), we have Ni+LG,iTCTR-1CLG,i=Mi+1. This completes the proof.Lemma3.2 is useful for breaking up big matrices of (3.3) into small matrices. We now exploit the recursive equations for two following quantities:(3.30)βj∣k=LG,jNj-1([C̃NTΠN-1Ỹm,N-10]+[C̃jTG̃jT]Rj-1Ỹm,j-1)+LB,jŨj-1+LS,jỸj-1,
for 1≤j≤N(3.31)γj∣k=ΞjNj-1([C̃NTΠN-1Ỹm,N-10]+[C̃jTG̃jT]Rj-1Ỹm,j-1)+Lb,jŨj-1+Ls,jỸj-1,
for N-h≤j≤N, where Lb,j and Ls,j are defined in a form (2.24), and Ỹj-1, Ũj-1, LB,j, LS,j, and Ỹm,j-1 are given by(3.32)Ỹj-1≜[yk-NT⋯yk-N+j-1T]T,Ũj-1≜[uk-NT⋯uk-N+j-1T]T,LB,j≜[Asj-1BAsj-2B⋯B],LS,j≜[Asj-1GSR-1Asj-2GSR-1⋯GSR-1],Ỹm,j-1≜(I-S̃j)Ỹj-1-B̃jŨj-1.
Note that β0∣k=C̃NTΠN-1((I-S̃N)Yk-1-B̃NUk-1), γN-h∣k=βN-h∣k, 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 Lemma3.2, we will obtain a recursive equation for βj∣k, which is introduced in the following theorem.Theorem 3.3.
On0≤j≤N, βj∣k in (3.30) can be computed as follows:
(3.33)βj+1∣k=Asβj∣k+Buk-N+j+GSR-1yk-N+j+AsPjCT(R+CPjCT)-1(yk-N+j-Cβj|k),
where Pj is given by
(3.34)Pj+1=AsPjAsT+GQsGT-AsPjCT(R+CPjCT)-1CPjAsT,
and initial conditions are set to P0=(C̃NTΠN-1C̃N)-1 and β0|k=P0C̃NTΠN-1Ỹm,N-1 computed in Section 3.1.Proof.
First, we will obtain a closed-form ofPj in (3.34) and then, using the closed-form of Pj, 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 Pj in (3.34) is of the form:
(3.35)Pj=LG,jNj-1LG,jT.
By an induction method, we will prove (3.35). For the first step, we check (3.35) for j=1. Given the initial value P0, we transform P1 to the form (3.35) by (3.34) as follows:
(3.36)P1=[AsG][P0-1+CTR-1C00Qs-1]-1[AsG]T.
Equation (3.38) can be written in terms of LG,i and Ni as P1=LG,1N1-1LG,1T. Now, we check Pj+1 under the assumption that Pj=LG,jNj-1LG,jT. From (3.34), we have
(3.37)Pj+1=AsLG,jMj+1-1LG,jTAsT+GQsGT=LG,j+1Nj+1-1LG,j+1T.
Thus, we can see that Pj in (3.34) can be represented as the form (3.35) in terms of LG,j and Nj. 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
(3.38)βj∣k=LG,jNj-1Tj+LB,jŨj-1+LS,jỸj-1,
where Tj is given by
(3.39)Tj≜[P0-10]β0|k+[C̃jTG̃jT]Rj-1Ỹm,j-1.
As in a derivation of Pj 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 P0=(C̃NTΠN-1C̃N)-1 and β0∣k=P0C̃NTΠN-1Ỹm,N-1:
(3.40)β1∣k=As(P0-1+CTR-1C)-1(P0-1β0∣k+CTR-1yk-N)+Buk-N+GSR-1yk-N,
where, in terms of LG,i, Ni, and Ti, β1|k can be written as β1|k=LG,1N1-1T1+LB,1Ũ0+LS,1Ỹ0. 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=LG,jNj-1Tj+LB,jŨj-1+LS,jỸj-1. First, note that we have the following relations:
(3.41){As-AsPjCT(R+CPjCT)-1C}LG,jNj-1Tj=AsLG,jMj+1-1Tj,AsPjCT(R+CPjCT)-1yk-N+j=AsLG,jMj+1-1LG,jTCTR-1yk-N+j.
Substituting (3.41) into (3.33) yields
(3.42)βj+1∣k=AsLG,jMj+1-1Tj+AsLG,jMj+1-1LG,jTCTR-1yk-N+j-AsLG,jMj+1-1LG,jTCTR-1C(LB,jŨj-1+LS,jỸj-1)+LB,j+1Ũj+LS,j+1Ỹj,=[AsLG,jG][Mj+1-100Qs]{[P0-10]β0∣k+[C̃j+1T[G̃j+1o0]T]Rj+1-1Ỹm,j}+LB,j+1Ũj+LS,j+1Ỹj=LG,j+1Nj+1-1Tj+1+LB,j+1Ũj+LS,j+1Ỹj,
where G and Qs 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 P0=(C̃NTΠN-1C̃N)-1 and β0∣k=P0C̃NTΠN-1Ỹm,N-1. β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.


### 3.2.2. 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.Theorem 3.4.
OnN-h≤j≤N-1, γj|k in (3.31) can be computed as follows:
(3.43)γj+1∣k=γj∣k+KjCT(CPjCT+R)-1(yk-N+j-Cβj∣k),
where γN-h∣k=βN-h|k, βj∣k and Pj are obtained from (3.33) and (3.34), respectively, and 𝒦j is given by
(3.44)Kj=Kj-1(I-CT(CPj-1CT+R)-1CPj-1)AT,
with the initial condition 𝒦N-h=PN-h.Proof.
Using Lemma3.2, γj+1∣k in (3.31) for j≥N-h can be represented as
(3.45)γj+1∣k=γj∣k+ΞjNj-1LG,jTCT(CPjCT+R)-1(yk-N+j-Cβj∣k).
If we denote ΞjNj-1LG,jT by 𝒦j, we have only to prove (3.44). Since ΞN-h=LG,N-h, 𝒦N-h is equal to PN-h. Using Lemma 3.2, we can represent 𝒦j in a recursive form as follows:
(3.46)Kj=Kj-1(I-CT(CPjCT+R)-1CPj-1)AT.
This completes the proof.It is noted that the recursive equation (3.43) in Theorem 3.4 is a fixed-point smoother of the state xk-h, which runs with a recursive equation (3.33) in Theorem 3.3. Variables for recursive equations in Sections 3.1 and 3.2 are visualized in Figure 2. Starting from P̂1=0 and α1∣k=0, we compute P̂N and αN∣k recursively in a backward time. From P̂N and αN∣k, we compute P0 and β0∣k, from which we drive the forward recursive equation for βj∣k to get γN∣k=x̂k-h∣k.Figure 2
An iterative form of the LMS RH estimator.


## 3.2.1. Recursive Equation forβj|k on 0≤j≤N

Using Lemma3.2, we will obtain a recursive equation for βj∣k, which is introduced in the following theorem.Theorem 3.3.
On0≤j≤N, βj∣k in (3.30) can be computed as follows:
(3.33)βj+1∣k=Asβj∣k+Buk-N+j+GSR-1yk-N+j+AsPjCT(R+CPjCT)-1(yk-N+j-Cβj|k),
where Pj is given by
(3.34)Pj+1=AsPjAsT+GQsGT-AsPjCT(R+CPjCT)-1CPjAsT,
and initial conditions are set to P0=(C̃NTΠN-1C̃N)-1 and β0|k=P0C̃NTΠN-1Ỹm,N-1 computed in Section 3.1.Proof.
First, we will obtain a closed-form ofPj in (3.34) and then, using the closed-form of Pj, 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 Pj in (3.34) is of the form:
(3.35)Pj=LG,jNj-1LG,jT.
By an induction method, we will prove (3.35). For the first step, we check (3.35) for j=1. Given the initial value P0, we transform P1 to the form (3.35) by (3.34) as follows:
(3.36)P1=[AsG][P0-1+CTR-1C00Qs-1]-1[AsG]T.
Equation (3.38) can be written in terms of LG,i and Ni as P1=LG,1N1-1LG,1T. Now, we check Pj+1 under the assumption that Pj=LG,jNj-1LG,jT. From (3.34), we have
(3.37)Pj+1=AsLG,jMj+1-1LG,jTAsT+GQsGT=LG,j+1Nj+1-1LG,j+1T.
Thus, we can see that Pj in (3.34) can be represented as the form (3.35) in terms of LG,j and Nj. 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
(3.38)βj∣k=LG,jNj-1Tj+LB,jŨj-1+LS,jỸj-1,
where Tj is given by
(3.39)Tj≜[P0-10]β0|k+[C̃jTG̃jT]Rj-1Ỹm,j-1.
As in a derivation of Pj 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 P0=(C̃NTΠN-1C̃N)-1 and β0∣k=P0C̃NTΠN-1Ỹm,N-1:
(3.40)β1∣k=As(P0-1+CTR-1C)-1(P0-1β0∣k+CTR-1yk-N)+Buk-N+GSR-1yk-N,
where, in terms of LG,i, Ni, and Ti, β1|k can be written as β1|k=LG,1N1-1T1+LB,1Ũ0+LS,1Ỹ0. 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=LG,jNj-1Tj+LB,jŨj-1+LS,jỸj-1. First, note that we have the following relations:
(3.41){As-AsPjCT(R+CPjCT)-1C}LG,jNj-1Tj=AsLG,jMj+1-1Tj,AsPjCT(R+CPjCT)-1yk-N+j=AsLG,jMj+1-1LG,jTCTR-1yk-N+j.
Substituting (3.41) into (3.33) yields
(3.42)βj+1∣k=AsLG,jMj+1-1Tj+AsLG,jMj+1-1LG,jTCTR-1yk-N+j-AsLG,jMj+1-1LG,jTCTR-1C(LB,jŨj-1+LS,jỸj-1)+LB,j+1Ũj+LS,j+1Ỹj,=[AsLG,jG][Mj+1-100Qs]{[P0-10]β0∣k+[C̃j+1T[G̃j+1o0]T]Rj+1-1Ỹm,j}+LB,j+1Ũj+LS,j+1Ỹj=LG,j+1Nj+1-1Tj+1+LB,j+1Ũj+LS,j+1Ỹj,
where G and Qs 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 P0=(C̃NTΠN-1C̃N)-1 and β0∣k=P0C̃NTΠN-1Ỹm,N-1. β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.


## 3.2.2. 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.Theorem 3.4.
OnN-h≤j≤N-1, γj|k in (3.31) can be computed as follows:
(3.43)γj+1∣k=γj∣k+KjCT(CPjCT+R)-1(yk-N+j-Cβj∣k),
where γN-h∣k=βN-h|k, βj∣k and Pj are obtained from (3.33) and (3.34), respectively, and 𝒦j is given by
(3.44)Kj=Kj-1(I-CT(CPj-1CT+R)-1CPj-1)AT,
with the initial condition 𝒦N-h=PN-h.Proof.
Using Lemma3.2, γj+1∣k in (3.31) for j≥N-h can be represented as
(3.45)γj+1∣k=γj∣k+ΞjNj-1LG,jTCT(CPjCT+R)-1(yk-N+j-Cβj∣k).
If we denote ΞjNj-1LG,jT by 𝒦j, we have only to prove (3.44). Since ΞN-h=LG,N-h, 𝒦N-h is equal to PN-h. Using Lemma 3.2, we can represent 𝒦j in a recursive form as follows:
(3.46)Kj=Kj-1(I-CT(CPjCT+R)-1CPj-1)AT.
This completes the proof.It is noted that the recursive equation (3.43) in Theorem 3.4 is a fixed-point smoother of the state xk-h, which runs with a recursive equation (3.33) in Theorem 3.3. Variables for recursive equations in Sections 3.1 and 3.2 are visualized in Figure 2. Starting from P̂1=0 and α1∣k=0, we compute P̂N and αN∣k recursively in a backward time. From P̂N and αN∣k, we compute P0 and β0∣k, from which we drive the forward recursive equation for βj∣k to get γN∣k=x̂k-h∣k.Figure 2
An iterative form of the LMS RH estimator.


## 4. 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(4.1)xi+1=[1.5400+2δi-0.73790.7379δi]xi+[0.49210.7594]ui+[11]wi,yi=[1+δi1+δi]xi+vi,
where δi is an uncertain model parameter given as(4.2)δi={0.1,200≤i≤220,0,otherwise.
The system noise covariance Q and the measurement noise covariance R are set to 0.012 and 0.0272, 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 fixed-lag 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 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.Figure 3
Estimation errors of LMS RH and Kalman estimators when temporary uncertainties exist.


## 5. Conclusions

In this paper, we proposed a receding horizon (RH) estimator based on the mean-square-error 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 anya 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].



---

*Source: 631759-2012-03-05.xml*