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

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

Nestor N. Deniz, Marina H. Murillo, Guido Sanchez, Lucas M. Genzelis, Leonardo L. Giovanini
###### 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.

## 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 filter is designed minimizing the 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 ), 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 denotes the set of integers in the interval denotes the set of integers greater or equal to . Boldface symbols denote sequences of finite or infinite length, i.e., , respectively. We denote as the finite sequence given at time . By we denote the Euclidean norm of a vector . Let denote the supreme norm of the sequence . A function is of class if is continuous, strictly increasing and . If is also unbounded, it is of class . A function is of class if is non increasing and . A function is of class if is of class for each fixed , and of class for each fixed .

The following inequalities hold for all

 γ(a1+a2+…+an)≤γ(na1)+…+γ(nan),β(a1+a2+…+an,k)≤β(na1,k)+…+β(nan,k). (1)

The preceding inequalities hold since is included in the sequence and functions are non-negative strictly increasing functions.

Bounded sequences: A sequence is bounded if is finite. The set of bounded sequences is denoted as for some

Convergent sequences: A bounded infinite sequence is convergent if as . Let denote the set of convergent sequences :

 Cw\coloneqq{w∈W(wmax)|w is convergent}

Analogously, is defined for the sequence .

### 2.2 Problem statement

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

 xk+1=f(xk,wk),x0=x0yk=h(xk)+vk, (2)

where are the state, process noise, measurement and estimation residuals vectors, respectively. The process disturbance and estimation residuals are unknown but assumed to be bounded, i.e, for some . are compact and convex sets with the null vector belongs to them. In the following we assume that is continuous, locally Lipschitz on and is continuous. The solution to the system (2) at time is denoted by , with initial condition and process disturbance sequence . Furthermore, the initial condition is unknown, but a prior knowledge is assumed to be available and its error is assumed to be bounded, i.e., , .

The solution of the estimation problem aims to find at time an estimate of the current state minimizing a performance metric using by the MHE. At each sampling time , given the previous measurements , the following optimization problem is solved

 min^xk−N|k,^wj|kΨ\coloneqqΓk−N|k(^xk−N|k)+k−1∑j=k−Nℓ(^wj|k,^vj|k)s.t.⎧⎪ ⎪⎨⎪ ⎪⎩^xj+1|k=f(^xj|k,^wj|k),j∈Z[k−N,k−1]yj=h(^xj|k)+^vj|k,^xj|k∈X,^wj|k∈W, ^vj|k∈V, (3)

where is the optimal estimated and is the optimal process noise estimate at sample based on measurements available at time . The process noise and are the optimization variables. The stage cost penalizes the estimated process noise sequence and the estimation residuals , while penalizes the prior estimated . The adequate choice of and , and their parameters, allows to ensure the robust stability of the estimator [muller2017nonlinear]. While the estimation window is not full, , problem (3) can be reformulated and solved as a full information problem

 min^xk−N|k,^wΨ\coloneqqΓ0|k(^x0|k)+k−1∑j=0ℓ(^wj|k,^vj|k)\pars.t.⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩^xj+1|k=f(^xj|k,^wj|k),j∈Z[0,k−1]yj=h(^xj|k)+^vj|k,\par^xj|k∈X, ^wj|k∈W, ^vj|k∈V,

as increases this problem becomes (3) for all .

In previous works, the robust stability of MHE has been achieved by modifying the standard least-square cost function through the inclusion of a –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 with a time–varying prior weight of the form

 Γk−N|k(^xk−N|k)=∥^xk−N|k−¯xk−N∥P−1k−N|k, (4)

