# Enhancement of the Moving Horizon Estimation Performance Based on an Adaptive Estimation Algorithm

Authors: Steve Alan Talla Ouambo; Alexandre Teplaira Boum; Adolphe Moukengue Imano; Jean-Pierre Corriou
Journal: Journal of Control Science and Engineering
 (2021)

---

## Abstract

Although moving horizon estimation (MHE) is a very efficient technique for estimating parameters and states of constrained dynamical systems, however, the approximation of the arrival cost remains a major challenge and therefore a popular research topic. The importance of the arrival cost is such that it allows information from past measurements to be introduced into current estimates. In this paper, using an adaptive estimation algorithm, we approximate and update the parameters of the arrival cost of the moving horizon estimator. The proposed method is based on the least-squares algorithm but includes a variable forgetting factor which is based on the constant information principle and a dead zone which ensures robustness. We show by this method that a fairly good approximation of the arrival cost guarantees the convergence and stability of estimates. Some simulations are made to show and demonstrate the effectiveness of the proposed method and to compare it with the classical MHE.


---

## Body

## 1. Introduction

Since the 1960s, the state estimation of systems by a moving horizon approach has been highlighted. The dazzling success of this estimation technique is based on its ability to explicitly take into account the constraints on the estimate of variables of a dynamic system [1–3]. In [4], the problem of estimating the state of a linear system on a peer-to-peer network of linear sensors is discussed. The proposed approach is fully distributed and scalable and allows for taking into account the constraints on noise and state variables by resorting to the moving horizon estimation paradigm. Each node in the network calculates its local state estimate by minimizing a cost function defined over a sliding window of fixed size. The cost function includes a fused arrival cost which is computed in a distributed way by performing a consensus on the local arrival costs. In [5], the authors have proposed a novel type of iterative partition-based moving horizon estimators, the estimates of which approach those of a centralized moving horizon estimator as the number of iterations increases. It uses a deterministic setting and can handle known inputs as well as bounds on the estimated state. They have derived conditions on the system, its partitions, and the scalar regularization parameter, which guarantee convergence towards the optimal centralized state estimate as well as stability of the estimation error dynamics, even with a finite number of iterations at each time step. The authors, in [6], have investigated on a problem of event-triggered state estimation for networked linear systems. They consider that the stochastic system disturbances and noise are bounded and moving horizon estimation (MHE) is used to handle these constraints and establish an event-based state estimation mechanism that aims to provide good state estimates while reducing the frequencies of both the evaluation of the state estimator and networked communication between the plant and the estimator. In [7], the authors focus on distributed moving horizon estimation (DMHE) for a class of two-time-scale nonlinear systems described in the framework of singularly perturbed systems. By taking advantage of the time-scale separation property, a two-time-scale system is first decomposed into a reduced-order fast system and a reduced-order slow system. The slow system is further decomposed into several interconnected slow subsystems. In the proposed distributed state estimation scheme, a local estimator is designed for each slow subsystem and for the reduced-order fast system.In the moving horizon estimation approach, the state estimate is determined by solving an optimization problem online that minimizes the sum of the squared errors. At a sampling time, when a new measurement is obtainable, the older measuring within the estimation window is discarded and the horizon finite optimization problem is solved again to obtain the new state estimate [8–10]. To take into account the historical data outside the estimation window, a variable that summarizes the information from this data is included in the objective function of the MHE optimization problem called arrival cost. This arrival cost is crucial for the stability and performance of the moving horizon estimation; therefore, its approximation is an active research topic.To reduce the size of the estimation window and the size of the optimization problem, a good approximation of the arrival cost is necessary, ensuring good performance and good robustness. The arrival cost is commonly estimated using a weighted norm of states at the beginning of the estimated horizon window [8, 11, 12]. In [13], the authors proposed to use the Kalman filter approach for updating the term of the arrival cost (in particular the weighting matrix) for the case of linear systems, leading to an inadequate approximation of the weighting matrix leading to relatively poor estimates due to a Gaussian distribution of the Kalman filter. To overcome this problem, the authors in [14] proposed an iterative scheme for updating the arrival cost using the information on the active or inactive constraints from the previous iteration and a quadratic approximation. The idea brings the quadratic approximation as close as possible to the optimal solution of the arrival cost. The fundamental assumption behind the updating scheme is that the constraints should remain unchanged after a certain moment despite the increase of the estimation horizon. For a large estimation window, this assumption is functional giving good results; however, if after smoothing some constraints become inactive, less good results are obtained.In [15], the author uses nonlinear filters based on sampling, including the unscented Kalman filter based on deterministic sampling, the particle filter based on random sampling, and the cell filter based on the aggregate Markov chain to update the arrival cost of the moving horizon estimation and demonstrate the advantages and performance of those filters over the traditional extended Kalman filter approach. However, the computation time is relatively long and better performance is also obtained due to a fairly large horizon length.In more recent works, as in [16], the authors proposed a moving horizon estimation algorithm based on multiple estimation windows. The major advantage of this algorithm is that it reduces the size of the optimization problem while maintaining the performance properties of the estimator and its stability by taking advantage of the inactivity of the constraints. In [17], the authors developed a two-step design strategy, namely, decoupling step and convergence step for networked linear systems with unknown inputs under dynamic quantization effects. In the decoupling step, the decoupling parameter of the moving horizon estimator is designed based on certain assumptions on system parameters and quantization parameters. In the convergence step, by employing a special observability decomposition scheme, the convergence parameters of the moving horizon estimator are achieved such that the estimation error dynamics is ultimately bounded. The conventional moving horizon estimation strategy turns out to be unable to guarantee satisfactory performance as the estimation error depends on external disturbances. In [18], a moving horizon estimation strategy is formulated and the corresponding analytical solution of the state estimation is derived by using the completing the square technique. Based on the lifting approach, the dynamics of the estimation error covariance is analyzed and the desired estimator parameters are obtained by solving a set of linear matrix inequalities.In this work, the weighting matrix is updated using an adaptive estimation algorithm, which is the least-squares algorithm including adead zone and a variable forgetting factor [19, 20]. The dead zone here is used to obtain the robustness; the formula of a forgetting factor based on the principle of constant information can be applied to the process model with unmodeled process uncertainties or external disturbances. Uncertainties can cause the estimates to differ when the standard recursive least-squares algorithm is used, but to prevent this, the inclusion of the dead zone in the estimation algorithm is required. Therefore, the weighting matrix will be calculated in closed loop unlike standard estimation techniques. The effectiveness of the proposed approach is evaluated by carrying out simulation studies on benchmark problems found in the literature.The main contributions of this work aim to (i) show that an adaptive estimation method can be implemented to approximate the arrival cost and have much better efficiency than the traditional MHE, (ii) also show that the proposed scheme guarantees robustness to the estimator for sudden disturbances, and (iii) finally show that the bounds of estimates are crucial when updating the arrival cost. The paper is organized as follows. In Section2, we present the moving horizon estimation problem, followed by the arrival cost update mechanism in Section 3. Section 4 presents a stability analysis; the results are discussed in Section 5, and finally, Section 6 will be devoted to the conclusions.


