Least-Mean-Square Receding Horizon Estimation

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

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 time 𝑘, the least-mean-square (LMS) RH estimator is to estimate the state 𝑥_𝑘−ℎ at the time 𝑘−ℎ from the inputs 𝑢_𝑘−⋅ and the outputs 𝑦_𝑘−⋅ over the recent finite horizon [𝑘−𝑁,𝑘−1], which can be written as𝑥_{𝑘−ℎ∣𝑘}=𝐄𝑥_𝑘−ℎ∣𝑢_𝑘−𝑁⋯𝑢_𝑘−1,𝑦_𝑘−𝑁⋯𝑦_𝑘−1, (1.1) where ℎ and 𝑁, 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 𝑘−ℎ 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, 𝐄[(𝑥_{𝑘−ℎ∣𝑘}−𝑥_𝑘−ℎ)^𝑇(𝑥_{𝑘−ℎ∣𝑘}−𝑥_𝑘−ℎ)∣𝑢_𝑘−𝑁⋯𝑢_𝑘−1,𝑦_𝑘−𝑁⋯𝑦_𝑘−1]. For ℎ≥2, ℎ=1, and ℎ≤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 𝑁 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 𝑁. 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 (D is 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 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.

2. LMS RH Estimator

Consider a linear discrete-time state space model with an external control input:𝑥_𝑖+1=𝐴𝑥_𝑖+𝐵𝑢_𝑖+G𝑤_𝑖,𝑦_𝑖=𝐶𝑥_𝑖+𝑣_𝑖, (2.1) where 𝑥_𝑖∈ℜ^𝑛, 𝑢_𝑖∈ℜ^𝑙, and 𝑦_𝑖∈ℜ^𝑞 are the state, the input, and the output, respectively, and the system noise 𝑤_𝑖∈ℜ^𝑝 and the measurement noise 𝑣_𝑖∈ℜ^𝑞 are assumed to be zero-mean white Gaussian with𝐄𝑤_𝑖𝑣_𝑖𝑤^𝑇_𝑖𝑣^𝑇_𝑖=𝑆𝑄𝑆_𝑇𝑅,𝑄≥0,𝑅>0, (2.2) where nonzero 𝑆 often happens when I/O models are converted to state space models. It is also assumed that (𝐴,𝐶) of the system is observable. Through this paper, the current time is denoted by 𝑘. For mathematical tractability, the state space model (2.1) is equivalently changed to𝑥_𝑖+1=𝐴_𝑠𝑥_𝑖+𝐵𝑢_𝑖+𝐺𝑤_𝑠,𝑖+𝐺𝑆𝑅⁻¹𝑦_𝑖,𝑦_𝑖=𝐶𝑥_𝑖+𝑣_𝑖, (2.3) where 𝐴_𝑠≜𝐴−𝐺𝑆𝑅⁻¹𝐶 and 𝑤_𝑠,𝑖≜𝑤_𝑖−𝑆𝑅⁻¹𝑣_𝑖. It can be easily shown that 𝑤_𝑠,𝑖 and 𝑣_𝑖 in (2.3) are not correlated. In other words, we have𝐄𝑤_𝑠,𝑖𝑣_𝑖𝑤^𝑇_𝑠,𝑖𝑣^𝑇_𝑖=𝑄−𝑆𝑅⁻¹𝑆^𝑇𝑂𝑂𝑅, (2.4) where 𝑄−𝑆𝑅⁻¹𝑆^𝑇≥0. It is noted that off-diagonal blocks are filled with zeros.

For simple representations, several variables are defined as𝐶_𝑁≜𝐶𝐶𝐴_𝑠𝐶𝐴²_𝑠⋮𝐶𝐴_𝑠^𝑁−1,𝑌_𝑘−1≜𝑦_𝑘−𝑁𝑦_{𝑘−𝑁+1}𝑦_{𝑘−𝑁+2}⋮𝑦_𝑘−1𝐵_𝑁,≜00⋯00𝐶𝐵0⋯00𝐶𝐴_𝑠𝐵𝐶𝐵⋯00⋮⋮⋮⋮⋮𝐶𝐴_𝑠^𝑁−2𝐵𝐶𝐴_𝑠^𝑁−3𝑄𝐵⋯𝐶𝐵0_𝑁,≜𝑄−𝑆𝑅⁻¹𝑆^𝑇⊕⋯⊕^{𝑁0𝑥0𝑓𝑒37}𝑄−𝑆𝑅⁻¹𝑆^𝑇𝑅_𝑁,≜𝑅⊕𝑅⊕⋯^{𝑁0𝑥0𝑓𝑒37}⊕𝑅, (2.5)(2.6)(2.7)(2.8) where ⊕ denotes the direct sum of the matrices. Additionally, 𝑊_𝑘−1,𝑉_𝑘−1,  and 𝑈_𝑘−1 are defined by replacing 𝑦 in (2.5) with 𝑤_𝑠, 𝑣, and 𝑢, respectively. 𝐺_𝑁 and 𝑆_𝑁 are also defined by replacing 𝐵 in 𝐵_𝑁 with 𝐺 and 𝐺𝑆𝑅⁻¹, respectively.

First, we consider the case of 0≤ℎ≤𝑁−1 in (1.1). A prediction problem for ℎ≤−1 will be discussed later on.

By using the defined variables, the state 𝑥_𝑘−ℎ to be estimated is represented in terms of the initial state 𝑥_𝑘−𝑁 on the horizon, inputs, and system noise on the recent finite horizon [𝑘−𝑁,  𝑘−1] as𝑥_𝑘−ℎ=𝐴_𝑠^𝑁−ℎ𝑥_𝑘−𝑁+𝐿_𝑔,𝑁𝑊_𝑘−1+𝐿_𝑏,𝑁𝑈_𝑘−1+𝐿_𝑠,𝑁𝑌_𝑘−1, (2.9) where 𝐿_𝑔,𝑁 is given by𝐿_𝑔,𝑁≜𝐴_𝑠^{𝑁−ℎ−1}𝐺𝐴_𝑠^{𝑁−ℎ−2}𝐺⋯𝐺𝑂⋯𝑂^{ℎ0𝑥0𝑓𝑒37}, (2.10) and 𝐿_𝑏,𝑁 and 𝐿_𝑠,𝑁 are obtained by replacing 𝐺 in 𝐿_𝑔,𝑁 with 𝐵 and 𝐺𝑆𝑅⁻¹, respectively. From (2.9), the conditional expectation 𝐄[𝑥_𝑘−ℎ|𝑈_𝑘−1,𝑌_𝑘−1] in (1.1) can be represented as 𝐄𝑥_𝑘−ℎ∣𝑈_𝑘−1,𝑌_𝑘−1=𝐴_𝑠^𝑁−ℎ𝐄𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1+𝐿_𝑔,𝑁𝐄𝑊_𝑘−1∣𝑈_𝑘−1,𝑌_𝑘−1+𝐿_𝑏,𝑁𝑈_𝑘−1+𝐿_𝑠,𝑁𝑌_𝑘−1. (2.11) Note that if 𝐄[𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1] and 𝐄[𝑊_𝑘−1∣𝑈_𝑘−1,𝑌_𝑘−1] are known, 𝐄[𝑥_𝑘−ℎ∣𝑈_𝑘−1,𝑌_𝑘−1] can be obtained. In order to compute 𝐄[𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1] and 𝐄[𝑊_𝑘−1∣𝑈_𝑘−1,𝑌_𝑘−1], we first try to find the relation among 𝑥_𝑘−𝑁, 𝑊_𝑘−1, 𝑈_𝑘−1, and 𝑌_𝑘−1. For this purpose, the system (2.3) is represented in a batch form on the recent finite horizon [𝑘−𝑁,𝑘−1] as follows:𝑌_𝑘−1=𝐶_𝑁𝑥_𝑘−𝑁+𝐵_𝑁𝑈_𝑘−1+𝐺_𝑁𝑊_𝑘−1+𝑆_𝑁𝑌_𝑘−1+𝑉_𝑘−1, (2.12) from which we obtain𝑊_𝑘−1𝑉_𝑘−1=−𝐺𝐼0_𝑁−𝐶_𝑁𝑊_𝑘−1𝑥_𝑘−𝑁+0𝑌_𝑘−1−𝑆_𝑁𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1. (2.13) We can see from (2.13) that 𝑊_𝑘−1 and 𝑉_𝑘−1 are linearly, more correctly affinely, dependent on 𝑊_𝑘−1 and 𝑥_𝑘−𝑁. The joint probability density function of 𝑊_𝑘−1 and 𝑥_𝑘−𝑁, that is, 𝑓_𝑤𝑥(𝑊_𝑘−1,𝑥_𝑘−𝑁), can be expressed as 𝑘𝑓_𝑤𝑣(𝑊_𝑘−1,𝑉_𝑘−1∣𝑈_𝑘−1,𝑌_𝑘−1)=𝑘𝑓_𝑤𝑣(𝑊_𝑘−1,𝑌_𝑘−1−𝐶_𝑁𝑥_𝑘−𝑁−𝐵_𝑁𝑈_𝑘−1−𝐺_𝑁𝑊_𝑘−1−𝑆_𝑁𝑌_𝑘−1) for an appropriate scaling factor 𝑘. How to choose 𝑘 will be discussed later on. By using the probability density function of 𝑊_𝑘−1 and 𝑉_𝑘−1 given as 𝑓_𝑤𝑣𝑊_𝑘−1,𝑉_𝑘−1=1(2𝜋)^{(𝑝+𝑞)𝑁}det𝑄_𝑁det𝑅_𝑁−1×exp2𝑊_𝑘−1𝑉_𝑘−1^𝑇𝑄_𝑁00𝑅_𝑁⁻¹𝑊_𝑘−1𝑉_𝑘−1, (2.14)𝐄[𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1] can be computed from𝐄𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1=𝑥_𝑘−𝑁𝑓_𝑤𝑥𝑊_𝑘−1,𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1𝑑𝑊_𝑘−1𝑑𝑥_𝑘−𝑁=𝑥_𝑘−𝑁𝑘𝑓_𝑤𝑣𝑊_𝑘−1,𝑌_𝑘−1−𝐶_𝑁𝑥_𝑘−𝑁−𝐵_𝑁𝑈_𝑘−1−𝐺_𝑁𝑊_𝑘−1−𝑆_𝑁𝑌_𝑘−1𝑑𝑊_𝑘−1𝑑𝑥_𝑘−𝑁, (2.15) where a constant 𝑘 is chosen to satisfy the following normalization condition:𝑓_𝑤𝑥𝑊_𝑘−1,𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1𝑑𝑊_𝑘−1𝑑x_𝑘−𝑁=𝑘𝑓_𝑤𝑣𝑊_𝑘−1,𝑌_𝑘−1−𝐶_𝑁𝑥_𝑘−𝑁−𝐵_𝑁𝑈_𝑘−1−𝐺_𝑁𝑊_𝑘−1−𝑆_𝑁𝑌_𝑘−1𝑑𝑊_𝑘−1𝑑𝑥_𝑘−𝑁=1. (2.16) In the same way, 𝐄[𝑊_𝑘−1∣𝑈_𝑘−1,𝑌_𝑘−1] can be obtained. It is noted that 𝑓_𝑤𝑣(𝑊_𝑘−1,𝑌_𝑘−1−𝐶_𝑁𝑥_𝑘−𝑁−𝐵_𝑁𝑈_𝑘−1−𝐺_𝑁𝑊_𝑘−1−𝑆_𝑁𝑌_𝑘−1) is an exponential function with the following exponent:𝑌^𝑇_𝑘−1𝐼−𝑆_𝑁^𝑇𝑈^𝑇_𝑘−1𝑥^𝑇_𝑘−𝑁𝑊^𝑇_𝑘−1𝑅_𝑁⁻¹−𝑅_𝑁⁻¹𝐵_𝑁−𝑅_𝑁⁻¹𝐺_𝑁−𝑅_𝑁⁻¹𝐶_𝑁∗𝐵^𝑇_𝑁𝑅_𝑁⁻¹𝐵_𝑁𝐵^𝑇_𝑁𝑅_𝑁⁻¹𝐺_𝑁𝐵^𝑇_𝑁𝑅_𝑁⁻¹𝐶_𝑁∗∗𝑊_1,1𝑊_1,2∗∗∗𝑊_2,2𝑌^𝑇_𝑘−1𝐼−𝑆_𝑁^𝑇𝑈^𝑇_𝑘−1𝑥^𝑇_𝑘−𝑁𝑊^𝑇_𝑘−1^𝑇, (2.17) where ∗ denotes symmetric parts and 𝑊_1,1, 𝑊_1,2, and 𝑊_2,2 are given by𝑊_1,1=𝐶^𝑇_𝑁𝑅_𝑁⁻¹𝐶_𝑁,𝑊_1,2=𝐶^𝑇_𝑁𝑅_𝑁⁻¹𝐺_𝑁,𝑊_2,2=𝐺^𝑇_𝑁𝑅_𝑁⁻¹+𝑄_𝑁⁻¹. (2.18) By completing the square of the terms of the exponent and recalling the integration of Gaussian functions from −∞ to ∞, we can compute 𝐄[𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1] and 𝐄[𝑊_𝑘−1∣𝑈_𝑘−1,𝑌_𝑘−1]. To begin with, we introduce the following completion of squares:𝑎𝑏^𝑇𝛽𝛼𝛽_𝑇𝛾𝑎𝑏=𝑏+𝛾⁻¹𝛽^𝑇𝑎^𝑇𝛾𝑏+𝛾⁻¹𝛽^𝑇𝑎+𝑎^𝑇𝛼−𝛽𝛾⁻¹𝛽^𝑇𝑎 (2.19) to obtain∫𝑏exp[^𝑎_𝑏]^𝑇_𝛽^𝛼𝛽_𝑇𝛾[^𝑎_𝑏]𝑑𝑏∫exp[^𝑎_𝑏]^𝑇_𝛽^𝛼𝛽_𝑇𝛾[^𝑎_𝑏]𝑑𝑏=−𝛾⁻¹𝛽^𝑇𝑎, (2.20) for some vectors 𝑎 and 𝑏, 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:𝑌𝑎⟵^𝑇𝑘−1𝐼−𝑆_𝑁^𝑇𝑈^𝑇_𝑘−1𝑥,𝑏⟵^{𝑇𝑘−𝑁}𝑊^𝑇_𝑘−1𝑅,𝛼⟵^𝑁−1−𝑅_𝑁⁻¹𝐵_𝑁∗𝐵^𝑇_𝑁𝑅_𝑁⁻¹𝐵_𝑁,𝛽⟵−𝑅_𝑁⁻¹𝐺_𝑁−𝑅_𝑁⁻¹𝐶_𝑁𝐵^𝑇_𝑁𝑅_𝑁⁻¹𝐺_𝑁𝐵^𝑇_𝑁𝑅_𝑁⁻¹𝐶_𝑁𝑊,𝛾⟵^1,1𝑊_1,2∗𝑊_2,2, (2.21) we have𝐄𝑥_𝑘−𝑁∣𝑈_𝑘−1,𝑌_𝑘−1𝐄𝑊_𝑘−1∣𝑈_𝑘−1,𝑌_𝑘−1=𝑊_1,1𝑊_1,2𝑊^𝑇_1,2𝑊_2,2⁻¹𝐶_𝑁𝐺_𝑁×𝑅_𝑁⁻¹𝐼−𝑆_𝑁𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1, (2.22) where 𝑊_1,1, 𝑊_1,2, and 𝑊_2,2 are given by (2.18). It is noted that if (𝐴,𝐶) of the system (2.1) is observable, 𝑊_1,1 is positive definite and 𝑊_2,2−𝑊^𝑇_1,2𝑊⁻¹_1,1𝑊_1,2 is also positive definite, which implies that the block matrix including 𝑊_1,1, 𝑊_1,2, and 𝑊_2,2 is guaranteed to be positive definite and hence nonsingular. Substituting (2.22) into (2.11), we can obtain 𝐄[𝑥_𝑘−ℎ∣𝑈_𝑘−1,𝑌_𝑘−1]. Through a long and tedious algebraic calculation, we have the covariance matrix of the mean-square-error 𝐄[(𝑥_{𝑘−ℎ∣𝑘}−𝑥_𝑘−ℎ)(𝑥_{𝑘−ℎ∣𝑘}−𝑥_𝑘−ℎ)^T] as follows:𝑃=Ξ_𝑁𝑊_1,1𝑊_1,2𝑊^𝑇_1,2𝑊_2,2⁻¹Ξ^𝑇_𝑁, (2.23) where Ξ_𝑁 is defined byΞ_𝑗≜𝐴_𝑠^𝑁−ℎ𝐴_𝑠^{𝑁−ℎ−1}𝐺𝐴_𝑠^{𝑁−ℎ−2}𝐺⋯𝐺𝑂𝑂⋯^{𝑗−𝑁+ℎ0𝑥0𝑓𝑒37}𝑂, (2.24) for 𝑁−ℎ≤𝑗≤𝑁.

What we have done until now can be summarized in the following theorem.

Theorem 2.1. Suppose that (𝐴,𝐶) of the system (2.1) is observable. Then, the LMS RH estimator (1.1) is given by 𝑥_{𝑘−ℎ|𝑘}=Ξ_𝑁𝑊_1,1𝑊_1,2𝑊^𝑇_1,2𝑊_2,2⁻¹𝐶_𝑁𝐺_𝑁𝑅_𝑁⁻¹𝐼−𝑆_𝑁𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1+𝐿_𝑏,𝑁𝑈_𝑘−1+𝐿_𝑠,𝑁𝑌_𝑘−1, (2.25) where 𝑥_{𝑘−ℎ∣𝑘} denotes 𝐄[𝑥_𝑘−ℎ∣𝑈_𝑘−1,𝑌_𝑘−1] and 𝑊_1,1, 𝑊_1,2, and 𝑊_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 𝐴. 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, 𝑄 and 𝐵_𝑁 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 [𝑘−𝑁,𝑘−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, 𝑥_{𝑘−ℎ|𝑘}=𝑥_𝑘−ℎ at the fixed-lag (ℎ≥1) or the current (ℎ=0) time.

Proof. The LMS RH estimator (2.25) can be rearranged in terms of 𝑌_𝑘−1 and 𝑈_𝑘−1, that is, 𝐻𝑌_𝑘−1+𝐿𝑈_𝑘−1 for some gain matrices 𝐻 and 𝐿. Substituting 𝑌_𝑘−1 in (2.12) into 𝐻𝑌_𝑘−1+𝐿𝑈_𝑘−1, using (2.9), and removing all noise, we have 𝐻𝑌_𝑘−1+𝐿𝑈_𝑘−1=𝐻𝐶_𝑁−𝐴_𝑠^𝑁−ℎ𝑥_𝑘−𝑁+𝑥_𝑘−ℎ+𝐻𝐵_𝑁−𝐿_𝑏,𝑁+𝐿𝑈_𝑘−1+𝐻𝑆_𝑁−𝐿_𝑠,𝑁𝑌_𝑘−1. (2.26) It can be easily seen that 𝐻𝐶_𝑁=𝐴_𝑠^𝑁−ℎ and terms associated with inputs 𝑈_𝑘−1 and outputs 𝑌_𝑘−1 become zero. It follows then that we have 𝐻𝑌_𝑘−1+𝐿𝑈_𝑘−1=𝑥_𝑘−ℎ. 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 ℎ≤−1 can be easily solved from the previous results for ℎ=0. Since 𝑥_𝑘−ℎ is represented as𝑥_𝑘−ℎ=𝐴^−ℎ𝑥_𝑘+^−ℎ−1_𝑖=0𝐴^𝑖𝐵𝑢_𝑘+𝑖+^−ℎ−1_𝑖=0𝐴^𝑖𝐺𝑤_𝑘+𝑖, (2.27) we have𝐄𝑥_𝑘−ℎ∣𝑈_𝑘−1,𝑌_𝑘−1=𝐴^−ℎ𝐄𝑥_𝑘∣𝑈_𝑘−1,𝑌_𝑘−1+^−ℎ−1_𝑖=0𝐴^𝑖𝐵𝑢_𝑘+𝑖, (2.28) which means that 𝐄[𝑥_𝑘−ℎ∣𝑈_𝑘−1,𝑌_𝑘−1] can be obtained from 𝐄[𝑥_𝑘∣𝑈_𝑘−1,𝑌_𝑘−1]. Setting ℎ in (2.25) to zero, we can compute 𝐄[𝑥_𝑘∣𝑈_𝑘−1,𝑌_𝑘−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:𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁⁻¹=𝑊_1,1−𝑊_1,2𝑊⁻¹_2,2𝑊^𝑇_1,2⁻¹,𝐶^𝑇_𝑁Π_𝑁⁻¹=𝐶^𝑇_𝑁𝑅_𝑁⁻¹−𝑊_1,2𝑊⁻¹_2,2𝐺^𝑇_𝑁𝑅_𝑁⁻¹, (3.1) we obtain𝐶^𝑇_𝑁Π_𝑁⁻¹0+𝐶^𝑇_𝑁𝐺^𝑇_𝑁𝑅_𝑁⁻¹=2𝐶^𝑇_𝑁𝑅_𝑁⁻¹−𝑊_1,2𝑊⁻¹_2,2𝐺^𝑇_𝑁𝑅_𝑁⁻¹𝐺^𝑇_𝑁𝑅_𝑁⁻¹=2𝐼−𝑊_1,2𝑊⁻¹_2,2𝐶0𝐼𝐺^𝑇𝑁𝑇_𝑁𝑅_𝑁⁻¹,𝑊_1,1+𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁𝑊_2,2𝑊^𝑇_1,2𝑊_2,2⁻¹2𝐼−𝑊_1,2𝑊⁻¹_2,2=10𝐼2𝐼12𝑊_1,2𝑊⁻¹_2,2𝑊0𝐼^1,1+𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁𝑊_1,2𝑊^𝑇_1,2𝑊_2,2⁻¹=𝑊_1,1𝑊_1,2𝑊^𝑇_1,2𝑊_2,2⁻¹. (3.2) Note that nonsingularity of 𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁 is guaranteed for the observability of (𝐴,𝐶). From (3.2), it can be easily seen that the LMS RH (2.25) can be represented as𝑥_{𝑘−ℎ|𝑘}=Ξ_𝑁𝑊_1,1+𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁𝑊_1,2𝑊^𝑇_1,2𝑊_2,2⁻¹×𝐶^𝑇_𝑁Π_𝑁⁻¹𝐼−𝑆_𝑁𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−10+𝐶^𝑇_𝑁𝐺^𝑇_𝑁𝑅_𝑁⁻¹𝐼−𝑆_𝑁𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1+𝐿_𝑠,𝑁𝑌_𝑘−1+𝐿_𝑏,𝑁𝑈_𝑘−1. (3.3)

Now, we consider two recursive equations for the batch form (3.3). One is for obtaining 𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁 and 𝐶^𝑇_𝑁Π_𝑁⁻¹((𝐼−𝑆_𝑁)𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1). The other is for (3.3) given 𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁 and 𝐶^𝑇_𝑁Π_𝑁⁻¹((𝐼−𝑆_𝑁)𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−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 obtain 𝐶^𝑇_𝑁Π⁻¹_N𝐶_𝑁 and 𝐶^𝑇_𝑁Π_𝑁⁻¹((𝐼−𝑆_𝑁)𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1) in (3.3) will be discussed. These values will be computed in a backward time. (𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁)⁻¹𝐶^𝑇_𝑁Π_𝑁⁻¹((𝐼−𝑆_𝑁)𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1) and (𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁)⁻¹ will be denoted by 𝛽_0|𝑘 and 𝑃₀, respectively, for consistency with the next section. The following theorem provides the main result.

Theorem 3.1. 𝑃₀=(𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁)⁻¹ and 𝛽_0|𝑘=𝑃₀𝐶^𝑇_𝑁Π_𝑁⁻¹((𝐼−𝑆_𝑁)𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1) can be computed recursively as follows: 𝑃₀=𝐶^𝑇𝑅⁻¹𝐶+𝑃_𝑁⁻¹,𝛽_0|𝑘=𝑃₀𝐶^𝑇𝑅⁻¹𝑦_𝑘−𝑁+𝛼_𝑁|𝑘, (3.4)(3.5) where 𝑃_𝑁 and 𝛼_𝑁∣𝑘 are obtained from 𝑃_𝑗+1=𝐴^𝑇_𝑠𝐶^𝑇𝑅⁻¹𝐶𝐴_𝑠+𝐴^𝑇_𝑠𝑃_𝑗𝐴_𝑠−𝐴^𝑇_𝑠𝐶^𝑇𝑅⁻¹𝐶+𝑃_𝑗𝐺×𝑄_𝑠⁻¹+𝐺^𝑇𝐶^𝑇𝑅⁻¹𝐶+𝑃_𝑗𝐺⁻¹𝐺^𝑇𝐶^𝑇𝑅⁻¹𝐶+𝑃_𝑗𝐴_𝑠,𝛼_{𝑗+1∣𝑘}=𝐴^𝑇_𝑠−𝐴^𝑇_𝑠𝐶^𝑇𝑅⁻¹𝐶+𝑃_𝑗𝐺𝑄_𝑠⁻¹+𝐺^𝑇𝐶^𝑇𝑅⁻¹𝐶+𝑃_𝑗𝐺⁻¹𝐺^𝑇×𝛼_𝑗∣𝑘+𝐶^𝑇𝑅⁻¹𝑦_𝑘−𝑗−𝑃_𝑗+𝐶^𝑇𝑅⁻¹𝐶𝐵𝑢_{𝑘−𝑗−1}+𝐺𝑆𝑅⁻¹𝑦_{𝑘−𝑗−1}, (3.6)(3.7) for 1≤𝑗≤𝑁−1, 𝑄_𝑠=𝑄−𝑆𝑅⁻¹𝑆^𝑇,  and 𝑃₁ and 𝛼_1∣𝑘 are zero matrices with appropriate dimensions.

Proof. Before going into a main proof, we introduce some variables 𝐶_𝑗, 𝑁_𝑗, Π_𝑗, and 𝑃_𝑗 as 𝐶_𝑗≜𝐶𝐴_𝑠⋮𝐶𝐴_𝑠^𝑗−3𝐶𝐴_𝑠^𝑗−2𝐶𝐴_𝑠^𝑗−1=𝐶𝐶_𝑗−1𝑁𝐴,_𝑗≜C𝐺0⋯0𝐶𝐴_𝑠𝐺𝐶𝐺⋯0⋮⋮⋮⋮𝐶𝐴_𝑠^𝑗−2𝐺𝐶𝐴_𝑠^𝑗−3=𝐺⋯𝐶𝐺𝐶𝐺0𝐶_𝑗−1𝐺𝑁_𝑗−1Π_𝑗,≜𝑁_𝑗𝑄_𝑗−1𝑁^𝑇_𝑗+𝑅_𝑗−1∈ℜ^{𝑞(𝑗−1)×𝑞(𝑗−1)}𝑃_𝑗,≜𝐶^𝑇_𝑗Π_𝑗⁻¹𝐶_𝑗, (3.8)(3.9)(3.10)(3.11) for 2≤𝑗≤𝑁. In terms of Π_𝑗 in (3.10) and 𝑃_𝑗 in (3.11), 𝛽_0∣𝑘 and 𝑃₀ can be represented as 𝛽_0∣𝑘=𝑃₀𝐶^𝑇𝑅⁻¹𝑦_𝑘−𝑁+𝐶^𝑇_𝑁Π_𝑁⁻¹𝐼^𝑜−𝑆^𝑜_𝑁𝑌_𝑁−𝐵^𝑜_𝑁𝑈_𝑁=𝑃₀𝐶^𝑇𝑅⁻¹𝑦_𝑘−𝑁+𝛼_𝑁∣𝑘,𝑃₀=𝐶^𝑇𝑅⁻¹𝐶+𝑃_𝑁⁻¹, (3.12) where 𝐼^𝑜, 𝑆^𝑜_𝑗, and 𝐵^𝑜_𝑗 are obtained by removing the first row blocks of 𝐼, 𝑆_𝑗, and 𝐵_𝑗, respectively, and 𝑌_𝑗, 𝑈_𝑗, and 𝛼_𝑗|𝑘 are given by 𝑌_𝑗≜𝑦^𝑇_𝑘−𝑗⋯𝑦^𝑇_𝑘−1^𝑇,𝑈_𝑗≜𝑢^𝑇_𝑘−𝑗⋯𝑢^𝑇_𝑘−1^𝑇,𝛼_𝑗∣𝑘≜𝐶^𝑇_𝑗Π_𝑗⁻¹𝐼^𝑜−𝑆^𝑜_𝑗𝑌_𝑗−𝐵^𝑜_𝑗𝑈_𝑗, (3.13)(3.14) for 2≤𝑗≤𝑁. In order to obtain 𝑃₀ and 𝛽_0|𝑘, we have only to know 𝑃_𝑁 and 𝛼_𝑁|𝑘. Now, in order to get 𝑃_𝑁 and 𝛼_𝑁|𝑘, we try to find recursive equations for 𝑃_𝑗 and 𝛼_𝑗|𝑘 on 1≤𝑗≤𝑁. By using recursions in (3.8) and (3.9), we have the following equality: Π⁻¹_𝑗+1=Δ_𝑗⁻¹−Δ_𝑗⁻¹𝐶𝐶_𝑗𝐺𝑄_𝑠⁻¹+𝐺^𝑇𝐶𝐶_𝑗^𝑇Δ_𝑗⁻¹𝐶𝐶_𝑗𝐺⁻¹×𝐺^𝑇𝐶𝐶_𝑗^𝑇Δ_𝑗⁻¹, (3.15) where Δ_𝑗 is given by Δ_𝑗≜0Π𝑅0_𝑗. (3.16) Pre- and postmultiplying (3.15) by 𝐶^𝑇_𝑗+1=𝐴^𝑇_𝑠𝐶𝐶_𝑗^𝑇,𝐶_𝑗+1=𝐶𝐶_𝑗𝐴_𝑠, (3.17) respectively, we have (3.6). Pre- and postmultiplying (3.15) by 𝐶^𝑇_𝑗+1=𝐴^𝑇_𝑠𝐶𝐶_𝑗^𝑇,𝐼^𝑜−𝑆^𝑜_𝑗+1𝑌_𝑗+1−𝐵^𝑜_𝑗+1𝑈_𝑗+1=𝑦_𝑘−𝑗𝐼^𝑜−𝑆^𝑜_𝑗𝑌_𝑗−𝐵^𝑜_𝑗𝑈_𝑗−𝐶𝐶_𝑗𝐵𝑢_{𝑘−𝑗−1}+𝐺𝑆𝑅⁻¹𝑦_{𝑘−𝑗−1}, (3.18)(3.19) respectively, we have (3.7). Note that (3.6) and (3.7) hold for 𝑖≥2. From (3.11) and (3.14), 𝑃₂ and 𝛼_2∣𝑘 can be written as 𝑃₂=𝐶^𝑇₂Π₂⁻¹𝐶₂=𝐴^𝑇_𝑠𝐶^𝑇𝐶𝐺𝑄_𝑠𝐺^𝑇𝐶^𝑇+𝑅⁻¹𝐶𝐴^𝑇_𝑠,𝛼_2∣𝑘=𝐶^𝑇₂Π₂⁻¹𝐼^𝑜−𝑆^𝑜₂𝑌₂−𝐵^𝑜₂𝑈₂=𝐴^𝑇_𝑠𝐶^𝑇𝐶𝐺𝑄_𝑠𝐺^𝑇𝐶^𝑇+𝑅⁻¹𝑦_𝑘−1−𝐶𝐺𝑆𝑅⁻¹𝑦_𝑘−2−𝐶𝐵𝑢_𝑘−2. (3.20)(3.21) If 𝑃₁ and 𝛼_1∣𝑘 are set to zero matrices with appropriate dimensions, 𝑃₂ in (3.20) and 𝛼_2∣𝑘 in (3.21) can be calculated from (3.6) and (3.7). Thus, we can say that (3.6) and (3.7) hold for 𝑖≥1 and are initiated with 𝑃₁=0 and 𝛼_1∣𝑘=0. After obtaining 𝑃_𝑁 from (3.6), we can calculate 𝑃₀ from (3.4). 𝛽_0∣𝑘 in (3.5) can be obtained from 𝛼_𝑁|𝑘 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 𝑃₀=(𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁)⁻¹ and 𝛽_0∣𝑘=𝑃₀𝐶^𝑇_𝑁Π_𝑁⁻¹((𝐼−𝑆_𝑁)𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1) are given. 𝑃₀ and 𝛽_0|𝑘 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.

𝐿_𝐺,𝑖, 𝑀_𝑖, and 𝑁_𝑖 are defined as𝐿_𝐺,𝑖≜𝐴^𝑖_𝑠𝐴_𝑠^𝑖−1𝐺𝐴_𝑠^𝑖−2𝐺⋯𝐴_𝑠𝑀𝐺𝐺,_𝑖≜𝐶^𝑇_𝑖𝑅_𝑖⁻¹𝐶_𝑖+𝑃₀⁻¹𝐶^𝑇_𝑖𝑅_𝑖⁻¹𝐺^𝑜_𝑖𝐺_𝑖^𝑜𝑇𝑅_𝑖⁻¹𝐶_𝑖𝐺_𝑖^𝑜𝑇𝑅_𝑖⁻¹𝐺^𝑜_𝑖+𝑄⁻¹_𝑖−1𝑁_𝑖,≜𝑀_𝑖00𝑄_𝑠⁻¹=𝐶^𝑇_𝑖𝑅_𝑖⁻¹𝐶_𝑖+𝑃₀⁻¹𝐶^𝑇_𝑖𝑅_𝑖⁻¹𝐺_𝑖𝐺^𝑇_𝑖𝑅_𝑖⁻¹𝐶_𝑖𝐺^𝑇_𝑖𝑅_𝑖⁻¹𝐺_𝑖+𝑄_𝑖⁻¹, (3.22)(3.23)(3.24) where 2≤𝑖≤𝑁 and 𝐺^𝑜_𝑖 is the matrix obtained by removing the last zero column block from 𝐺_𝑖, that is, 𝐺_𝑖=𝐺^𝑜_𝑖0. In particular, 𝐿_𝐺,1, 𝑀₁, and 𝑁₁ are defined as [𝐴_𝑠𝐺], 𝐶^𝑇𝑅⁻¹𝐶+𝑃₀⁻¹ and 𝑀₁⊕𝑄_𝑠⁻¹. The following lemma shows how variables 𝐿_𝐺,𝑖, 𝑁_𝑖, and 𝑀_𝑖 are related to one another.

Lemma 3.2. The following relation is satisfied: 𝑁_𝑖+𝐿^𝑇_𝐺,𝑖𝐶^𝑇𝑅⁻¹𝐶𝐿_𝐺,𝑖=𝑀_𝑖+1, (3.25) where 𝐿_𝐺,𝑖,  𝑁_𝑖,  and 𝑀_𝑖 are defined in (3.22)–(3.24), respectively.

Proof. If Γ_𝑖 is defined as Γ_𝑖≜𝐴_𝑠^𝑖−1𝐺𝐴_𝑠^𝑖−2𝐺⋯𝐴_𝑠𝐺𝐺, (3.26)𝐿_𝐺,𝑖 can be represented as 𝐿_𝐺,𝑖=𝐴^𝑖_𝑠Γ_𝑖, (3.27) from which we have 𝐿^𝑇_𝐺,𝑖𝐶^𝑇𝑅⁻¹𝐶𝐿_𝐺,𝑖=𝐴_𝑠^𝑇𝑖𝐶^𝑇𝑅⁻¹𝐶𝐴^𝑖_𝑠𝐴_𝑠^𝑇𝑖𝐶^𝑇𝑅⁻¹𝐶Γ_𝑖Γ^𝑇_𝑖𝐶^𝑇𝑅⁻¹𝐶𝐴^𝑖_𝑠Γ^𝑇_𝑖𝐶^𝑇𝑅⁻¹𝐶Γ_𝑖. (3.28)
The four block elements in 𝑀_𝑖+1 can be expressed recursively as 𝐶^𝑇_𝑖+1𝑅⁻¹_𝑖+1𝐶_𝑖+1=𝐶^𝑇_𝑖𝑅_𝑖⁻¹𝐶_𝑖+𝐴_𝑠^𝑇𝑖𝐶^𝑇𝑅⁻¹𝐶𝐴^𝑖_𝑠,𝐶^𝑇_𝑖+1𝑅⁻¹_𝑖+1𝐺^𝑜_𝑖+1=𝐶^𝑇_𝑖𝑅_𝑖⁻¹𝐺^𝑜_𝑖0+𝐴_𝑠^𝑇𝑖𝐶^𝑇𝑅⁻¹𝐶Γ_𝑖,𝐺^𝑜𝑇_𝑖+1𝑅⁻¹_𝑖+1𝐺^𝑜_𝑖+1=𝐺_𝑖^𝑜𝑇𝑅_𝑖⁻¹𝐺^𝑜_𝑖000+Γ^𝑇_𝑖𝐶^𝑇𝑅⁻¹𝐶Γ_𝑖. (3.29) Using (3.28) and (3.29), we have 𝑁_𝑖+𝐿^𝑇_𝐺,𝑖𝐶^𝑇𝑅⁻¹𝐶𝐿_𝐺,𝑖=𝑀_𝑖+1. 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:𝛽_𝑗∣𝑘=𝐿_𝐺,𝑗𝑁_𝑗⁻¹𝐶^𝑇_𝑁Π_𝑁⁻¹𝑌_{𝑚,𝑁−1}0+𝐶^𝑇_𝑗𝐺^𝑇_𝑗𝑅_𝑗⁻¹𝑌_{𝑚,𝑗−1}+𝐿_𝐵,𝑗𝑈_𝑗−1+𝐿_𝑆,𝑗𝑌_𝑗−1, (3.30) for 1≤𝑗≤𝑁𝛾_𝑗∣𝑘=Ξ_𝑗𝑁_𝑗⁻¹𝐶^𝑇_𝑁Π_𝑁⁻¹𝑌_{𝑚,𝑁−1}0+𝐶^𝑇_𝑗𝐺^𝑇_𝑗𝑅_𝑗⁻¹𝑌_{𝑚,𝑗−1}+𝐿_𝑏,𝑗𝑈_𝑗−1+𝐿_𝑠,𝑗𝑌_𝑗−1, (3.31) for 𝑁−ℎ≤𝑗≤𝑁, where 𝐿_𝑏,𝑗 and 𝐿_𝑠,𝑗 are defined in a form (2.24), and 𝑌_𝑗−1, 𝑈_𝑗−1, 𝐿_𝐵,𝑗, 𝐿_𝑆,𝑗, and 𝑌_{𝑚,𝑗−1} are given by𝑌_𝑗−1≜𝑦^𝑇_𝑘−𝑁⋯𝑦^𝑇_{𝑘−𝑁+𝑗−1}^𝑇,𝑈_𝑗−1≜𝑢^𝑇_𝑘−𝑁⋯𝑢^𝑇_{𝑘−𝑁+𝑗−1}^𝑇,𝐿_𝐵,𝑗≜𝐴_𝑠^𝑗−1𝐵𝐴_𝑠^𝑗−2,𝐿𝐵⋯𝐵_𝑆,𝑗≜𝐴_𝑠^𝑗−1𝐺𝑆R⁻¹𝐴_𝑠^𝑗−2𝐺𝑆𝑅⁻¹⋯𝐺𝑆𝑅⁻¹,𝑌_{𝑚,𝑗−1}≜𝐼−𝑆_𝑗𝑌_𝑗−1−𝐵_𝑗𝑈_𝑗−1. (3.32) Note that 𝛽_0∣𝑘=𝐶^𝑇_𝑁Π_𝑁⁻¹((𝐼−𝑆_𝑁)𝑌_𝑘−1−𝐵_𝑁𝑈_𝑘−1), 𝛾_{𝑁−ℎ∣𝑘}=𝛽_{𝑁−ℎ∣𝑘}, and 𝑥_{𝑘−ℎ∣𝑘}=𝛾_𝑁∣𝑘. However, 𝛽_𝑗∣𝑘 is recursively computed from 𝑗=0 to 𝑗=𝑁 and 𝛾_𝑗∣𝑘 from 𝑗=𝑁−ℎ to 𝑗=𝑁 with the initial value 𝛾_{𝑁−ℎ∣𝑘}=𝛽_{𝑁−ℎ∣𝑘}. Even though 𝛽_𝑗∣𝑘 does not look useful on [𝑘−ℎ,𝑘], it is still used to compute 𝛾_𝑗∣𝑘 on that horizon. Now, we try to find out recursive equations for 𝛽_𝑗∣𝑘 and 𝛾_𝑗∣𝑘 in what follows.

3.2.1. Recursive Equation for 𝛽_𝑗|𝑘 on 0≤𝑗≤𝑁

Using Lemma 3.2, we will obtain a recursive equation for 𝛽_𝑗∣𝑘, which is introduced in the following theorem.

Theorem 3.3. On 0≤𝑗≤𝑁, 𝛽_𝑗∣𝑘 in (3.30) can be computed as follows: 𝛽_{𝑗+1∣𝑘}=𝐴_𝑠𝛽_𝑗∣𝑘+𝐵𝑢_{𝑘−𝑁+𝑗}+𝐺𝑆𝑅⁻¹𝑦_{𝑘−𝑁+𝑗}+𝐴_𝑠𝑃_𝑗𝐶^𝑇𝑅+𝐶𝑃_𝑗𝐶^𝑇−1𝑦_{𝑘−𝑁+𝑗}−𝐶𝛽_𝑗|𝑘, (3.33) where 𝑃_𝑗 is given by 𝑃_𝑗+1=𝐴_𝑠𝑃_𝑗𝐴^𝑇_𝑠+𝐺𝑄_𝑠𝐺^𝑇−𝐴_𝑠𝑃_𝑗𝐶^𝑇𝑅+𝐶𝑃_𝑗𝐶^𝑇−1𝐶𝑃_𝑗𝐴^𝑇_𝑠, (3.34) and initial conditions are set to 𝑃₀=(𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁)⁻¹ and 𝛽_0|𝑘=𝑃₀𝐶^𝑇_𝑁Π_𝑁⁻¹𝑌_{𝑚,𝑁−1} computed in Section 3.1.

Proof. First, we will obtain a closed-form of 𝑃_𝑗 in (3.34) and then, using the closed-form of 𝑃_𝑗, we show that 𝛽_𝑗∣𝑘 in (3.30) can be computed recursively from (3.33).
Using the defined variables (3.22) and (3.24), we assume that 𝑃_𝑗 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 𝑗=1. Given the initial value 𝑃₀, we transform 𝑃₁ to the form (3.35) by (3.34) as follows: 𝑃₁=𝐴_𝑠𝐺𝑃₀⁻¹+𝐶^𝑇𝑅⁻¹𝐶00𝑄_𝑠⁻¹⁻¹𝐴_𝑠𝐺^𝑇. (3.36) Equation (3.38) can be written in terms of 𝐿_𝐺,𝑖 and 𝑁_𝑖 as 𝑃₁=𝐿_𝐺,1𝑁₁⁻¹𝐿^𝑇_𝐺,1. Now, we check 𝑃_𝑗+1 under the assumption that 𝑃_𝑗=𝐿_𝐺,𝑗𝑁_𝑗⁻¹𝐿^𝑇_𝐺,𝑗. From (3.34), we have 𝑃_𝑗+1=𝐴_𝑠𝐿_𝐺,𝑗𝑀⁻¹_𝑗+1𝐿^𝑇_𝐺,𝑗𝐴^𝑇_𝑠+𝐺𝑄_𝑠𝐺^𝑇=𝐿_𝐺,𝑗+1𝑁⁻¹_𝑗+1𝐿^𝑇_𝐺,𝑗+1. (3.37) Thus, we can see that 𝑃_𝑗 in (3.34) can be represented as the form (3.35) in terms of L_𝐺,𝑗 and 𝑁_𝑗. Using this result, we are in a position to show that 𝛽_𝑗∣𝑘 in (3.30) can be computed recursively from (3.33). We can rewrite 𝛽_𝑗|𝑘 in (3.30) as 𝛽_𝑗∣𝑘=𝐿_𝐺,𝑗𝑁_𝑗⁻¹𝑇_𝑗+𝐿_𝐵,𝑗𝑈_𝑗−1+𝐿_𝑆,𝑗𝑌_𝑗−1, (3.38) where 𝑇_𝑗 is given by 𝑇_𝑗≜𝑃₀⁻¹0𝛽_0|𝑘+𝐶^𝑇_𝑗𝐺^𝑇_𝑗𝑅_𝑗⁻¹𝑌_{𝑚,𝑗−1}. (3.39) As in a derivation of 𝑃_𝑗 in a batch form, we show by an induction method that 𝛽_𝑗|𝑘 in (3.38) is equivalent to the one in (3.33). First, we show that 𝛽_1|𝑘 can be obtained from 𝑃₀=(𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁)⁻¹ and 𝛽_0∣𝑘=𝑃₀𝐶^𝑇_𝑁Π_𝑁⁻¹𝑌_{𝑚,𝑁−1}: 𝛽_1∣𝑘=𝐴_𝑠𝑃₀⁻¹+𝐶^𝑇𝑅⁻¹𝐶⁻¹𝑃₀⁻¹𝛽_0∣𝑘+𝐶^𝑇𝑅⁻¹𝑦_𝑘−𝑁+𝐵𝑢_𝑘−𝑁+𝐺𝑆𝑅⁻¹𝑦_𝑘−𝑁, (3.40) where, in terms of 𝐿_𝐺,𝑖, 𝑁_𝑖, and 𝑇_𝑖, 𝛽_1|𝑘 can be written as 𝛽_1|𝑘=𝐿_𝐺,1𝑁₁⁻¹𝑇₁+𝐿_𝐵,1𝑈₀+𝐿_𝑆,1𝑌₀. Next, we show that 𝛽_{𝑗+1∣𝑘} in the form (3.38) can be obtained from 𝛽_𝑗∣𝑘 of the batch form, that is, 𝛽_𝑗|𝑘=𝐿_𝐺,𝑗𝑁_𝑗⁻¹𝑇_𝑗+𝐿_𝐵,𝑗𝑈_𝑗−1+𝐿_𝑆,𝑗𝑌_𝑗−1. First, note that we have the following relations: 𝐴_𝑠−𝐴_𝑠𝑃_𝑗𝐶^𝑇𝑅+𝐶𝑃_𝑗𝐶^𝑇−1𝐶𝐿_𝐺,𝑗𝑁_𝑗⁻¹𝑇_𝑗=𝐴_𝑠𝐿_𝐺,𝑗𝑀⁻¹_𝑗+1𝑇_𝑗,𝐴_𝑠𝑃_𝑗𝐶^𝑇𝑅+𝐶𝑃_𝑗𝐶^𝑇−1𝑦_{𝑘−𝑁+𝑗}=𝐴_𝑠𝐿_𝐺,𝑗𝑀⁻¹_𝑗+1𝐿^𝑇_𝐺,𝑗𝐶^𝑇𝑅⁻¹𝑦_{𝑘−𝑁+𝑗}. (3.41) Substituting (3.41) into (3.33) yields 𝛽_{𝑗+1∣𝑘}=𝐴_𝑠𝐿_𝐺,𝑗𝑀⁻¹_𝑗+1𝑇_𝑗+𝐴_𝑠𝐿_𝐺,𝑗𝑀⁻¹_𝑗+1𝐿^𝑇_𝐺,𝑗𝐶^𝑇𝑅⁻¹𝑦_{𝑘−𝑁+𝑗}−𝐴_𝑠𝐿_𝐺,𝑗𝑀⁻¹_𝑗+1𝐿^𝑇_𝐺,𝑗𝐶^𝑇𝑅⁻¹𝐶𝐿_𝐵,𝑗𝑈_𝑗−1+𝐿_𝑆,𝑗𝑌_𝑗−1+𝐿_𝐵,𝑗+1𝑈_𝑗+𝐿_𝑆,𝑗+1𝑌_𝑗,=𝐴_𝑠𝐿_𝐺,𝑗𝐺𝑀⁻¹_𝑗+100𝑄_𝑠𝑃₀⁻¹0𝛽_0∣𝑘+𝐶^𝑇_𝑗+1𝐺^𝑜_𝑗+10^𝑇𝑅⁻¹_𝑗+1𝑌_𝑚,𝑗+𝐿_𝐵,𝑗+1𝑈_j+𝐿_𝑆,𝑗+1𝑌_𝑗=𝐿_𝐺,𝑗+1𝑁⁻¹_𝑗+1𝑇_𝑗+1+𝐿_𝐵,𝑗+1𝑈_𝑗+𝐿_𝑆,𝑗+1𝑌_𝑗, (3.42) where 𝐺 and 𝑄_𝑠 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 𝑃₀=(𝐶^𝑇_𝑁Π_𝑁⁻¹𝐶_𝑁)⁻¹ and 𝛽_0∣𝑘=𝑃₀𝐶^𝑇_𝑁Π_𝑁⁻¹𝑌_{𝑚,𝑁−1}. 𝛽_𝑗|𝑘 on 1≤𝑗≤𝑁 provides initial values and inputs for a recursive equation of 𝛾_𝑗∣𝑘 (3.31), which will be investigated in what follows.

3.2.2. Recursive Equation for 𝛾_𝑗∣𝑘 on 𝑁−ℎ≤𝑗≤𝑁

Here, we discuss the recursive equation for 𝛾_𝑗∣𝑘 on 𝑁−ℎ≤𝑗≤𝑁. As mentioned before, the recursive equation for 𝛾_𝑗∣𝑘 starts from 𝑗=𝑁−ℎ with the initial value 𝛾_{𝑁−ℎ∣𝑘}=𝛽_{𝑁−ℎ∣𝑘}. On 𝑁−ℎ≤𝑗≤𝑁, 𝛾_𝑗∣𝑘 can be recursively computed with the help of 𝛽_𝑗∣𝑘. However, 𝛾_𝑁∣𝑘=𝑥_{𝑘−ℎ∣𝑘} is what we want to find out finally. The recursive equation for 𝛾_𝑗∣𝑘 is given in the following theorem.

Theorem 3.4. On 𝑁−ℎ≤𝑗≤𝑁−1, 𝛾_𝑗|𝑘 in (3.31) can be computed as follows: 𝛾_{𝑗+1∣𝑘}=𝛾_𝑗∣𝑘+𝒦_𝑗𝐶^𝑇𝐶𝑃_𝑗𝐶^𝑇+𝑅⁻¹𝑦_{𝑘−𝑁+𝑗}−𝐶𝛽_𝑗∣𝑘, (3.43) where 𝛾_{𝑁−ℎ∣𝑘}=𝛽_{𝑁−ℎ|𝑘}, 𝛽_𝑗∣𝑘 and 𝑃_𝑗 are obtained from (3.33) and (3.34), respectively, and 𝒦_𝑗 is given by 𝒦_𝑗=𝒦_𝑗−1𝐼−𝐶^𝑇𝐶𝑃_𝑗−1𝐶^𝑇+𝑅⁻¹𝐶𝑃_𝑗−1𝐴^𝑇, (3.44) with the initial condition 𝒦_𝑁−ℎ=𝑃_𝑁−ℎ.

Proof. Using Lemma 3.2, 𝛾_{𝑗+1∣𝑘} in (3.31) for 𝑗≥𝑁−ℎ can be represented as 𝛾_{𝑗+1∣𝑘}=𝛾_𝑗∣𝑘+Ξ_𝑗𝑁_𝑗⁻¹𝐿^𝑇_𝐺,𝑗𝐶^𝑇𝐶𝑃_𝑗𝐶^𝑇+𝑅⁻¹𝑦_{𝑘−𝑁+𝑗}−𝐶𝛽_𝑗∣𝑘. (3.45) If we denote Ξ_𝑗𝑁_𝑗⁻¹𝐿^𝑇_𝐺,𝑗 by 𝒦_𝑗, we have only to prove (3.44). Since Ξ_𝑁−ℎ=𝐿_{𝐺,𝑁−ℎ}, 𝒦_𝑁−ℎ is equal to 𝑃_𝑁−ℎ. Using Lemma 3.2, we can represent 𝒦_𝑗 in a recursive form as follows: 𝒦_𝑗=𝒦_𝑗−1𝐼−𝐶^𝑇𝐶𝑃_𝑗𝐶^𝑇+𝑅⁻¹𝐶𝑃_𝑗−1𝐴^𝑇. (3.46) 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 𝑥_𝑘−ℎ, 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 𝑃₁=0 and 𝛼_1∣𝑘=0, we compute 𝑃_𝑁 and 𝛼_𝑁∣𝑘 recursively in a backward time. From 𝑃_𝑁 and 𝛼_𝑁∣𝑘, we compute 𝑃₀ and 𝛽_0∣𝑘, from which we drive the forward recursive equation for 𝛽_𝑗∣𝑘 to get 𝛾_𝑁∣𝑘=𝑥_{𝑘−ℎ∣𝑘}.

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𝑥_𝑖+1=1.5400+2𝛿_𝑖−0.73790.7379𝛿_𝑖𝑥_𝑖+𝑢0.49210.7594^𝑖+11𝑤_𝑖,𝑦_𝑖=1+𝛿_𝑖1+𝛿_𝑖𝑥_𝑖+𝑣_𝑖, (4.1) where 𝛿_𝑖 is an uncertain model parameter given as𝛿_𝑖=0.1,200≤𝑖≤220,0,otherwise. (4.2) The system noise covariance 𝑄 and the measurement noise covariance 𝑅 are set to 0.01² and 0.027², respectively. The memory size and the fixed-lag size are taken as 𝑁=10 and ℎ=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 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].

Acknowledgments

This work was supported by the Ministry of Knowledge Economy, Rebublic of Korea (ID: 10040132). The material in thid paper was presented in part at the International Symposium on Intelligent Informatics, Qingdao, China, May 2011.