whose parameters are recursively updated using the information available at time ([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 is given by

 ϵk−N=yk−N−^yk−N|k,Nk=[1+^xTk−N|k−1Pk−N−1^xk−N|k−1]σ|ϵk−N|22αk=1−1Nk,Wk=⎡⎣I−Pk−N−1^xk−N|k−1^xTk−N|k−11+^xTk−N|k−1Pk−N−1^xk−N|k−1⎤⎦Pk−N−1,Pk−N={1αkWkif 1αkTr(Wk)≤c,Wkotherwise, (5)

where and , where denotes the process noise variance. The prior knowledge of the window is updated using a smoothed estimate ([findeisen1997moving])

 ¯xk−N=^xk−N|k−1. (6)

The optimization problem (3) can be reformulated in terms of the initial condition and the estimated process noises and the residuals along the entire trajectory as follows

 min^x0|k,^wΨ\coloneqqk−1∑j=k−Nℓ(^wj|k,^vj|k)+k−N−1∑j=1αk−N−jkℓ(^wj|k,^vj|k)+αk−NkΓ0|k(^x0|k)\pars.t.⎧⎪ ⎪⎨⎪ ⎪⎩^xj+1|k=f(^xj|k,^wj|k),j∈Z[0,k−1],α∈(0,1]yj=h(^xj|k)+^vj|k,^xj|k∈X, ^wj|k∈W, ^vj|k∈V,

This formulation of problem (3) allows to explicitly see the effect of past data on the current state estimate . In this formulation it is easy to see the exponential averaging of these data. Allowing change in time, the past data has different affects on the current estimates depending on .

Before proceeding to the development of the main results, we state the main properties and assumptions about the prior weighting .

The updating mechanism (5) is a time-varying filter whose inputs are and the initial condition . It generates recursively a real-time estimation of by updating with an exponential time-averaging of . The updating mechanism (5) only use data and it does not rely on a model of the system. The sequence 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 is a continuous function lower bounded by and upper bounded by such that:

 γ––p(|^xk−N|k−¯xk−N|)≤Γk−N(^xk−N|k)Γk−N(^xk−N|k)≤¯γp(|^xk−N|k−¯xk−N|) (7)

for all and

 (8)

where and .

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

 |P−10|ra≤Γk−N(^xk−N|k)≤|P−1∞|ra. (9)
###### Definition 1.

The system (2) is incrementally input/output-to-state stable if there exist functions and such that for every two initial states , , and any two disturbances sequences the following holds for all :

 (10)

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 , since 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 and satisfies the following inequality

 β(r,s)≤cβrps−q (11)

for some , and and .

###### Assumption 3.

The stage cost is a continuous function bounded by such that the following inequalities are satisfied

 γ––w(w)+γ––v(v)≤ℓ(w,v)≤¯γw(w)+¯γv(v). (12)

Functions and from Definition 1 are related with the bounds of stage cost and through the following inequalities

 (13)

for 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 and for , with prior estimate for . The moving horizon state estimator given by equation (3) with adaptive prior weight is robustly globally asymptotically stable (RGAS) if there exists functions and , such that for all , all , the following is satisfied for all

 |xk−^xk|≤Φ(|x0−¯x0|,k)+πw(∥w∥[0,k−1])+πv(∥v∥[0,k−1]). (14)

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., ), then as .

## 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 for nonlinear detectable systems under bounded disturbances. Furthermore, a function exist such that (14) is valid with this and for all estimation horizon .

###### Theorem 1.

Consider an i-IOSS system (2) with disturbances , . Assume that the arrival cost weight matrix of the MHE problem is updated using the adaptive algorithm (5). Moreover, Assumptions 1, 2 and 3 are fulfilled and initial condition is unknown, but a prior estimate is available. Then, the MHE estimator (3) is .

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

 Ψ∗N=Ψ(^x∗k−N|k,^w∗[k−N,k−1])=Γk−N(^x∗k−N|k)+k−1∑j=k−Nℓ(^w∗j|k,^v∗j|k),

which is bounded (Assumptions 1 and 3) and for all by

 Ψ∗N ≤¯¯¯γp(|^x∗k−N|k−¯xk−N|)+N¯¯¯γw(|^w∗j|k|)+N¯γv(|^v∗j|k|), Ψ∗N ≥γ––p(|^x∗k−N|k−¯xk−N|)+Nγ––w(|^w∗j|k|)+Nγ––v(|^v∗j|k|).

Due optimality, the following inequalities hold

 (15)

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

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

 |^xk−N|k−¯xk−N|≤γ––−1p(3¯γp(|xk−N−¯xk−N|))+γ––−1p(3N¯γw(∥w∥))+γ––−1p(3N¯γv(∥v∥)),≤31a|P−10|(|P−1∞|1a|xk−N−¯xk−N|+N1a¯γ1aw(∥w∥)+N1a¯γ1av(∥v∥)).

Analogously, bounds for and can be found

 |^wj|k|≤γ––−1w(3N¯γp(|xk−N−¯xk−N|))+γ––−1w(3¯γw(∥w∥))+γ––−1w(3¯γv(∥v∥)),|^vj|k|≤γ––−1v(3N¯γp(|xk−N−¯xk−N|))+γ––−1v(3¯γw(∥w∥))+γ––−1v(3¯γv(∥v∥)). (16)

Next, let us consider some sample . Assuming that system (2) is i-IOSS with and for all . Since we obtain

 |xk−^xk|k|≤β(|xk−N−^xk−N|k|,N)+γ1(∥∥wj−^wj|k∥∥)+γ2(∥∥vj−^vj|k∥∥). (17)

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

Using Assumptions 1 and 2, function is bounded by

Taking in account that is a symmetric positive definite matrix for all , then , where denotes the maximal eigenvalue of matrix . Denoting as the minimal eigenvalue of matrix and taking in account that , the maximum conditioning number of matrix can be defined as , then can be bounded by

 β(|xk−N−^xk−N|k|,N)≤cβ18p|P−10|(¯γpaw(∥w∥)+¯γpav(∥v∥))+cβNq(2p+CpP−118p)|xk−N−¯xk−N|p. (18)

The first term in the right side of this equation is bounded due the assumption that , while the second term are finite constants. To extend the validness of (18) to the full estimation horizon, an extension of the function at the beginning of the estimation, , is required.

The second term in the right hand side of equation (17), can be bounded by the following inequality

 γ1(∥∥wj−^wj|k∥∥)≤γ1(∥w∥+∥∥^wj|k∥∥)≤γ1(∥w∥+γ––−1w(3N¯γp(|xk−N−¯xk−N|))+γ––−1w(3¯γw(∥w∥))+γ––−1w(3¯γv(∥v∥))).

Recalling Assumption 3, the reader can verify the following inequality

 (19)

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

 (20)

Once an upper bound for the three terms of equation (17) were found, defining and , equation (17) can be rewritten as follows

 (21)

Defining the functions and for all and as follows

 ¯β(r,s)\coloneqqrζsη(CρP−1(cβ18p+λα1min(P−10)(c13α1+c23α2))+cβ2p), (22) (23) (24)

equation (21) can be written as follows

 |xk−^xk|k|≤¯β(|xk−N−¯xk−N|,N)+ϕw(∥w∥)+ϕv(∥v∥). (25)

To guarantee the validity of previous results on the entire time horizon we must extend the definition of . Because of , and for , it is sufficient to define for some to extend the definition of . We would like to determinate the decreasing rate for the function samplings time in the future. In order to do that, let define the constants

 μ∈R>0,δ>2+μ1+μ

and

 rmax\coloneqq δ(1+μ)(ϕw(∥w∥)+ϕv(∥v∥))}

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

 N≥(δζrζ−1maxCρP−1(cβ18p+λα1min(P−10)(c13α1+c23α2)+cβ2p))1η (26)

Adopting an estimator with a window length greater or equal to such that

 ¯β(δr,N)≤(NN)ηr, (27)

the effects of the initial conditions will vanish with a decreasing rate . As , the estimation will entry to the bounded set defined by the noises of the system

 X(w,v)\coloneqq{|xk+j−^xk+j|k+j|≤δ(1+μ)(ϕw(∥w∥)+ϕv(∥v∥))}. (28)

This set define the minimum size region of error space that the error can achieve by removing the effect of errors in initial conditions (). Equation (27) establish a trade off between speed of convergence and window length, which is related with the size of .

For any MHE with adaptive arrival cost and window length two situations can be considered

• The estimator removed the effects of on such that , and

• The estimator has not removed the effects of on such that ,

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

 |xk+N−^xk+N|k+N|≤¯β(|xk−¯xk|,k)+ϕw(∥w∥)+ϕv(∥v∥),≤|xk−¯xk|δ(NN)η+ϕw(∥w∥)+ϕv(∥v∥),≤(2+μ)(ϕw(∥w∥)+ϕv(∥v∥)),≤δ(1+μ)(ϕw(∥w∥)+ϕv(∥v∥)). (29)

This equation implies the fact that the estimation error .

In the other case, when the estimation error is outside of , equations (25) and (27) are recalled again and the following inequalities hold

 |xk+N−^xk+N|k+N|≤|xk−¯xk|δ(NN)η+ϕw(∥w∥)+ϕv(∥v∥),≤|xk−¯xk|δ(NN)η+|xk−¯xk|δ(1+μ)(NN)η,≤|xk−¯xk|(NN)η(2+μδ(1+μ)). (30)

Since , then we have

 θ\coloneqq(NN)η(2+μδ(1+μ))<1. (31)

Equations (30) and (31) reveal a contractive behaviour of the estimation error with as contraction factor. For some finite time the estimation error will decrease until .

In an equivalent formulation, equ