## 2. Preliminaries and Problem Statement

Notation 1.
In the following, we will consider the following notations,ξ∈ℝn is a column vector and ξT∈ℝn will be its transpose. ξj|ℓ denotes a finite sequence of vectors over a given index up to ℓ. Let B∈ℝn×n be a n×m matrix and BT is its transpose. If B∈ℝn×n, then B−1 denotes its inverse. Let P∈ℝn×n be a symmetric real matrix, P is said to be positive definite if ΞTPΞ>0 and ∀Ξ∈ℝn. The Euclidean matrix or vector norm is denoted by ⋅, for instance, ξ∈ℝn and P∈ℝn×n are positive definite; it follows that ξp=ξTPξ.Consider the linear time invariant discrete system:(1)xℓ+1=Axℓ+Gwℓ,yℓ=Cxℓ+vℓ,where xℓ∈Xθ⊆ℝn is the state vector, wℓ∈Wθ⊆ℝp is the process noise vector, yℓ∈Yθ⊆ℝm is the measurement vector, and vℓ∈Vθ⊆ℝm is the measurement noise vector. The process noise wℓ and the measurement noise vℓ are assumed to be bounded and unknown eventually, i.e., w∈WθϑWθmax and v∈VθϑVθmax for ϑWθmax,ϑVθmax∈ℝ+2, with Xθ, Yθ, Wθ, and Vθ being compact and convex sets considered closed with 0∈Xθ, 0∈Yθ, 0∈Wθ, and 0∈Vθ.Using all measures for the full information estimator to estimate the system statesx^ℓ increases the size of the optimization problem over time and leads to heavy calculation.To overcome this problem, the moving horizon estimation was developed to take into account a certain amount of data and possibly to make dynamic updates of certain parameters in an estimation window moving in time.A penalty term called arrival cost is essential in the formulation of the MHE because it takes into account past measurements to estimate current states. The problem of the MHE aims to determine, at a timeℓ, an estimate x^ℓ|ℓ of the current state and process noises which solve the optimization problem.Forℓ≤N, the full information problem arises:(2)minw^ℓ,x^0|ℓΨ=x^0|ℓ−x˜0P0−12+∑j=0ℓ−1w^j|ℓQ−12+∑j=0ℓv^j|ℓR−12,s.t.x^j+1|ℓ=Ax^j|ℓ+w^j|ℓ,yj=Cx^j|ℓ+v^j|ℓ,x^j|ℓ∈Xθ,w^j|ℓ∈Wθ,v^j|ℓ∈Vθ.However, for allℓ>N, the optimization problem becomes(3)minw^ℓ,x^ℓ−N|ℓΨ=Γℓ−N|ℓx^ℓ−N|ℓ+∑j=ℓ−Nℓ−1w^j|ℓQ−12+∑j=ℓ−N+1ℓv^j|ℓR−12,.s.t.x^j+1|ℓ=Ax^j|ℓ+w^j|ℓ,yj=Cx^j|ℓ+v^j|ℓ,x^j|ℓ∈Xθ,w^j|ℓ∈Wθ,v^j|ℓ∈Vθ,where the optimal estimated is x^ℓ−j|ℓ, w^j|ℓ is the process noise estimate based on available measurements at time ℓyℓ−j, and v^j|ℓ=yj−Cx^j|ℓ are the estimation residuals. To ensure a fairly good robust stability of the MHE, a judicious and adequate choice of Γℓ−N|ℓ⋅ and its parameters is necessary [21].Usually, a relatively long estimation window is used to overcome a poor estimate caused by an inadequate cost function approximation [22]; this case is strongly encountered in the conventional MHE. It can also be noted that, for a short estimation window, an approximation of the optimal cost function is required to properly incorporate the past information into the problem, in case of which a poor estimate will be revealed.In many dynamical systems, in particular linear systems, a Gaussian distribution is assumed and the arrival cost function is updated using a Kalman filter [23]. Besides, the MHE is equivalent to the Kalman filter for unconstrained linear systems with indeed the same Gaussian distribution. In the following, we will present the proposed mechanism for updating the arrival cost parameters of the MHE.


