Robust stability of moving horizon estimation for nonlinear systems with bounded disturbances using adaptive arrival cost

TL;DR

This paper proves the robust stability and convergence of a moving horizon estimator with adaptive arrival cost for nonlinear systems under bounded disturbances, supported by simulations comparing it to existing methods.

Contribution

It introduces a novel adaptive arrival cost approach that ensures robust stability and convergence in nonlinear moving horizon estimation.

Findings

Proves global asymptotic stability with bounded disturbances

Demonstrates convergence to true state with vanishing disturbances

Shows improved estimator performance through simulations

Abstract

In this paper, the robust stability and convergence to the true state of moving horizon estimator based on an adaptive arrival cost are established for nonlinear detectable systems. Robust global asymptotic stability is shown for the case of non-vanishing bounded disturbances whereas the convergence to the true state is proved for the case of vanishing disturbances. Several simulations were made in order to show the estimator behaviour under different operational conditions and to compare it with the state of the art estimation methods.

Tables

Table 1. Table 1: Example 1 averaged MSE over 300 trials.

	FIE MAX	ADAP	MAX	EKF
$x_0$	0.02040	0.02176	0.02206	0.02296
$x_1$	0.00135	0.00151	0.00156	0.00154

Table 2. Table 2: Example 1 aver. MSE over 300 trials with variable measurement noise.

	ADAP	MAX	FIE MAX
$x_0$	0.02068	0.03067	0.00761
$x_1$	0.00290	0.00335	0.00068

Table 3. Table 3: Example 2 averaged MSE over 300 trials and different horizon size.

	N=2		N=5		N=10
	ADAP	MAX	ADAP	MAX	ADAP	MAX	FIE
$x_0$	0.18808	0.58652	0.03367	0.04615	0.00171	0.00772	0.00024
$x_1$	0.23037	0.66768	0.04074	0.05077	0.00285	0.00951	0.00120

Equations

$$\begin{array}{rl}\gamma \left(a_1+a_2+\dots +a_n\right)& \le \gamma \left(n{a}_1\right)+\dots +\gamma \left(n{a}_n\right),\\ \beta \left(a_1+a_2+\dots +a_n,k\right)& \le \beta \left(n{a}_1,k\right)+\dots +\beta \left(n{a}_n,k\right).\end{array}$$

$$\begin{array}[]{rcl}\mathcal{C}_{w}\coloneqq\{\boldsymbol{w}\in\mathcal{W}\left(w_{\max}\right)|\,\boldsymbol{w}\text{ is convergent}\}\end{array}$$

$$\begin{array}{rl}{x}_{k+1}& =f\left(x_k,w_k\right),x_0={\mathtt{x}}_0\\ y_k& =h\left(x_k\right)+v_k,\end{array}$$