## 3. Arrival Cost Update Algorithm

For the constrained MHE problem, it is difficult to find an exact analytical expression of the arrival cost; therefore, an approximation of the arrival cost used for the unconstrained problem will be used; this quadratic approximation is of the following form:(4)Γℓ−N=x^ℓ−N−x˜ℓ−NPℓ−N−12,where Pℓ−N is the weighting matrix and x˜ℓ−N is the initial state. Both define the approximation of the arrival cost function. The proposed approach for updating the arrival cost in this work, based on an adaptive algorithm, is formulated as an update of the following elements: x˜ℓ−N and Pℓ−N.The update of the initial statex˜ℓ−N will be made by a smooth update, which implies that once the estimate leaves the estimation horizon, the state of the constraints does not change [14, 24]:(5)x˜ℓ−N=x^ℓ−N|ℓ−1∗.Regarding the update of the weighting matrixPℓ−N, it is necessary to notice from the outset that, in a stochastic interpretation of the MHE problem, the weighting matrix can be seen as the covariance of the state xℓ−N [8]. For unconstrained linear systems, analytical solutions exist to compute the weighting matrix; on the contrary, for constrained linear systems, a deduction of the approximate solutions is made, which are fundamentally and based on a recursive process of updating useful information from the system at a given time. In the proposed adaptive estimation method, an update based on a recursive process relying on states (or not) and measurements is used [19, 25]:(6)Pℓ−N−1=λℓPℓ−N−1−1+αℓλℓx^ℓ−N|ℓ−1x^ℓ−N|ℓ−1T.WithP0>0, where λℓ is the forgetting factor with 0<λℓ≤1, αℓ is a function which varies between 0 and 1 called the dead zone. Equation (6) produces a recursively real-time estimate of Pℓ−N, by updating Pℓ−N−1 (previous estimate) with a dynamic temporal average over time of matrix x^ℓ−N|ℓ−1x^ℓ−N|ℓ−1T. We can qualify this updating mechanism as a dynamic temporal filter overtime with P0, the initial condition and x^ℓ−N|ℓ−1x^ℓ−N|ℓ−1T are the input.Lemma 1.
If a matrixZ>0, then there exists an invertible matrix A such that ATZA=I, which is the identity matrix. Moreover, AAT=Z−1 which exists.Proof.
See Lemma 7 in [26].
Note that the forgetting factorλℓ and the function αℓ each have a different role. When λℓ≈1, the algorithm tends to keep the old data in Pℓ−N−1; in this case, all information is included, improving the estimation problem. On the contrary, when λℓ≪1, the algorithm tends to discard the old information in Pℓ−N−1, considering recent data to estimate Pℓ−N [27]. This case allows avoiding many old data that could negatively affect the estimator performance. The dead zone function αℓ is used to ensure robustness when estimating.
By the matrix inversion Lemma, we can rewrite (6) as follows:(7)Pℓ−N=1λℓPℓ−N−1−αℓPℓ−N−1x^ℓ−N|ℓ−1x^ℓ−N|ℓ−1TPℓ−N−1λℓ+x^ℓ−N|ℓ−1TPℓ−N−1x^ℓ−N|ℓ−1.
A feedback loop is introduced into the MHE problem by this way to computePℓ−N; this way makes it possible to ensure the convergence of the estimates x^ℓ−N|ℓ and w^jj=ℓ−N,…,N and to adjust Pℓ−N−1 according to the amount of available information. The adaptive algorithm used proposes a variable forgetting factor varying between 0 and 1 and under certain conditions to be respected and a dead zone, both allowing Pℓ−N to be updated as follows [19]:(8)λℓ=N0+Δℓ2N0+M1ℓ,ifΔℓ2M1ℓ<1,1,otherwise,(9)αℓ=1,ifΔℓ21+μℓλℓ≤eℓ2,0,otherwise,where N0 is a positive design parameter chosen which is quite large: M1ℓ=eℓ2/1+μℓ, μℓ=x^ℓ−N|ℓ−1TPℓ−N−1x^ℓ−N|ℓ−1, and eℓ=yℓ−N−y^ℓ−N, where Δℓ=x^ℓ−N|ℓ−1d1+d2/nℓ−1, with nℓ−1=max1,x^ℓ−N|ℓ−1 and d1 and d2 are also design parameters, which are chosen positive.
Any sequenceλℓ is bounded by some σ∗>0. Thus, the size of dead zone Δℓ21+μℓ/λℓ is also bounded because μℓ is bounded too.
The forgetting factor proposed in [19] is based on the constant information principle, which is adequate for external disturbances and unmodeled system uncertainties.
The following algorithm summarizes the proposed moving horizon estimation estimator based on the adaptive estimation method for the arrival cost update.Algorithm 1:x˜0, P0,Q,R>0.
d1,d2,N0>0 and N≥2Initializationx˜ℓ−N←x˜0Pℓ−N←P0forℓ=0,1,2…,doGet measurementsyℓifℓ<Ntheny←y0,y1,…,yℓTSolveΨℓ with x˜0, P0 and yElsey←yℓ−N,…,yℓTSolveΨℓN with x˜ℓ−N, Pℓ−N and yx˜ℓ−N←x^ℓ−N+1|ℓUpdatingPℓ−N using (7)–(9)end ifGet current estimatesx^j|ℓj=ℓ−N,…,ℓend for


## 4. Stability Analysis

The stability of the estimator requires that the prediction error converges to zero for the following nominal system:(10)xℓ+1=Axℓ,yℓ=Cxℓ,with no measurement and state noise. In [8], the authors require that the system has to evolve according to the constraints because, for the constrained systems’ case, a poor range of choice of constraints can prevent convergence to the real values of the system.Lemma 2.
The estimation algorithm assumes thatsupδℓ≤d1 and supwℓ≤d2, with δℓ as the system uncertainty vector; then, the following properties are valid:(i)
Pℓ converges and is uniformly bounded(ii)
supeℓ2−1+KΔℓ2≤0, ∀K≥0, with eℓ=yℓ−y^ℓProof.
By recursive calculation and the matrix inversion lemma, we get from (7)(11)Pℓ−1=σ0,ℓP0−1+∑i=1ℓαiσ0,ix^ℓ|i−1x^ℓ|i−1T.
It follows that(12)Pℓ≤1+v0N0P0,which implies the uniform boundedness of Pℓ; by the property 0<σ∗≤σ0,ℓ, it is easily proven from (11) that Pℓ uniformly converges.Remark 1.
The matrixPℓ owns a steady-state solution, i.e.,(13)limℓ⟶∞Pℓ=P∞.Proof.
The update expression ofPℓ−N given by equation (7) can be seen as a particular case of Ricatti equation:(14)λℓPℓ−N=Pℓ−N−1−αℓPℓ−N−1x^ℓ−N|ℓ−1x^ℓ−N|ℓ−1TPℓ−N−1λℓ+x^ℓ−N|ℓ−1TPℓ−N−1x^ℓ−N|ℓ−1=Pℓ−N−1−αℓPℓ−N−1x^ℓ−N|ℓ−1λℓ+x^ℓ−N|ℓ−1TPℓ−N−1x^ℓ−N|ℓ−1−1x^ℓ−N|ℓ−1TPℓ−N−1.
It is stable if and only if the eigenvalues of the transfer matrix are strictly inside the unit circle of the complex plane:(15a)Λ=I−αℓPℓ−N−1x^ℓ−N|ℓ−1x^ℓ−N|ℓ−1Tλℓ+x^ℓ−N|ℓ−1TPℓ−N−1x^ℓ−N|ℓ−1,(15b)Λ=I−Λ˜,since(16)Pℓ−N−1x^ℓ−N|ℓ−1x^ℓ−N|ℓ−1T≤Pℓ−N−1x^ℓ−N|ℓ−1x^ℓ−N|ℓ−1T,x^ℓ−N|ℓ−1TPℓ−N−1x^ℓ−N|ℓ−1≤x^ℓ−N|ℓ−1TPℓ−N−1x^ℓ−N|ℓ−1.
Hence,(17)Λ˜≤1.
All eigenvalues of matrixΛ˜ are greater than 0 and bounded by 1; then, the eigenvalues of (4) are also greater than 0 and bounded by 1. Hence, the matrix Pℓ, therefore, has a steady-state solution P∞.Theorem 1.
Letvℓ=x˜ℓ∗TPℓ−1x˜ℓ∗ be a Lyapunov function for the moving horizon estimator with the adaptive arrival cost, where x˜ℓ∗=xℓ−x^ℓ and the matrix Pℓ−1 is given by equation (6), i.e., Pℓ−1=λℓPℓ−1−1+αℓ/λℓx^ℓ−N|ℓ−1x^ℓ−N|ℓ−1T.Proof.
λℓvℓ−1≤vℓ−1; we have(18)vℓ−λℓvℓ−1=x˜ℓ∗TPℓ−1x˜ℓ∗−λℓx˜ℓ−1∗TPℓ−1x˜ℓ−1∗=Ax˜ℓ−1∗TPℓ−1Ax˜ℓ−1∗−λℓx˜ℓ−1∗TPℓ−1−1x˜ℓ−1∗=x˜ℓ−1∗TATPℓ−1Ax˜ℓ−1∗−λℓx˜ℓ−1∗TPℓ−1−1x˜ℓ−1∗=x˜ℓ−1∗TATλℓPℓ−1−1+αℓλℓx^ℓ−N|ℓ−1x^ℓ−N|ℓ−1TAx˜ℓ−1∗−λℓx˜ℓ−1∗TPℓ−1−1x˜ℓ−1∗vℓ−λℓvℓ−1=λℓx˜ℓ−1∗TATPℓ−1−1A+αℓλℓATx^ℓ−N|ℓ−1x^ℓ−N|ℓ−1TA−Pℓ−1−1x˜ℓ−1∗,Pℓ,and ATx^ℓ−N|ℓ−1x^ℓ−N|ℓ−1TA are symmetric positive definite and symmetric positive semidefinite, respectively, for all ℓ; then, if (10) is observable, it follows that the MHE with the proposed adaptive arrival cost update is an asymptotically stable observer.
Regarding bounded stability, we can note that the arrival cost satisfies(19)λminPℓ−N−1x^ℓ−N|ℓ−x˜ℓ−N≤x^ℓ−N|ℓ−x˜ℓ−NPℓ−N−1≤λmaxPℓ−N−1x^ℓ−N|ℓ−x˜ℓ−NλminQ−1w^≤w^Q−1≤λmaxQ−1w^λminR−1v^≤v^R−1≤λmaxR−1v^.
The estimator with adaptive arrival cost update is robust globally asymptotically stable under bounded disturbances; besides, iflimℓ⟶∞vℓ=0 and limℓ⟶∞wℓ=0, then limℓ⟶∞x^ℓ−x˜ℓ=0.Proof.
cf. Theorem 7 in [28].