$$\begin{array}[]{c}\underset{\hat{x}_{k-N|k},\boldsymbol{\hat{w}}_{j|k}}{\operatorname{min}}\,\Psi\coloneqq\Gamma_{k-N|k}\left(\hat{x}_{k-N|k}\right)+\sum\limits_{j=k-N}^{k-1}\ell\left(\hat{w}_{j|k},\hat{v}_{j|k}\right)\\ \text{s.t.}\left\{\begin{array}[]{l}\begin{array}[]{rll}\hat{x}_{j+1|k}=&f\left(\hat{x}_{j|k},\hat{w}_{j|k}\right),&j\in\mathbb{Z}_{\left[k-N,k-1\right]}\\ y_{j}=&h\left(\hat{x}_{j|k}\right)+\hat{v}_{j|k},&\\ \end{array}\\ \hat{x}_{j|k}\in\mathcal{X},\,\hat{w}_{j|k}\in\mathcal{W},\ \hat{v}_{j|k}\in\mathcal{V},\end{array}\right.\end{array}$$

$$\begin{array}[]{c}\underset{\hat{x}_{k-N|k},\boldsymbol{\hat{w}}}{\operatorname{min}}\,\Psi\coloneqq\Gamma_{0|k}\left(\hat{x}_{0|k}\right)+\sum\limits_{j=0}^{k-1}\ell\left(\hat{w}_{j|k},\hat{v}_{j|k}\right)\\[14.22636pt] \text{s.t.}\left\{\begin{array}[]{l}\begin{array}[]{rll}\hat{x}_{j+1|k}=&f\left(\hat{x}_{j|k},\hat{w}_{j|k}\right),&j\in\mathbb{Z}_{\left[0,k-1\right]}\\ y_{j}=&h\left(\hat{x}_{j|k}\right)+\hat{v}_{j|k},&\\ \end{array}\\ \hat{x}_{j|k}\in\mathcal{X},\ \hat{w}_{j|k}\in\mathcal{W},\ \hat{v}_{j|k}\in\mathcal{V},\end{array}\right.\end{array}$$

$${\Gamma }_{k-N|k}\left({\hat{x}}_{k-N|k}\right)=\Vert {\hat{x}}_{k-N|k}-{\bar{x}}_{k-N}{\Vert }_{P_{k-N|k}^{-1}},$$

$$\begin{split}\epsilon_{k-N}&=y_{k-N}-\hat{y}_{k-N|k},\\ N_{k}&=\left[1+\hat{x}_{k-N|k-1}^{T}\,P_{k-N-1}\hat{x}_{k-N|k-1}\right]\frac{\sigma}{|{\epsilon_{k-N}}|_{2}^{2}}\\ \alpha_{k}&=1-\frac{1}{N_{k}},\\ W_{k}&=\left[I-\frac{P_{k-N-1}\hat{x}_{k-N|k-1}\hat{x}_{k-N|k-1}^{T}}{1+\hat{x}_{k-N|k-1}^{T}P_{k-N-1}\hat{x}_{k-N|k-1}}\right]P_{k-N-1},\\ P_{k-N}&=\left\{\begin{array}[]{ccc}\frac{1}{\alpha_{k}}W_{k}&&\text{if }\frac{1}{\alpha_{k}}Tr\left({W_{k}}\right)\leq c,\\ W_{k}&&\text{otherwise},\end{array}\right.\end{split}$$

$${\bar{x}}_{k-N}={\hat{x}}_{k-N|k-1}.$$

$$\begin{array}[]{c}\underset{\hat{x}_{0|k},\boldsymbol{\hat{w}}}{\operatorname{min}}\,\Psi\coloneqq\sum\limits_{j=k-N}^{k-1}\ell\left(\hat{w}_{j|k},\hat{v}_{j|k}\right)+\\ \sum\limits_{j=1}^{k-N-1}\alpha_{k}^{k-N-j}\ell\left(\hat{w}_{j|k},\hat{v}_{j|k}\right)+\alpha_{k}^{k-N}\Gamma_{0|k}\left(\hat{x}_{0|k}\right)\\[14.22636pt] \text{s.t.}\left\{\begin{array}[]{l}\begin{array}[]{rll}\hat{x}_{j+1|k}=&f\left(\hat{x}_{j|k},\hat{w}_{j|k}\right),&j\in\mathbb{Z}_{\left[0,k-1\right]},\alpha\in(0,1]\\ y_{j}=&h\left(\hat{x}_{j|k}\right)+\hat{v}_{j|k},&\\ \end{array}\\ \hat{x}_{j|k}\in\mathcal{X},\ \hat{w}_{j|k}\in\mathcal{W},\ \hat{v}_{j|k}\in\mathcal{V},\end{array}\right.\end{array}$$

$$\begin{array}{rl}{\underline{\gamma }}_p\left(|{\hat{x}}_{k-N|k}-{\bar{x}}_{k-N}|\right)& \le {\Gamma }_{k-N}\left({\hat{x}}_{k-N|k}\right)\\ {\Gamma }_{k-N}\left({\hat{x}}_{k-N|k}\right)& \le {\bar{\gamma }}_p\left(|{\hat{x}}_{k-N|k}-{\bar{x}}_{k-N}|\right)\end{array}$$

$${\underline{\gamma }}_p\left(r\right)\ge {\underline{c}}_p{r}^a,{\bar{\gamma }}_p\left(r\right)\le {\bar{c}}_p{r}^a.$$

$$|P_0^{-1}|r^a\le {\Gamma }_{k-N}\left({\hat{x}}_{k-N|k}\right)\le |P_{\infty }^{-1}|r^a.$$

$$\begin{array}{rl}|x(k,z_1,{\mathbf{w}}_1)-x(k,z_2,{\mathbf{w}}_2)|& \le \max \{\beta \left(|z_1-z_2|,k\right),\\ & {\gamma }_1\left(\Vert {\mathbf{w}}_1-{\mathbf{w}}_2\Vert \right),\\ & {\gamma }_2\left(\Vert \mathbf{h}\left({\mathbf{x}}_1\right)-\mathbf{h}\left({\mathbf{x}}_2\right)\Vert \right)\}\\ & \le \beta (|z_1-z_2|,k)+\\ & {\gamma }_1\left(\Vert {\mathbf{w}}_1-{\mathbf{w}}_2\Vert \right)+\\ & {\gamma }_2\left(\Vert {\mathbf{v}}_1-{\mathbf{v}}_2\Vert \right)\end{array}$$

$$\beta (r,s)\le c_{\beta }{r}^p{s}^{-q}$$

$${\underline{\gamma }}_w\left(w\right)+{\underline{\gamma }}_v\left(v\right)\le \ell \left(w,v\right)\le {\bar{\gamma }}_w\left(w\right)+{\bar{\gamma }}_v\left(v\right).$$

$${\gamma }_1\left(3{\underline{\gamma }}_w^{-1}\left(r\right)\right)\le c_1{r}^{{\alpha }_1},{\gamma }_2\left(3{\underline{\gamma }}_v^{-1}\left(r\right)\right)\le c_2{r}^{{\alpha }_2}$$

$$\begin{array}{rl}|x_k-{\hat{x}}_k|& \le \Phi \left(|x_0-{\bar{x}}_0|,k\right)+{\pi }_w\left({\Vert \mathbf{w}\Vert }_{[0,k-1]}\right)+\\ & {\pi }_v\left({\Vert \mathbf{v}\Vert }_{[0,k-1]}\right).\end{array}$$

$$\begin{array}{rl}{\Psi }_N^{\ast }& =\Psi \left({\hat{x}}_{k-N|k}^{\ast },{\hat{\mathbf{w}}}^{\ast }\left[k-N,k-1\right]\right)\\ & ={\Gamma }_{k-N}\left({\hat{x}}_{k-N|k}^{\ast }\right)+\sum _{j=k-N}^{k-1}\ell \left({\hat{w}}_{j|k}^{\ast },{\hat{v}}_{j|k}^{\ast }\right),\end{array}$$

$${\Psi }_N^{\ast }$$

$$\begin{array}{rl}\Psi \left({\hat{x}}_{k-N|k}^{\ast },{\hat{\mathbf{w}}}^{\ast }\right)& \le \Psi \left(x_{k-N},\mathbf{w}\right),\\ & \le {\bar{\gamma }}_p\left(|x_{k-N}-{\bar{x}}_{k-N}|\right)+\\ & N{\bar{\gamma }}_w\left(\Vert \mathbf{w}\Vert \right)+N{\bar{\gamma }}_v\left(\Vert \mathbf{v}\Vert \right),\end{array}$$

$$\begin{array}{r}|{\hat{x}}_{k-N|k}-{\bar{x}}_{k-N}|\le {\underline{\gamma }}_p^{-1}({\bar{\gamma }}_p\left(|x_{k-N}-{\bar{x}}_{k-N}|\right)+\\ N{\bar{\gamma }}_w\left(\Vert \mathbf{w}\Vert \right)+N{\bar{\gamma }}_v\left(\Vert \mathbf{v}\Vert \right)).\end{array}$$

$$\begin{array}{rl}|{\hat{x}}_{k-N|k}-{\bar{x}}_{k-N}|& \le {\underline{\gamma }}_p^{-1}\left(3{\bar{\gamma }}_p\left(|x_{k-N}-{\bar{x}}_{k-N}|\right)\right)+\\ & {\underline{\gamma }}_p^{-1}\left(3N{\bar{\gamma }}_w\left(\Vert \mathbf{w}\Vert \right)\right)+{\underline{\gamma }}_p^{-1}\left(3N{\bar{\gamma }}_v\left(\Vert \mathbf{v}\Vert \right)\right),\\ & \le \frac{3^{\frac{1}{a}}}{|P_0^{-1}|}(|P_{\infty }^{-1}{|}^{\frac{1}{a}}|x_{k-N}-{\bar{x}}_{k-N}|+\\ & N^{\frac{1}{a}}{\bar{\gamma }}_w^{\frac{1}{a}}\left(\Vert \mathbf{w}\Vert \right)+N^{\frac{1}{a}}{\bar{\gamma }}_v^{\frac{1}{a}}\left(\Vert \mathbf{v}\Vert \right)).\end{array}$$

$$\begin{array}{rl}|{\hat{w}}_{j|k}|& \le {\underline{\gamma }}_w^{-1}\left(\frac{3}{N}{\bar{\gamma }}_p\left(|x_{k-N}-{\bar{x}}_{k-N}|\right)\right)+\\ & {\underline{\gamma }}_w^{-1}\left(3{\bar{\gamma }}_w\left(\Vert \mathbf{w}\Vert \right)\right)+{\underline{\gamma }}_w^{-1}\left(3{\bar{\gamma }}_v\left(\Vert \mathbf{v}\Vert \right)\right),\\ |{\hat{v}}_{j|k}|& \le {\underline{\gamma }}_v^{-1}\left(\frac{3}{N}{\bar{\gamma }}_p\left(|x_{k-N}-{\bar{x}}_{k-N}|\right)\right)+\\ & {\underline{\gamma }}_v^{-1}\left(3{\bar{\gamma }}_w\left(\Vert \mathbf{w}\Vert \right)\right)+{\underline{\gamma }}_v^{-1}\left(3{\bar{\gamma }}_v\left(\Vert \mathbf{v}\Vert \right)\right).\end{array}$$

$$\begin{array}{rl}|x_k-{\hat{x}}_{k|k}|& \le \beta \left(|x_{k-N}-{\hat{x}}_{k-N|k}|,N\right)+{\gamma }_1\left(\Vert {\mathbf{w}}_{\mathbf{j}}-{\hat{\mathbf{w}}}_{\mathbf{j}|\mathbf{k}}\Vert \right)\\ & +{\gamma }_2\left(\Vert {\mathbf{v}}_{\mathbf{j}}-{\hat{\mathbf{v}}}_{\mathbf{j}|\mathbf{k}}\Vert \right).\end{array}$$

$$\begin{array}{r}\beta \left(|x_{k-N}-{\hat{x}}_{k-N|k}|,N\right)\le \beta \left(2|x_{k-N}-{\bar{x}}_{k-N}|,N\right)+\\ \beta \left(2|{\hat{x}}_{k-N|k}-{\bar{x}}_{k-N}|,N\right)\\ \le \beta \left(2|x_{k-N}-{\bar{x}}_{k-N}|,N\right)+\\ \beta (\frac{23^{\frac{1}{a}}|P_{\infty }^{-1}{|}^{\frac{1}{a}}}{|P_0^{-1}|}|x_{k-N}-{\bar{x}}_{k-N}|+\frac{23^{\frac{1}{a}}{N}^{\frac{1}{a}}}{|P_0^{-1}|}{\bar{\gamma }}_w^{\frac{1}{a}}\left(\Vert \mathbf{w}\Vert \right)+\\ \frac{23^{\frac{1}{a}}{N}^{\frac{1}{a}}}{|P_0^{-1}|}{\bar{\gamma }}_v^{\frac{1}{a}}\left(\Vert \mathbf{v}\Vert \right),N)\\ \le \beta \left(2|x_{k-N}-{\bar{x}}_{k-N}|,N\right)+\\ \beta \left(\frac{63^{\frac{1}{a}}|P_{\infty }^{-1}{|}^{\frac{1}{a}}}{|P_0^{-1}|}|x_{k-N}-{\bar{x}}_{k-N}|,N\right)+\\ \beta \left(\frac{63^{\frac{1}{a}}{N}^{\frac{1}{a}}}{|P_0^{-1}|}{\bar{\gamma }}_w^{\frac{1}{a}}\left(\Vert \mathbf{w}\Vert \right),N\right)+\beta \left(\frac{63^{\frac{1}{a}}{N}^{\frac{1}{a}}}{|P_0^{-1}|}{\bar{\gamma }}_v^{\frac{1}{a}}\left(\Vert \mathbf{v}\Vert \right),N\right).\end{array}$$

$$\begin{array}{r}\beta \left(|x_{k-N}-{\hat{x}}_{k-N|k}|,N\right)\le \frac{c_{\beta }{2}^p}{N^q}|x_{k-N}-{\bar{x}}_{k-N}{|}^p+\\ \frac{c_{\beta }{6}^p{3}^{\frac{p}{a}}|P_{\infty }^{-1}{|}^{\frac{p}{a}}}{|P_0^{-1}{|}^p{N}^q}|x_{k-N}-{\bar{x}}_{k-N}{|}^p+\\ \frac{c_{\beta }{6}^p{3}^{\frac{p}{a}}{N}^{\frac{p}{a}}}{|P_0^{-1}|}{\bar{\gamma }}_w^{\frac{p}{a}}\left(\Vert \mathbf{w}\Vert \right)+\frac{c_{\beta }{6}^p{3}^{\frac{p}{a}}{N}^{\frac{p}{a}}}{|P_0^{-1}|}{\bar{\gamma }}_v^{\frac{p}{a}}\left(\Vert \mathbf{v}\Vert \right)\\ \le \frac{c_{\beta }{2}^p}{N^q}|x_{k-N}-{\bar{x}}_{k-N}{|}^p+\\ {\left(\frac{|P_{\infty }^{-1}|}{|P_0^{-1}|}\right)}^p\frac{c_{\beta }{6}^p{3}^{\frac{p}{a}}}{N^q}|x_{k-N}-{\bar{x}}_{k-N}{|}^p+\\ \frac{c_{\beta }{6}^p{3}^{\frac{p}{a}}{N}^{\frac{p}{a}-q}}{|P_0^{-1}|}{\bar{\gamma }}_w^{\frac{p}{a}}\left(\Vert \mathbf{w}\Vert \right)+\frac{c_{\beta }{6}^p{3}^{\frac{p}{a}}{N}^{\frac{p}{a}-q}}{|P_0^{-1}|}{\bar{\gamma }}_v^{\frac{p}{a}}\left(\Vert \mathbf{v}\Vert \right).\end{array}$$

$$\begin{array}{r}\beta \left(|x_{k-N}-{\hat{x}}_{k-N|k}|,N\right)\le \frac{c_{\beta }18^p}{|P_0^{-1}|}\left({\bar{\gamma }}_w^{\frac{p}{a}}\left(\Vert \mathbf{w}\Vert \right)+{\bar{\gamma }}_v^{\frac{p}{a}}\left(\Vert \mathbf{v}\Vert \right)\right)+\\ \frac{c_{\beta }}{N^q}\left(2^p+{\mathbb{C}}_{P^{-1}}^{p}18^p\right)|x_{k-N}-{\bar{x}}_{k-N}{|}^p.\end{array}$$

$$\begin{array}{rl}{\gamma }_1\left(\Vert {\mathbf{w}}_{\mathbf{j}}-{\hat{\mathbf{w}}}_{\mathbf{j}|\mathbf{k}}\Vert \right)& \le {\gamma }_1\left(\Vert \mathbf{w}\Vert +\Vert {\hat{\mathbf{w}}}_{\mathbf{j}|\mathbf{k}}\Vert \right)\\ & \le {\gamma }_1(\Vert \mathbf{w}\Vert +{\underline{\gamma }}_w^{-1}\left(\frac{3}{N}{\bar{\gamma }}_p\left(|x_{k-N}-{\bar{x}}_{k-N}|\right)\right)+\\ & {\underline{\gamma }}_w^{-1}\left(3{\bar{\gamma }}_w\left(\Vert \mathbf{w}\Vert \right)\right)+{\underline{\gamma }}_w^{-1}\left(3{\bar{\gamma }}_v\left(\Vert \mathbf{v}\Vert \right)\right)).\end{array}$$

$$\begin{array}{r}{\gamma }_1\left(\Vert {\mathbf{w}}_{\mathbf{j}}-{\hat{\mathbf{w}}}_{\mathbf{j}|\mathbf{k}}\Vert \right)\le \frac{c_1{3}^{{\alpha }_1}|P_{\infty }^{-1}{|}^{{\alpha }_1}}{N^{{\alpha }_1}}|x_{k-N}-{\bar{x}}_{k-N}{|}^{a{\alpha }_1}+\\ c_1{3}^{{\alpha }_1}{\bar{\gamma }}_v^{{\alpha }_1}\left(\Vert \mathbf{v}\Vert \right)+{\gamma }_1\left(3\left(\Vert \mathbf{w}\Vert +{\underline{\gamma }}_w^{-1}\left(3{\bar{\gamma }}_w\left(\Vert \mathbf{w}\Vert \right)\right)\right)\right).\end{array}$$

$$\begin{array}{r}{\gamma }_2\left(\Vert {\mathbf{v}}_{\mathbf{j}}-{\hat{\mathbf{v}}}_{\mathbf{j}|\mathbf{k}}\Vert \right)\le \frac{c_2{3}^{{\alpha }_2}|P_{\infty }^{-1}{|}^{{\alpha }_2}}{N^{{\alpha }_2}}|x_{k-N}-{\bar{x}}_{k-N}{|}^{a{\alpha }_2}+\\ c_2{3}^{{\alpha }_2}{\bar{\gamma }}_w^{{\alpha }_2}\left(\Vert \mathbf{w}\Vert \right)+{\gamma }_2\left(3\left(\Vert \mathbf{v}\Vert +{\underline{\gamma }}_v^{-1}\left(3{\bar{\gamma }}_v\left(\Vert \mathbf{v}\Vert \right)\right)\right)\right).\end{array}$$

Full text

1 Introduction

State estimation plays a fundamental role in feedback control, system monitoring and system optimization because noisy measurements is the only information available from the system. Several methods have been developed for accomplishing such task (see jazwinski2007stochastic ; crassidis2004optimal ; among others). All these methods have been developed upon the assumption on the knowledge of noises and model of the system, as well as, the absence of constraints.

In practice, these assumptions are not easily satisfied and research efforts were focused on approaches that do not relay on such requirements (see li1997linear , sayed2001framework , blanchini2008set , among others). For example, an $H_{\infty}$ filter is designed minimizing the $H_{\infty}$ norm of the mapping between disturbances and estimation error. In li1997linear , el1997robust and hu2009improved an approach that solves a least-square estimation problem is introduced. Both methods are based on the adequate selection of the uncertainty model instead of relying on statistical assumptions on noises. In these approaches, uncertainty models are formulated based on the available information of the system. In the same way, robust estimation algorithms based on as min-max robust filtering, set-valued estimation and guaranteed cost paradigm, have attracted the attention of the research community (see sayed2001framework , zhu2002design ).

Building on the success on moving horizon control, moving horizon estimation (MHE) has attracted attention of researchers since the pioneering work of jazwinski1968limited (see also schweppe1973uncertain , rao2001constrained and rao2003constrained ). The interest in such estimation methods stems from the possibility of dealing with limited amount of data, instead of using all the information available from the beginning, and the ability to incorporate constraints. In recent years, both theoretical properties of various MHE schemes as well as efficient computational methods for real-time implementation have been studied (see alessandri2005robust , alessandri2008moving , alessandri2012min , garcia2016new , sartipizadeh2016computationally , sanchez2017adaptive ). In particular, it is of interest to establish robust stability and estimate convergence properties. In recent years several results have been obtained for different algorithms, advancing from idealistic assumptions (observability and no disturbances) to realistic situations (detectability and bounded disturbances).

For nonlinear observable systems, rao2003constrained established the asymptotic stability of the estimation error for the standard cost function. Furthermore, if the disturbances are asymptotically vanishing the estimation error is robust asymptotically stable and it asymptotically converges to zero (rawlings2009model -rawlings2012optimization ). alessandri2008moving and alessandri2010advances proposed an estimation scheme, based on least-square cost function of the estimation residuals, that guaranteed the boundedness of estimation error for observable systems subject to bounded additive disturbances. Finally, for the general case of nonlinear detectable systems subject to bounded disturbances, ji2016robust and muller2017nonlinear showed the robust global asymptotic stability (RGAS) and convergence of estimation error in case of bounded or vanishing disturbances, respectively. In these works, the least-square objective function was modified by adding a max-term. ji2016robust established RGAS for the full information estimator while muller2017nonlinear established RGAS and convergence for the moving horizon estimator. Furthermore, for a particular choice of the weights of the objective function, muller2017nonlinear established these results for the least-squares type objective function.

This paper introduces the RGAS and convergence analysis for the moving horizon estimator based on adaptive arrival cost proposed in sanchez2017adaptive in the practical case of nonlinear detectable systems subject to bounded disturbances. To establish robust stability properties for MHE it is crucial that the prior weighting in the cost function is chosen properly. In various schemes the necessary assumptions in the prior weighting are difficult to verify (rao2003constrained , rawlings2009model ), while in others can be verified a prior muller2017nonlinear . In the MHE scheme analysed in this work, the assumption on the prior weighting can be verified a prior by design. Furthermore, the disturbances gains become uniform (i.e., they are valid independent of $N$), allowing to extend the stability analysis to full information estimators with least-square type cost functions.

The rest of the paper is organized as follows: Section 2 introduces the notation, definitions and properties that will be used through the paper. Section 3 presents the main result and shows its connections with previous stability analysis. Section 4 discusses simple examples, previously used in the literature, with the purpose of illustrating the concepts and also in order to show the difference with others MHE algorithms. Finally, Section 5 presents conclusions.

2 Preliminaries and setup

2.1 Notation

Let $\mathbb{Z}_{\left[a,b\right]}$ denotes the set of integers in the interval $\left[a,b\right]\subseteq\mathbb{R}\textnormal{, and }\mathbb{Z}_{\geq a}$ denotes the set of integers greater or equal to $a$. Boldface symbols denote sequences of finite or infinite length, i.e., $\boldsymbol{w}\coloneqq\{w_{k_{1}},\ldots,w_{k_{2}}\}\textnormal{ for some }k_{1},k_{2}\in\mathbb{Z}_{\geq 0}\textnormal{ and }k_{1}<k_{2}$, respectively. We denote $x_{j|k}$ as the finite sequence $\boldsymbol{x}$ given at time $k\in\mathbb{Z}_{\geq 0}\textnormal{ and }j\in\left[k_{1},k_{2}\right]$. By $|{x}|$ we denote the Euclidean norm of a vector $x\in\mathbb{R}^{n}$. Let $\left\lVert\boldsymbol{x}\right\rVert\coloneqq\sup_{k\in\mathbb{Z}_{\geq 0}}|{x_{k}}|$ denote the supreme norm of the sequence $\boldsymbol{x}\textnormal{ and }\left\lVert\boldsymbol{x}\right\rVert_{\left[a,b\right]}\coloneqq\sup_{k\in\mathbb{Z}_{\left[a,b\right]}}|{x_{k}}|$ . A function $\gamma:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is of class $\mathcal{K}$ if $\gamma$ is continuous, strictly increasing and $\gamma\left(0\right)=0$ . If $\gamma$ is also unbounded, it is of class $\mathcal{K}_{\infty}$. A function $\zeta:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is of class $\mathcal{L}$ if $\zeta\left(k\right)$ is non increasing and $\lim_{k\rightarrow\infty}\zeta\left(k\right)=0$. A function $\beta:\mathbb{R}_{\geq 0}\times\mathbb{Z}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is of class $\mathcal{KL}$ if $\beta\left(\cdot,k\right)$ is of class $\mathcal{K}$ for each fixed $k\in\mathbb{Z}_{\geq 0}$, and $\beta\left(r,\cdot\right)$ of class $\mathcal{L}$ for each fixed $r\in\mathbb{R}_{\geq 0}$.

The following inequalities hold for all $\beta\in\mathcal{KL},\;\gamma\in\mathcal{K}\textnormal{ and }a_{j}\in\mathbb{R}_{\geq 0}\textnormal{ with }j\in\mathbb{Z}_{\left[1,n\right]}$

[TABLE]

The preceding inequalities hold since $\max\{a_{j}\}$ is included in the sequence $\{a_{1},a_{2},\ldots,a_{n}\}$ and $\mathcal{K}$ functions are non-negative strictly increasing functions.

Bounded sequences: A sequence $\boldsymbol{w}$ is bounded if $\left\lVert\boldsymbol{w}\right\rVert$ is finite. The set of bounded sequences $\boldsymbol{w}$ is denoted as $\mathcal{W}\left(w_{\max}\right)\coloneqq\{w:\left\lVert\boldsymbol{w}\right\rVert\leq w_{\max}\}$ for some $w_{\max}\in\mathbb{R}_{\geq 0}$

Convergent sequences: A bounded infinite sequence $\boldsymbol{w}$ is convergent if $|{w_{k}}|\rightarrow 0$ as $k\rightarrow\infty$. Let denote the set of convergent sequences $\mathcal{C}$:

[TABLE]

Analogously, $\mathcal{C}_{v}$ is defined for the sequence $\boldsymbol{v}$.

2.2 Problem statement

Let us consider the state estimation problem for nonlinear discrete time systems of the form

[TABLE]

where $x_{k}\in\mathcal{X}\subseteq\mathbb{R}^{n},w_{k}\in\mathcal{W}\subseteq\mathbb{R}^{p},y_{k}\in\mathcal{Y}\subseteq\mathbb{R}^{m},v_{k}\in\mathcal{V}\subseteq\mathbb{R}^{m}$ are the state, process noise, measurement and estimation residuals vectors, respectively. The process disturbance $w_{k}$ and estimation residuals $v_{k}$ are unknown but assumed to be bounded, i.e, $\boldsymbol{w}\in\mathcal{W}\left(w_{\max}\right)\textnormal{ and }\boldsymbol{v}\in\mathcal{V}\left(v_{\max}\right)$ for some $w_{\max},v_{\max}\in\mathbb{R}_{\geq 0}$. $\mathcal{X},\mathcal{Y},\mathcal{W}\textnormal{ and }\mathcal{V}$ are compact and convex sets with the null vector $\mathbf{0}$ belongs to them. In the following we assume that $f:\mathbb{R}^{n}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{n}$ is continuous, locally Lipschitz on $x_{k}$ and $h:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is continuous. The solution to the system (2) at time $k$ is denoted by $x\left(k;\mathtt{x}_{0},\boldsymbol{w}\right)$, with initial condition $\mathtt{x}_{0}$ and process disturbance sequence $\boldsymbol{w}$. Furthermore, the initial condition $\mathtt{x}_{0}$ is unknown, but a prior knowledge $\bar{x}_{0}$ is assumed to be available and its error is assumed to be bounded, i.e., $\bar{x}_{0}\in\mathcal{X}_{0}\coloneqq\{\bar{x}_{0}:|{\mathtt{x}_{0}-\bar{x}_{0}}|\leq e_{\max}\}$, $\mathcal{X}_{0}\subseteq\mathcal{X}$.

The solution of the estimation problem aims to find at time $k$ an estimate $\hat{x}_{k|k}$ of the current state $x_{k}$ minimizing a performance metric using by the MHE. At each sampling time $k$, given the previous $N$ measurements $\boldsymbol{y}:=\{y_{k-N},\ldots,y_{k-1}\}$, the following optimization problem is solved

[TABLE]

where $\hat{x}_{k-j|k}$ is the optimal estimated and $\hat{w}_{j|k}$ is the optimal process noise estimate at sample $k-j\quad j=0,1,\ldots,N$ based on measurements $y_{k-j}$ available at time $k$. The process noise $\boldsymbol{\hat{w}}_{j|k}\coloneqq\{\hat{w}_{k-N|k},\ldots,\hat{w}_{k-1|k}\}$ and $\hat{x}_{k-N|k}$ are the optimization variables. The stage cost $\ell\left(w_{j|k},v_{j|k}\right)$ penalizes the estimated process noise sequence $\boldsymbol{\hat{w}}_{j|k}$ and the estimation residuals $\boldsymbol{\hat{v}}_{j|k}=\boldsymbol{y}_{j}-h\left(\boldsymbol{\hat{x}}_{j|k}\right)$, while $\Gamma_{k-N}\left(\hat{x}_{k-N|k}\right)$ penalizes the prior estimated $\hat{x}_{k-N|k}$. The adequate choice of $\ell\left(\cdot\right)$ and $\Gamma_{k-N}\left(\cdot\right)$, and their parameters, allows to ensure the robust stability of the estimator muller2017nonlinear . While the estimation window is not full, $k\leq N$, problem (3) can be reformulated and solved as a full information problem

[TABLE]

as $k$ increases this problem becomes (3) for all $k\geq N$.

In previous works, the robust stability of MHE has been achieved by modifying the standard least-square cost function through the inclusion of a $max$–term (ji2016robust ; muller2017nonlinear ) or by a suitable choice of the cost’s function parameters (muller2017nonlinear ). Another mechanism to solve this problem is combining a suitable choice of the stage cost $\ell\left(\hat{w}_{j|k},\hat{v}_{j|k}\right)$ with a time–varying prior weight of the form

[TABLE]

whose parameters $\left(P^{-1}_{k-N|k},\bar{x}_{k-N}\right)$ are recursively updated using the information available at time $k$ (sanchez2017adaptive , 2017). The prior weighting is defined in this way to avoid the introduction of artificial cycling in the estimation process (see rawlings2009model ). In this approach, the prior weight matrix $P_{k-N|k}$ is given by

[TABLE]

where $\sigma,\,\sigma_{w},\,c,\,\lambda\in R_{>0},\,c>\lambda,\,P_{0}=\lambda I_{n\times n}$ and $\sigma\gg\sigma_{w}$, where $\sigma_{w}$ denotes the process noise variance. The prior knowledge of the window $\bar{x}_{k-N}$ is updated using a smoothed estimate (findeisen1997moving )

[TABLE]

The optimization problem (3) can be reformulated in terms of the initial condition $\hat{x}_{0}$ and the estimated process noises and the residuals along the entire trajectory as follows

[TABLE]

This formulation of problem (3) allows to explicitly see the effect of past data on the current state estimate $\hat{x}_{k|k}$. In this formulation it is easy to see the exponential averaging of these data. Allowing $\alpha$ change in time, the past data has different affects on the current estimates depending on $\hat{x}_{k|k}$.

Before proceeding to the development of the main results, we state the main properties and assumptions about the prior weighting $\Gamma_{k-N}$.

The updating mechanism (5) is a time-varying filter whose inputs are $\hat{x}_{k-N|k-1}\hat{x}_{k-N|k-1}^{T}$ and the initial condition $P_{0}$. It generates recursively a real-time estimation of $P_{k-N|k}$ by updating $P_{k-N-1|k-1}$ with an exponential time-averaging of $\hat{x}_{k-N|k-1}\hat{x}_{k-N|k-1}^{T}$. The updating mechanism (5) only use data and it does not rely on a model of the system. The sequence $P_{k|k}\quad k\geq 0$ is positive definite, it is decreasing in norm and it is bounded. The proof of these properties follows similar steps as in sanchez2017adaptive .

Assumption 1**.**

The prior weighting $\Gamma_{k-N}$ is a continuous function $\Gamma_{k-N}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ lower bounded by $\underline{\gamma}_{p}\in\mathcal{K}_{\infty}{}$ and upper bounded by $\bar{\gamma}_{p}\in\mathcal{K}_{\infty}$ such that:

[TABLE]

for all $\hat{x}\in\mathcal{X}$ and

[TABLE]

where $0\leq\underline{c}_{p}\leq\bar{c}_{p}$ and $a\in R_{\geq 1}$.

Given prior weighting updating scheme (5) inequality (7) satisfies [sanchez2017adaptive, ]

[TABLE]

Definition 1**.**

The system (2) is incrementally input/output-to-state stable if there exist functions $\beta\in\mathcal{KL}$ and $\gamma_{1},\gamma_{2}\in\mathcal{K}$ such that for every two initial states $z_{1}$, $z_{2}\in\mathbb{R}^{n}$, and any two disturbances sequences $\boldsymbol{w_{1}},\boldsymbol{w_{2}}$ the following holds for all $k\in\mathbb{Z}_{\geq 0}$:

[TABLE]

This definition combines the concepts of output-to-state-stability (OSS) and input-to-state-stability (ISS). As stated in sontag1997output , the notion of IOSS represents a natural combination of the ideas of strong observability and ISS, and it was called detectability in sontag1989some and strong unboundedness observability in jiang1994small . In addition, the existence of an observer for the system (2), which is incrementally input-output-to-state stable (i-IOSS) instead of IOSS (see Remark 24 in sontag1997output ), is assumed. Note that $\left\lVert\boldsymbol{h\left(x_{1}\right)-h\left(x_{2}\right)}\right\rVert=\left\lVert\boldsymbol{v_{1}-v_{2}}\right\rVert$, since $y_{k}=h\left(x_{k}\right)+v_{k}$ These assumptions will help us to bound the functions involved in the definition of i-IOSS and to relate them with the terms of the MHE cost function (stage cost and prior weight).

In the following sections the updating mechanism (5) and the assumption of i-IOSS sontag2008input will be used to prove robust stability of the proposed MHE in the presence of bounded disturbances and convergence to the true state in the case of convergent disturbances. Some assumptions about functions related to system (2) and Definition 1 will be helpful in the sequel.

Assumption 2**.**

The function $\beta(r,s)\in\mathcal{KL}$ and satisfies the following inequality

[TABLE]

for some $c_{\beta}\in\mathbb{R}_{\geq 0}$, $p\in\mathbb{R}_{\geq 0}$ and $q\in\mathbb{R}_{\geq 0}$ and $q\geq p$.

Assumption 3**.**

The stage cost $\ell\left(\cdot,\cdot\right):\mathbb{R}^{p}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a continuous function bounded by $\underline{\gamma}_{w},\underline{\gamma}_{v},\bar{\gamma}_{w},\bar{\gamma}_{v}$ $\in\mathcal{K}_{\infty}{}$ such that the following inequalities are satisfied $\forall w\in\mathcal{W}\textnormal{ and }v\in\mathcal{V}$

[TABLE]

Functions $\gamma_{1}$ and $\gamma_{2}$ from Definition 1 are related with the bounds of stage cost $\bar{\gamma}_{w},\underline{\gamma}_{w},\bar{\gamma}_{v}$ and $\underline{\gamma}_{v}$ through the following inequalities

[TABLE]

for $c_{1},c_{2},\alpha_{1},\alpha_{2}>0.$ Inequalities (11) to (13) were used in previous works (ji2016robust ; muller2017nonlinear ).

In this work, we claim that the proposed estimator holds the property of being robust global asymptotic stable, which is defined as follows.

Definition 2**.**

Consider the system described (2) subject to disturbances $\boldsymbol{w}\in\mathcal{W}\left(w_{\max}\right)$ and $\boldsymbol{v}\in\mathcal{V}\left(v_{\max}\right)$ for $w_{\max}\in\mathbb{R}_{\geq 0}$, $v_{max}\in\mathbb{R}_{\geq 0}$ with prior estimate $\bar{x}_{0}\in\mathcal{X}\left(e_{\max}\right)$ for $e_{\max}\in\mathbb{R}_{\geq 0}$. The moving horizon state estimator given by equation (3) with adaptive prior weight is robustly globally asymptotically stable (RGAS) if there exists functions $\Phi\in\mathcal{KL}$ and $\pi_{w}$, $\pi_{v}\in\mathcal{K}$ such that for all ${x}_{0}\in\mathcal{X}$, all $\bar{x}_{0}\in\mathcal{X}_{0}$, the following is satisfied for all $k\in\mathbb{Z}_{\geq 0}$

[TABLE]

We want to show that if system (2) is i-IOSS, then Assumptions 1, 2 and 3 are fulfilled and the proposed MHE estimator with adaptive arrival cost weight matrix is RGAS. Furthermore, if the process disturbance and measurement noise sequences are convergent (i.e., $\boldsymbol{w}\in\mathcal{C}_{w},\boldsymbol{v}\in\mathcal{C}_{v}$), then $\hat{x}_{k|k}\rightarrow x_{k}$ as $k\rightarrow\infty$.

3 Robust stability of moving horizon estimation under bounded disturbances

We are ready to derive the main result: RGAS of the proposed moving horizon estimator with a large enough estimation horizon $\mathcal{N}$ for nonlinear detectable systems under bounded disturbances. Furthermore, a $\mathcal{K}\mathcal{L}$ function exist such that (14) is valid with this $\Phi,\pi_{w}$ and $\pi_{v}$ for all estimation horizon $N\geq\mathcal{N}$.

Theorem 1**.**

Consider an i-IOSS system (2) with disturbances $\boldsymbol{w}\in\mathcal{W}\left(w_{\max}\right)$, $\boldsymbol{v}\in\mathcal{V}\left(v_{\max}\right)$. Assume that the arrival cost weight matrix of the MHE problem $\Gamma_{k-N}$ is updated using the adaptive algorithm (5). Moreover, Assumptions 1, 2 and 3 are fulfilled and initial condition $x_{0}$ is unknown, but a prior estimate $\bar{x}_{0}\in\mathcal{X}_{0}$ is available. Then, the MHE estimator (3) is $RGAS$.

Proof. The optimal cost of problem (3) is given by

[TABLE]

which is bounded (Assumptions 1 and 3) $\forall\,|{\hat{w}_{j|k}}|$ and $\forall\,|{\hat{v}_{j|k}}|$ for all $j\in\mathbb{Z}_{\left[k-N,k-1\right]}$ by

[TABLE]

Due optimality, the following inequalities hold $\forall k\in[k-N,k-1]$

[TABLE]

then, taking into account the lower and upper bounds we have

[TABLE]

By mean of Assumptions 1 and 3 the last inequality can be written as follows

[TABLE]

Analogously, bounds for $|{\hat{v}_{j|k}}|$ and $|{\hat{w}_{j|k}}|$ can be found

[TABLE]

Next, let us consider some sample $k\in\mathbb{Z}_{\geq N}$. Assuming that system (2) is i-IOSS with $z_{1}=x_{k-N},z_{2}=\hat{x}_{k-N|k},w_{1}=\{w_{j}\},w_{2}=\{\hat{w}_{j|k}\},v_{1}=\{v_{j}\}$ and $v_{2}=\{\hat{v}_{j|k}\}$ for all $j\in\mathbb{Z}_{\left[k-N,k-1\right]}$. Since $x(k)=x\left(N,z_{1},\boldsymbol{w_{1}}\right),\hat{x}(k)=\hat{x}_{k|k}=x\left(N,z_{2},\boldsymbol{w_{2}}\right)$ we obtain

[TABLE]

In order to get a finite upper bound for the estimation error, the three terms in the right hand side of equation (17) must be upper bounded. The first term can be written

[TABLE]

Using Assumptions 1 and 2, function $\beta(\cdot)$ is bounded by

[TABLE]

Taking in account that $P_{k}^{-1}$ is a symmetric positive definite matrix for all $k\in\mathbb{Z}_{\left[0,\infty\right)}$, then $|{P_{k}^{-1}}|\leq\lambda_{\max}\left(P_{k}^{-1}\right)$, where $\lambda_{\max}\left(P_{k}^{-1}\right)$ denotes the maximal eigenvalue of matrix $P_{k}^{-1}$. Denoting $\lambda_{\min}\left(P_{k}^{-1}\right)$ as the minimal eigenvalue of matrix $P_{k}^{-1}$ and taking in account that $|{P_{k}^{-1}}|\leq|{P_{k+1}^{-1}}|$, the maximum conditioning number of matrix $P_{k}^{-1}$ can be defined as $\mathbb{C}_{P^{-1}}\coloneqq\lambda_{\max}\left(P_{\infty}^{-1}\right)/\lambda_{\min}\left(P_{0}^{-1}\right)$, then $\beta\left(|{x_{k-N}-\hat{x}_{k-N|k}}|,\;N\right)$ can be bounded by

[TABLE]

The first term in the right side of this equation is bounded due the assumption that $|{x_{k-N}-\bar{x}_{k-N}}|\in\mathcal{X}_{0}\left(e_{\max}\right)$, while the second term are finite constants. To extend the validness of (18) to the full estimation horizon, an extension of the function $\beta$ at the beginning of the estimation, $N=0$, is required.

The second term in the right hand side of equation (17), can be bounded by the following inequality $\forall j\in\mathbb{Z}_{\left[k-N,k-1\right]}$

[TABLE]

Recalling Assumption 3, the reader can verify the following inequality

[TABLE]

In an equivalent manner, a bound for the third term in the right hand side of equation (17) can be found

[TABLE]

Once an upper bound for the three terms of equation (17) were found, defining $\zeta\coloneqq\max\{p,a\alpha_{1},a\alpha_{2}\},\eta\coloneqq\min\{q,\alpha_{1},\alpha_{2}\}$ and $\rho\coloneqq\max\{p,\alpha_{1},\alpha_{2}\}$, equation (17) can be rewritten as follows

[TABLE]

Defining the functions $\bar{\beta}\left(r,\;s\right),\phi_{w}\left(r\right)$ and $\phi_{v}\left(r\right)$ for all $r\geq 0$ and $s\;\in\mathbb{Z}_{\geq 1}$ as follows

[TABLE]

equation (21) can be written $\forall k\in\mathbb{Z}_{\left[1,N-1\right]}$ as follows

[TABLE]

To guarantee the validity of previous results on the entire time horizon we must extend the definition of $\beta\left(r,s\right)\text{ to }s=0$. Because of $\bar{\beta}\left(r,s\right)\in\mathcal{KL}$, $\bar{\beta}\left(r,0\right)\in\mathcal{KL}$ and $\bar{\beta}\left(r,0\right)\geq\bar{\beta}\left(r,k\right)$ for $k\in\mathbb{Z}_{\geq 1}$, it is sufficient to define $\bar{\beta}\left(r,0\right)\geq k_{\beta}\;\bar{\beta}\left(r,1\right)$ for some $k_{\beta}\in\mathbb{R}_{>1}$ to extend the definition of $\bar{\beta}\left(r,s\right)\text{for all }k\in\mathbb{Z}_{\geq 0}$. We would like to determinate the decreasing rate for the function $\bar{\beta}\left(r,s\right)$ $\mathcal{N}$ samplings time in the future. In order to do that, let define the constants

[TABLE]

and

[TABLE]

The minimum horizon length required to accomplish a decreasing rate $\delta$ will be given by

[TABLE]

Adopting an estimator with a window length greater or equal to $\mathcal{N}$ such that

[TABLE]

the effects of the initial conditions will vanish with a decreasing rate $\delta$. As $k\rightarrow\infty$, the estimation will entry to the bounded set $\mathcal{X}\left(w,v\right)\in\mathcal{X}$ defined by the noises of the system

[TABLE]

This set define the minimum size region of error space $\mathcal{X}$ that the error can achieve by removing the effect of errors in initial conditions ($e_{max}$). Equation (27) establish a trade off between speed of convergence and window length, which is related with the size of $\mathcal{X}\left(w,v\right)$.

For any MHE with adaptive arrival cost and window length $N\geq\mathcal{N}$ two situations can be considered

• •

The estimator removed the effects of $x_{0}$ on $\hat{x}_{k+j|k+j}$ such that $x_{k+j}-\hat{x}_{k+j|k+j}\in\mathcal{X}\left(w,v\right)$, and

• •

The estimator has not removed the effects of $x_{0}$ on $\hat{x}_{k+j|k+j}$ such that $x_{k+j}-\hat{x}_{k+j|k+j}\notin\mathcal{X}\left(w,v\right)$,

Assuming the first situation and recalling equations (25) and (27), the following inequalities hold

[TABLE]

This equation implies the fact that the estimation error $x_{k+j}-\hat{x}_{k+j|k+j}\in\mathcal{X}\left(w,v\right)\quad\forall j\in\mathbb{Z}_{\geq 0}$.

In the other case, when the estimation error is outside of $\mathcal{X}\left(w,v\right)$, equations (25) and (27) are recalled again and the following inequalities hold

[TABLE]

Since $\delta>\frac{2+\mu}{1+\mu}$, then $\forall N\geq\mathcal{N}$ we have

[TABLE]

Equations (30) and (31) reveal a contractive behaviour of the estimation error with $\theta$ as contraction factor. For some finite time $k^{*}$ the estimation error will decrease until $x_{k^{*}+j}-\hat{x}_{k^{*}+j|k^{*}+j}\in\mathcal{X}\left(w,v\right)$.

In an equivalent formulation, equations (29) and (30) put in evidence the existence of a positive invariant set and a Lyapunov like function for the proposed estimator. From equation (30), one can see that for the case that the estimation error belong to the set $\mathcal{X}\left(w,v\right)^{C}\cap\mathcal{X}$, the estimation error decreases in a factor of $\theta$ every $\mathcal{N}$ sampling time. Taking in account the general case in which $|{x_{k}-\hat{x}_{k|k}}|\in\mathcal{X}$ for $k\in\mathbb{Z}_{\geq\mathcal{N}}$, following the same procedure as in muller2017nonlinear , we could define $i\coloneqq\lfloor{\frac{k}{N}}\rfloor$ (where $\lfloor\cdot\rfloor$ denotes the floor function) and $j\coloneqq k\;\textit{mod}\;N$, therefore $k=iN+j$. Combining equations (29) and (30) and the fact that $|{x_{j}-\hat{x}_{j}}|\leq\delta r_{\max}$ for $j\in\mathbb{Z}_{\left[0,N-1\right]}$ one can obtain

[TABLE]

where

[TABLE]

Since $\bar{\beta}\left(r,s\right)\in\mathcal{KL}$, function $\Phi\left(|{x_{0}-\bar{x}_{0}}|,\;k\right)$ could increase in the steps $\mathbb{Z}_{\left[iN-1,iN\right]}$ for $i\geq 1$ (recall definition in Equation (22)). Therefore, define $\bar{\Phi}\left(|{x_{0}-\bar{x}_{0}}|,\;k\right)$ which is an upper bound for $\Phi\left(|{x_{0}-\bar{x}_{0}}|,\;k\right)$. Taking in account that noises at time $\geq k$ do not affect the estimation at time $k$, equation (32) can be rewritten as

[TABLE]

This equation is just equation (14) with

[TABLE]

therefore the estimator proposed in equations (3) is RGAS.

Finally, in order to prove that the estimation error $|{x_{k}-\hat{x}_{k|k}}|\rightarrow 0$ when $|{w_{k}}|\in\mathcal{C}_{w},|{v_{k}}|\in\mathcal{C}_{v}$, we must note that equation (25) holds for $\left\lVert\boldsymbol{w}\right\rVert_{\left[k-N,k-1\right]}$ instead of $\left\lVert\boldsymbol{w}\right\rVert$ and $\left\lVert\boldsymbol{v}\right\rVert_{\left[k-N,k-1\right]}$ instead $\left\lVert\boldsymbol{v}\right\rVert$ (it can be done omitting last step in equation (15)). From a qualitative point of view, taking in account that function $\bar{\Phi}\left(|{x_{0}-\bar{x}_{0}}|,\;k\right)\in\mathcal{KL}$ and sequences $\boldsymbol{w}$ and $\boldsymbol{v}$ are convergent, the right hand side of equation (33) tends to zero as $k\rightarrow\infty$. $\square$

The proof of Theorem 1 is constructive and provides an estimate of the estimation horizon $N$ required to guarantee RGAS of the MHE proposed in this work. The estimates $\mathcal{N}$ and functions $\Phi$, $\pi_{w}$ and $\pi_{v}$ can be quite conservative, since their derivation involved conservative estimates of noises, errors, stage costs and arrival cost.

Note that the minimum horizon necessary to guarantee RGAS $\mathcal{N}$ depends on $r_{max}$, which depends on the class of disturbances considered (upper bounds of noises and error), the initial value of the prior weighting matrix $P_{0}$ and the bounds of the stage cost. The minimum horizon length is independent of $\left\lVert\boldsymbol{w}\right\rVert$, $\left\lVert\boldsymbol{v}\right\rVert$, and the same $\mathcal{N}$ ensures RGAS for all bounded disturbances and bounded prior error, like the result obtained by muller2017nonlinear ). This implies that we can prove the RGAS property for full information estimator with least–square objective function.

Remark 1**.**

Functions $\phi_{w}$ and $\phi_{v}$ in equations (23) and (24), and hence $\pi_{w}$ and $\pi_{v}$ in equations (35) and (36), do not depend on the estimation horizon $N$ which means that the moving horizon estimator with adaptive arrival cost is RGAS with uniform gains given by (23) and (24).

4 Examples

The following examples will be used to illustrate the results presented in the previous sections and compare the performance of the estimators. The examples considered in this work are taken from muller2017nonlinear for a direct comparison of the results.

4.1 Example 1

The first example considers the system

[TABLE]

The stage cost is chosen as $\ell(w,v)=10w^{2}+10v^{2}$ and the horizon length is $N=10$. The prior weighting is chosen as $\Gamma(\chi)=0.1(\chi-\hat{x}(t|t))^{T}(\chi-\hat{x}(t|t))$ for the MAX estimator (muller2017nonlinear ) and $\Gamma_{t}(\chi)=(\chi-\hat{x}(t|t))^{T}\Pi_{k}^{-1}(\chi-\hat{x}(t|t))$ for the ADAP estimator (our method), where $\Pi_{0}=10I_{2}$ and $\Pi_{k}$ is obtained using equations (5) with $\sigma=0.2$ and $c=1e6$. The MAX estimator uses $\delta=1$, $\delta_{1}=\kappa^{N}$ with $\kappa=0.89^{2}$ and $\delta_{2}=1/N$ (see equation (3) of muller2017nonlinear ). The full information estimator (FIE MAX, see ji2016robust ) is configured with the same parameter used by muller2017nonlinear , maintaining the stage cost and prior weighting $\Gamma_{0}$, and $\delta=1$, $\delta_{1}=\kappa^{t}$ and $\delta_{2}=1/t$.

Table 1 shows the mean square estimation error of each estimator averaged over 300 trials. It can be seen that the proposed estimator average mean square estimation error is smaller than MAX ones and closer to FIE MAX. The main performance difference between ADAP and FIEMAX estimators is the inclusion of the max term in the last one, which allows to follow the sudden changes (see Figures 1 and 2).

Figures 1 and 2 shows simulation results with initial condition $x_{0}=[0.5,0]^{T}$ and prior estimate $\bar{x}_{0}=[0,0]^{T}$. The process and measurement disturbances $w$ and $v$ are sampled from an uniform distribution over the intervals $[-0.3,0.3]$ and $[-0.2,0.2]$, respectively. This figure shows that the estimators that use the max term are able of following the sudden changes, however in the remaining of the signal the MAX estimator is moving away of the FIEMAX while ADAP remains closer.

4.1.1 MHE in the presence of variable measurement noise

Now the MHE estimator is evaluated in the presence of time-varying measurement noise. The variance of the measurement noise is changed from $0.2$ to $1.0$ between times 20 and 40, then it returns to $0.2$.

Table 2 shows the average mean square error in the presence of variable measurement noise. In this case we can see that the behaviour of the proposed estimator is marginally affected by the variations of the measurement, while the mean square error of $x_{0}$ of other estimators increase significantly. These behaviours are due to the adaptation capabilities of the prior weighting updating mechanism, which is able of tracking the changes of noises, in the case of ADAP estimator, and the effect of the max term in MAX and FIEMAX estimators.

Figure 3 shows the evolution of the trace of $P^{-1}_{k-N}$ used in the prior weight of ADAP estimator in both examples. It can be seen that the trace of both matrices grow in similar way, however when the measurement noise changes its variance from $0.2$ to $1.0$ the trace of $P^{-1}_{k-N}$ increases its value (from $12.5$ to $22.5$) and them both traces have the same behaviour.

4.2 Example 2

As a second example, we consider a second order gas-phase irreversible reaction of the form $2A\rightarrow B$. This example has been considered in the context of moving horizon estimation in haseltine2005critical , ji2016robust and muller2017nonlinear . Assuming an isothermal reaction and that the ideal gas law holds, the system dynamics

[TABLE]

where $x=[x_{0},x_{1}]$, $x_{0}$ is the partial pressure of the reactant $A$, $x_{1}$ is the partial pressure of the product $B$, and $k=0.16$ is the reaction rate constant. The measured output of the system is the total pressure. The system is affected with additive process and measurement noise $w$ and $v$ drawn from normal distributions with zero mean and covariance $Q_{w}=0.001^{2}I_{2}$ and $R_{v}=0.1^{2}$, respectively. The stage cost and prior weighting are chosen as $\ell(w,v)=w^{T}Q_{w}^{-1}w+R_{v}^{-1}v^{2}$ and $\Gamma_{t}(\chi)=(\chi-\hat{x}(t|t))^{T}\Pi_{k}^{-1}(\chi-\hat{x}(t|t))$ with $\Pi_{0}=(1/36)I_{2}$, where $\Pi_{k}$ is determined by an extended Kalman filtering recursion in the case of the MAX estimator and the adaptive method in the case of the ADAP estimator with $\sigma=0.1$ and $c=1e6$. For the MAX estimator we use $\delta_{1}=1/N$, $\delta_{2}=1$ and $\delta=0$. In the case of the ADAP estimator, the stage cost weight matrices are chosen as $Q_{w}=0.001I_{2}$ and $R_{v}=0.1$. We use a multiple shooting strategy with a sampling time of $\Delta=0.1$ and we add the restrictions $x_{0}\geq 0$ and $x_{1}\geq 0$.

Table 3 shows the values of the mean squared error computed from the time $10$ (in order to neglect the initial transient error) up to the simulation end time and averaged over 300 trials for horizon sizes of $N=5$ and $N=10$.

Figures 4 and 5 show simulation results with $x_{0}=[3,1]^{T}$ and $\bar{x}_{0}=[0.1,4.5]^{T}$ and horizons of sizes $N=2,5\text{ and }10$, along with results for a full information estimator using the same parameters: the same stage cost $\ell(\cdot)$, prior weighting $\Gamma_{0}$, $\delta=0$, $\delta_{1}=1/t$ and $\delta_{2}=1$. These figures show that the behaviour of ADAP estimator hardly change with horizon length (only the startup behaviours show differences) and no offset in the estimates, while the behaviour of the MAX estimator changes significantly. In addition to the cycling effect caused by the use of the filtered estimate to update $\bar{x}_{k-N}$ [findeisen1997moving, ], the MAX estimator also exhibits offset in the estimate that depends on the estimation horizon length.

5 Conclusions

In this paper we established robust global asymptotic stability for moving horizon estimator with a least-square type cost function for nonlinear detectable (i-IOSS) systems in presence of bounded disturbances. It was also shown that the estimation error converges to zero in case that disturbances converge to zero. This was done for an estimator which uses a least-square type cost function whose arrival cost us updated using adaptive estimation methods. An advantage of this updating mechanism is that the required conditions on prior weighting are such that it can be chosen off-line. Furthermore, it introduces a feedback mechanism between the arrival cost weight and the estimation errors that automatically controls the amount of information used to compute it, which allows to shorten the estimation horizon.

The standard least-square type cost function is typically used in practical applications and RGAS has been proved in muller2017nonlinear . However, for this formulation, the disturbances gains depend on the estimation horizon. Hence, this result does not allow to establish robust global asymptotic stability for a full information estimator. We showed that changing the updating mechanism of arrival cost weight the disturbances gains becomes uniform, allowing to extend the stability analysis to full information estimators with least-square type cost functions.

Acknowledgment

The authors wish to thank the Consejo Nacional de Investigaciones Cientificas y Tecnicas (CONICET) from Argentina, for their support.