## 5. Simulation Results

To demonstrate the performance of the proposed method, we adopted discrete-time models used by [13, 29].


### 5.1. Example 1

(20)xℓ+1=0.990.2−0.10.3xℓ+01wℓ,yℓ=1−3xℓ+vℓ,where wℓ=sℓ with sℓ a sequence of independent, normally distributed random values with zero mean and covariance equal to 1, and vℓ is also a sequence of independent, normally distributed random values with zero mean and covariance equal to 0.1. The initial state x˜0 with zero mean is also normally distributed with a unit covariance. The constrained estimation problem here is formulated with Q=I, R=0.01, P0=0.5I, and x˜0=0.5,0.5T. We choose N=10 and d1=0.01 and d2=0.005 and N0=50 as design parameters.The proposed algorithm is compared with the classical MHE (MHEKF) which propagates the initial state and the arrival cost estimate using Kalman filter. The sum square estimation error is used as a benchmark:(21)ζi=∑l=0ℓxli−x^li2,where xi denotes the ith component of the vector x. Based on 70 trials, the average sum square estimation error was computed. The results are shown in Table 1, where it can be seen that the performance of MHEAD is superior to the classical MHE.Table 1
Average sum square estimation error (example 1).
MHEADMHEKFx144.0183.16x270.97125.10Figures1(a) and 1(b) show the states and their estimates x1 and x2 with the horizon length N=10.Figure 1
Comparison of estimators for model (20) with horizon length N=10. (a) Estimation of x1. (b) Estimation of x2. (c) x1−x^12. (d)x2−x^22.
(a)(b)(c)(d)We can see that despite the relatively short horizon length chosen, theMHEAD estimate is closer than MHEKF. The estimates of MHEAD tend to follow MHEKF for x1. On the contrary, for the state x2, we observe a divergent behavior. At every sample time, the square error is shown in Figures 1(c) and 1(d); we can observe that the error of the proposed method is lower than MHEKF, and we can also notice that the performance of MHEKF is quite degraded, probably, because of the relatively short horizon length.In Figures2(a) and 2(b), we can observe the evolution of the forgetting factor and the dead zone function, respectively.Figure 2
Evolution of the sequence of (a) forgetting factorλℓ and (b) dead zone function αℓ.
(a)(b)The forgetting factor is varying between 0 and 1, providing a good adaptation capability and allowing the proposed algorithm to insert pertinent information into the data and follow sudden variations. The dead zone function makes some variations between0 and 1 to allow a rapid adaptation of the weighting matrix when a sudden change occurs. We can note that when the estimation error is large, the dead zone function is equal to 1, and when the error is relatively low, the dead zone function is equal to 0.Figures3(a) and 3(b) allow appreciating the evolution of the trace of the weighting matrix Pℓ for both estimators.Figure 3
Evolution of tracePℓ. (a) Trace Pℓ−MHEKF. (b) Trace Pℓ−MHEAD.
(a)(b)TheMHEKF weighting matrix converges to the steady-state value and does not change anymore, while MHEAD adapts its weighting matrix according to variations of the process noise.A particular case is studied and presented in Figure4, where the process noise covariance σw varies with time. The covariance will take values 0.9 when 0≤ℓ<5, 0.55 when 5≤ℓ<20, 0.8 when 20≤ℓ<40, 1.3 when 40≤ℓ<60, and 0.4 when 60≤ℓ<80. We can see in Figures 4(a) and 4(b) that the proposed approach still performs efficiently better than the conventional MHE despite the variation with time of the covariance. We can also note that compared to the previous scenario, the estimate is less good for both estimators, but this can be explained by the fact that the estimators have no knowledge of this happening. Figures 4(c) and 4(d) show the trace of the matrix Pℓ used by MHEKF and MHEAD, respectively. We can observe that the MHEKF weighting matrix converges to a steady value at 20s and does not change anymore, while the proposed method tries to adapt its weighting matrix according to the variations of the process noise.Figure 4
Comparison of estimators for model (16) with horizon length N=10 and variable σw. (a) Estimation of x1. (b) Estimation of x2. (c) Trace Pℓ−MHEKF. (d) Trace Pℓ−MHEAD. (e) Forgetting factor. (f) Dead zone function.
(a)(b)(c)(d)(e)(f)It can be seen that the trace ofPℓ tends to stabilize when the error is small. On the contrary, the trace of the weighting matrix becomes smaller when the difference between the estimate and the real state grows, ensuring good adaptation despite the influence of varying disturbances.In Figures4(e) and 4(f), we can see the evolution of the forgetting factor and the dead zone, respectively.We can see that the evolution of the forgetting factor is independent of the dead zone, but it evolves according to the variation of the process noise. We also notice that the forgetting factor tries to adapt to the variations while keeping its value close to1. On the contrary, the behavior of the dead zone is such that when the estimation error is small, the dead zone function tends to 0 while when the error grows, the dead zone is equal to 1, so to allow a good adaptation and provide an acceptable estimation.Figure5 shows the performance of MHEAD with different values of the horizon length. We can notice that, for small values of N, the estimation error is not negligible, while for large values of N, the estimation error is very appreciable. Table 2 allows us to assess the results quantitatively.Figure 5
Performance ofMHEAD with N=3, N=6, and N=10 and variable σw. (a) Estimation of x1. (b) Estimation of x2.
(a)(b)Table 2
Comparison of the effect of the horizon size on the average sum square estimation error for statesx1 and x2.
MHEADx1x2N=377.5288.41N=667.0380.05N=1048.1171.15


### 5.2. Example 2

Still in the same logic of showing the performance of the proposed approach, the example of the following system is provided:(22)xℓ+1=Axℓ+Buℓ+Bddℓ+wℓyℓ=Cxℓ+vℓ,where(23)A=0.7081−128.29481.28290.00080.8513−0.01850.00378.45970.9154,B=1110.27840.53092.4259,Bd=15.5509−0.224423.2170,C=0.0085−0.64150.00640.00000.04230.0096.Here, disturbances are modeled with the same distribution as example 1. The constrained estimation problem is formulated withQ=I, R=20002, P0=10I, and x¯0=1.5−0.50T.uℓ=4Hℓ−0.15−4Hℓ−2.5 and dℓ=−37Hℓ−0.25+37Hℓ−0.7+37Hℓ−1.3−37Hℓ−1.7−40Hℓ−4.5 are the control input and the disturbance, respectively. H represents the Heaviside function. We also choose the horizon N=10 in this example and d1=0.00005, d2=0.001, and N0=5000 as design parameters. Figures 6(a)–6(c) show the states and their estimates x1, x2, and x3. Figures 6(d)–6(f) show the square error of states and their estimates. From these figures, we can see that the proposed method (MHEAD) gives a good estimate compared to the conventional MHE (MHEKF); this is because, in this example, the inputs of the system are bounded and persistent. We can simply say the MHEAD estimates can track properly the states. The divergence of the conventional MHE is due to the inclusion of a large value of the disturbance (sudden variation) at t=4.5s; on the contrary, the proposed algorithm adapts to slow or/and sudden changes. This can easily be explained by the fact that the crucial parameters of the proposed algorithm can follow slow or/and sudden changes or disturbances and provide good adaptability. Figures 7(a) and 7(b) show the evolution of the forgetting factor and the dead zone function, respectively. The same remarks as in the previous example are also considered in this case. It can be seen that the best performance of the proposed method occurred when the forgetting factor was equal to 1, and the dead zone function equal to 0, when the disturbances were less important. When the disturbances were considerable, the forgetting factor tried to adjust to them. Figures 8(a) and 8(b) allow appreciating the evolution of the trace of the weighting matrix Pℓ for both estimators. We can see that the evolution of the trace of the classical MHE weighting matrix has difficulty reaching the steady state unlike the previous example, which justifies the weak reconstruction; on the contrary, the evolution of the weighting matrix of MHEAD is smoother; that is, when the disturbances are less consequent, the evolution of the trace of the weighting matrix of MHEAD behaves as a classical constrained optimization problem. When the noise increases, the proposed estimator does not follow the system and the forgetting factor decreases and moves away from 1 to allow a fairly fast adaptation of Pℓ. Based on 70 trials, the average sum square estimation error was computed, and the results are shown in Table 3. Likewise, in Figure 9, we show the performance of MHEAD with different values of the horizon length. It can be seen that MHEAD errors do not change significantly with the different horizon sizes, showing that this algorithm seems to be a good method for the arrival cost approximation, even for short horizons. Table 4 allows us to assess the results quantitatively.Figure 6
Comparison of estimators for model (22) with horizon length N=10. (a) Estimation of x2. (b) Estimation of x2. (c) Estimation of x3. (d)x1−x^12. (e)x2−x^22. (f)x3−x^32.
(a)(b)(c)(d)(e)(f)Figure 7
Evolution of the sequence of (a) forgetting factorλℓ and (b) dead zone function αℓ.
(a)(b)Figure 8
Evolution of tracePℓ. (a) Trace Pℓ−MHEKF. (b) Trace Pℓ−MHEAD.
(a)(b)Table 3
Average sum square estimation error (example 2).
MHEADMHEKFx175.37185.10x2115.02201.50x3129.25198.22Figure 9
Performance ofMHEAD with N=3, N=6, and N=10: (a) estimation of x1; (b) estimation of x2; (c) estimation of x3.
(a)(b)(c)Table 4
Comparison of the effect of the horizon size on the average sum square estimation error for states:x1, x2, and x3.
MHEADx1x2x3N=378.01122.02133.41N=676.95119.55132.01N=1075.37115.02129.25


## 5.1. Example 1

(20)xℓ+1=0.990.2−0.10.3xℓ+01wℓ,yℓ=1−3xℓ+vℓ,where wℓ=sℓ with sℓ a sequence of independent, normally distributed random values with zero mean and covariance equal to 1, and vℓ is also a sequence of independent, normally distributed random values with zero mean and covariance equal to 0.1. The initial state x˜0 with zero mean is also normally distributed with a unit covariance. The constrained estimation problem here is formulated with Q=I, R=0.01, P0=0.5I, and x˜0=0.5,0.5T. We choose N=10 and d1=0.01 and d2=0.005 and N0=50 as design parameters.The proposed algorithm is compared with the classical MHE (MHEKF) which propagates the initial state and the arrival cost estimate using Kalman filter. The sum square estimation error is used as a benchmark:(21)ζi=∑l=0ℓxli−x^li2,where xi denotes the ith component of the vector x. Based on 70 trials, the average sum square estimation error was computed. The results are shown in Table 1, where it can be seen that the performance of MHEAD is superior to the classical MHE.Table 1
Average sum square estimation error (example 1).
MHEADMHEKFx144.0183.16x270.97125.10Figures1(a) and 1(b) show the states and their estimates x1 and x2 with the horizon length N=10.Figure 1
Comparison of estimators for model (20) with horizon length N=10. (a) Estimation of x1. (b) Estimation of x2. (c) x1−x^12. (d)x2−x^22.
(a)(b)(c)(d)We can see that despite the relatively short horizon length chosen, theMHEAD estimate is closer than MHEKF. The estimates of MHEAD tend to follow MHEKF for x1. On the contrary, for the state x2, we observe a divergent behavior. At every sample time, the square error is shown in Figures 1(c) and 1(d); we can observe that the error of the proposed method is lower than MHEKF, and we can also notice that the performance of MHEKF is quite degraded, probably, because of the relatively short horizon length.In Figures2(a) and 2(b), we can observe the evolution of the forgetting factor and the dead zone function, respectively.Figure 2
Evolution of the sequence of (a) forgetting factorλℓ and (b) dead zone function αℓ.
(a)(b)The forgetting factor is varying between 0 and 1, providing a good adaptation capability and allowing the proposed algorithm to insert pertinent information into the data and follow sudden variations. The dead zone function makes some variations between0 and 1 to allow a rapid adaptation of the weighting matrix when a sudden change occurs. We can note that when the estimation error is large, the dead zone function is equal to 1, and when the error is relatively low, the dead zone function is equal to 0.Figures3(a) and 3(b) allow appreciating the evolution of the trace of the weighting matrix Pℓ for both estimators.Figure 3
Evolution of tracePℓ. (a) Trace Pℓ−MHEKF. (b) Trace Pℓ−MHEAD.
(a)(b)TheMHEKF weighting matrix converges to the steady-state value and does not change anymore, while MHEAD adapts its weighting matrix according to variations of the process noise.A particular case is studied and presented in Figure4, where the process noise covariance σw varies with time. The covariance will take values 0.9 when 0≤ℓ<5, 0.55 when 5≤ℓ<20, 0.8 when 20≤ℓ<40, 1.3 when 40≤ℓ<60, and 0.4 when 60≤ℓ<80. We can see in Figures 4(a) and 4(b) that the proposed approach still performs efficiently better than the conventional MHE despite the variation with time of the covariance. We can also note that compared to the previous scenario, the estimate is less good for both estimators, but this can be explained by the fact that the estimators have no knowledge of this happening. Figures 4(c) and 4(d) show the trace of the matrix Pℓ used by MHEKF and MHEAD, respectively. We can observe that the MHEKF weighting matrix converges to a steady value at 20s and does not change anymore, while the proposed method tries to adapt its weighting matrix according to the variations of the process noise.Figure 4
Comparison of estimators for model (16) with horizon length N=10 and variable σw. (a) Estimation of x1. (b) Estimation of x2. (c) Trace Pℓ−MHEKF. (d) Trace Pℓ−MHEAD. (e) Forgetting factor. (f) Dead zone function.
(a)(b)(c)(d)(e)(f)It can be seen that the trace ofPℓ tends to stabilize when the error is small. On the contrary, the trace of the weighting matrix becomes smaller when the difference between the estimate and the real state grows, ensuring good adaptation despite the influence of varying disturbances.In Figures4(e) and 4(f), we can see the evolution of the forgetting factor and the dead zone, respectively.We can see that the evolution of the forgetting factor is independent of the dead zone, but it evolves according to the variation of the process noise. We also notice that the forgetting factor tries to adapt to the variations while keeping its value close to1. On the contrary, the behavior of the dead zone is such that when the estimation error is small, the dead zone function tends to 0 while when the error grows, the dead zone is equal to 1, so to allow a good adaptation and provide an acceptable estimation.Figure5 shows the performance of MHEAD with different values of the horizon length. We can notice that, for small values of N, the estimation error is not negligible, while for large values of N, the estimation error is very appreciable. Table 2 allows us to assess the results quantitatively.Figure 5
Performance ofMHEAD with N=3, N=6, and N=10 and variable σw. (a) Estimation of x1. (b) Estimation of x2.
(a)(b)Table 2
Comparison of the effect of the horizon size on the average sum square estimation error for statesx1 and x2.
MHEADx1x2N=377.5288.41N=667.0380.05N=1048.1171.15


## 5.2. Example 2

Still in the same logic of showing the performance of the proposed approach, the example of the following system is provided:(22)xℓ+1=Axℓ+Buℓ+Bddℓ+wℓyℓ=Cxℓ+vℓ,where(23)A=0.7081−128.29481.28290.00080.8513−0.01850.00378.45970.9154,B=1110.27840.53092.4259,Bd=15.5509−0.224423.2170,C=0.0085−0.64150.00640.00000.04230.0096.Here, disturbances are modeled with the same distribution as example 1. The constrained estimation problem is formulated withQ=I, R=20002, P0=10I, and x¯0=1.5−0.50T.uℓ=4Hℓ−0.15−4Hℓ−2.5 and dℓ=−37Hℓ−0.25+37Hℓ−0.7+37Hℓ−1.3−37Hℓ−1.7−40Hℓ−4.5 are the control input and the disturbance, respectively. H represents the Heaviside function. We also choose the horizon N=10 in this example and d1=0.00005, d2=0.001, and N0=5000 as design parameters. Figures 6(a)–6(c) show the states and their estimates x1, x2, and x3. Figures 6(d)–6(f) show the square error of states and their estimates. From these figures, we can see that the proposed method (MHEAD) gives a good estimate compared to the conventional MHE (MHEKF); this is because, in this example, the inputs of the system are bounded and persistent. We can simply say the MHEAD estimates can track properly the states. The divergence of the conventional MHE is due to the inclusion of a large value of the disturbance (sudden variation) at t=4.5s; on the contrary, the proposed algorithm adapts to slow or/and sudden changes. This can easily be explained by the fact that the crucial parameters of the proposed algorithm can follow slow or/and sudden changes or disturbances and provide good adaptability. Figures 7(a) and 7(b) show the evolution of the forgetting factor and the dead zone function, respectively. The same remarks as in the previous example are also considered in this case. It can be seen that the best performance of the proposed method occurred when the forgetting factor was equal to 1, and the dead zone function equal to 0, when the disturbances were less important. When the disturbances were considerable, the forgetting factor tried to adjust to them. Figures 8(a) and 8(b) allow appreciating the evolution of the trace of the weighting matrix Pℓ for both estimators. We can see that the evolution of the trace of the classical MHE weighting matrix has difficulty reaching the steady state unlike the previous example, which justifies the weak reconstruction; on the contrary, the evolution of the weighting matrix of MHEAD is smoother; that is, when the disturbances are less consequent, the evolution of the trace of the weighting matrix of MHEAD behaves as a classical constrained optimization problem. When the noise increases, the proposed estimator does not follow the system and the forgetting factor decreases and moves away from 1 to allow a fairly fast adaptation of Pℓ. Based on 70 trials, the average sum square estimation error was computed, and the results are shown in Table 3. Likewise, in Figure 9, we show the performance of MHEAD with different values of the horizon length. It can be seen that MHEAD errors do not change significantly with the different horizon sizes, showing that this algorithm seems to be a good method for the arrival cost approximation, even for short horizons. Table 4 allows us to assess the results quantitatively.Figure 6
Comparison of estimators for model (22) with horizon length N=10. (a) Estimation of x2. (b) Estimation of x2. (c) Estimation of x3. (d)x1−x^12. (e)x2−x^22. (f)x3−x^32.
(a)(b)(c)(d)(e)(f)Figure 7
Evolution of the sequence of (a) forgetting factorλℓ and (b) dead zone function αℓ.
(a)(b)Figure 8
Evolution of tracePℓ. (a) Trace Pℓ−MHEKF. (b) Trace Pℓ−MHEAD.
(a)(b)Table 3
Average sum square estimation error (example 2).
MHEADMHEKFx175.37185.10x2115.02201.50x3129.25198.22Figure 9
Performance ofMHEAD with N=3, N=6, and N=10: (a) estimation of x1; (b) estimation of x2; (c) estimation of x3.
(a)(b)(c)Table 4
Comparison of the effect of the horizon size on the average sum square estimation error for states:x1, x2, and x3.
MHEADx1x2x3N=378.01122.02133.41N=676.95119.55132.01N=1075.37115.02129.25


## 6. Conclusion

In this paper, an alternative method to approximate the arrival cost for the moving horizon estimation has been studied on linear systems. The proposed approach is based on the least-squares algorithm but includes a variable forgetting factor which is based on the constant information principle and a dead zone which ensures robustness. The use of CasADi toolbox [30] was used for the implementation of the optimization algorithm for real-time operation. The results obtained show that the adaptive estimation algorithm allows a good approximation of the arrival cost for the moving horizon estimation. This can easily be explained by the fact that the crucial parameters of the proposed method can follow slow or/and sudden changes or disturbances and provide good adaptability. Since the update mechanism does not depend on the model, the proposed method can be extended to nonlinear systems.



---

*Source: 3776506-2021-12-29.xml*