1 Introduction↩︎

Robust state estimation for nonlinear systems subject to noise is a problem of high practical relevance. However, if the underlying model is also uncertain and cannot accurately capture the real system behavior, the estimation error may even become unstable if this is not taken into account in the observer design. We address this issue by developing a robust moving horizon estimation (MHE) scheme for simultaneously estimating the unknown states and (constant) parameters of a nonlinear system.

MHE is an optimization-based technique and therefore naturally applicable to nonlinear, potentially constrained systems. Strong robust stability properties of MHE emerge under a mild detectability condition (namely, incremental input/output-to-state stability, i-IOSS), see, e.g., [1]–[4]. These guarantees, however, rely on an exact model of the real system and are therefore not necessarily valid in the case of (parametric) model uncertainties. To address this problem, a \(\min\)-\(\max\) MHE scheme was proposed earlier in [5], where at each time step a least-squares cost function is minimized for the worst case of the model uncertainties. Yet it is often advantageous to not only ensure robustness against model errors, but also to obtain an estimate of the uncertain parameters, since a good model is crucially required for, e.g., (high-performance) control, system monitoring, or fault detection. In this context, an MHE scheme was proposed in [6], treating the parameters as additional states (with constant dynamics), where the temporary loss of observability (due to lack of excitation) is handled by suitable regularization and adaptive weights. However, the robustness properties have not been analyzed, and the imposed conditions for guaranteed state and parameter convergence are not trivial to verify in practice. In [7], MHE under a non-uniform observability condition is considered, which is potentially also suitable to be used for joint state and parameter estimation. The results, however, rely on regularly persistent inputs and, in particular, no fallback strategy is provided in case a lack of excitation occurs in practice during estimation.

An alternative approach to simultaneous state and parameter estimation is provided by adaptive observers, which have been extensively studied in the literature, see, e.g., [8]. Theoretical guarantees usually involve a detectability/observability condition on the system states and a persistence of excitation (PE) condition to establish parameter convergence. Different system classes (usually neglecting disturbances) have been considered, e.g., linear time-varying (LTV) systems [9], Lipschitz nonlinear systems under a linear parameterization [10], nonlinearly parameterized systems [11], or systems in a certain nonlinear adaptive observer canonical form, cf., e.g., [12], [13]. Alternative approaches for systems in canonical forms can be found in, e.g., [14], where more general identifiers are used to estimate the dynamics.

Many results on state and parameter estimation rely on PE conditions that are uniform in time, which is usually restrictive and cannot be guaranteed a priori (except for, e.g., linear systems). To ensure practical applicability, it is essential to investigate weaker (especially non-uniform) excitation conditions. In this context, adaptive observer designs are proposed in, e.g., [15], [16], where boundedness of the state and parameter estimation error is guaranteed without excitation, and exponential stability in the presence of PE. Relaxed excitation conditions have recently received much attention in the context of (pure) parameter estimation of regression models. In [17], it was shown that weaker conditions than PE, however, generally only allow for non-uniform asymptotic stability of LTV systems (that describe the error dynamics of, e.g., simple least squares estimators under linear regression models), which is also consistent with earlier works, e.g., [18]. Using the dynamic regressor extension and mixing idea, exponential convergence could be established for linear regression models (and certain classes of nonlinear ones), merely assuming interval excitation (which is strictly weaker than uniform PE), cf., e.g., [19], [20].

Contribution: We propose an MHE scheme for joint state and parameter estimation for general nonlinear discrete-time systems subject to process disturbances and measurement noise. Our arguments are based on recent MHE results [1], [2], [4], where only state (but no parameter) estimation is considered. The MHE scheme avoids a uniform PE condition and instead uses online information about the current excitation of the parameters to suitably adjust the regularization term in the cost function (Section 3). We establish a bound on the state and parameter estimation error that is valid for all times (even if the parameters are never or only rarely excited), which improves the more often sufficient excitation is present. The bound specializes to a robust global exponential stability property under an additional uniform condition on the maximum duration of insufficient excitation. A method to online monitor the level of excitation that is applicable to general nonlinear systems is provided in Section 4, compare also [21] for further details. In combination, we provide a flexible MHE scheme for robust joint state and parameter estimation with intuitive error bounds that are valid independent of the parameter excitation. The numerical example in Section 5 illustrates that the proposed approach is able to efficiently compensate for phases of weak excitation and ensure reliable estimation results for all times.

In the preliminary conference version [22], we have established exponential convergence of MHE using a joint i-IOSS Lyapunov function for both the states and the parameters, and we have provided a constructive approach for the computation of joint i-IOSS Lyapunov functions that is applicable to a special class of nonlinear systems. In this paper, we show that the existence of such a joint i-IOSS Lyapunov function is equivalent to the detectability of the states and a uniform PE condition of the parameters (which is restrictive). Compared to [22], we consider the much more practically relevant case where the parameters may be insufficiently (or not at all) excited, and derive our results for a significantly larger class of nonlinear systems.

Notation: The set of integers is denoted by \(\mathbb{I}\), the set of all integers greater than or equal to \(a\) for any \(a \in \mathbb{I}\) by \(\mathbb{I}_{\geq a}\), and the set of integers in the interval \([a,b]\) for any \(a,b\in \mathbb{I}\) by \(\mathbb{I}_{[a,b]}\). The floor-function applied to some \(b\in\mathbb{R}_{\geq0}\) is defined as \(\lfloor b \rfloor = \max\{b'\in\mathbb{I}_{\geq0}:b'\leq b \}\). The \(n \times n\) identity matrix is denoted by \(I_n\) and the \(n \times m\) zero matrix by \(0_{n\times m}\), where we omit the indices if the dimension is unambiguous from the context. The weighted Euclidean norm of a vector \(x \in \mathbb{R}^n\) with respect to a positive definite matrix \(Q=Q^\top\) is defined as \(\|x\|_Q=\sqrt{x^\top Q x}\) with \(\|x\|=\|x\|_{I_n}\); the minimal and maximal eigenvalues of \(Q\) are denoted by \(\underline{\lambda}(Q)\) and \(\overline{\lambda}(Q)\), respectively. For two matrices \(A=A^\top\) and \(B=B^\top\), we write \(A \succeq B\) (\(A\succ B\)) if \(A-B\) is positive semi-definite (positive definite). For \(A,B\) positive definite, the maximum generalized eigenvalue (i.e., the largest scalar \(\lambda\) satisfying \(\det (A-\lambda B) = 0\)) is denoted by \(\overline{\lambda}(A,B)\).

2 Problem setup and preliminaries↩︎

We consider discrete-time systems in the form of \[\tag{1} \begin{align} x_{t+1} &= f(x_t,u_t,w_t,p),\tag{2}\\ y_t &= h(x_t,u_t,w_t,p),\tag{3} \end{align}\] where \(t\in\mathbb{I}_{\geq0}\) is the discrete time, \(x_t\in\mathbb{R}^n\) is the state at time \(t\), \(y_t\in\mathbb{R}^p\) is the (noisy) output measurement, \(u_t\in\mathbb{R}^m\) is a known input (e.g., the control input), \(w_t\in\mathbb{R}^q\) is an unknown generalized disturbance input (representing both process disturbances and measurement noise for the sake of conciseness, which covers the standard setting of independent process disturbance and measurement noise as a special case, cf. the simulation example in Section 5), and \(p\in\mathbb{R}^o\) is an unknown but constant parameter. The nonlinear continuous functions \(f:\mathbb{R}^n\times\mathbb{R}^m\times\mathbb{R}^q\times\mathbb{R}^o\rightarrow\mathbb{R}^n\) and \(h:\mathbb{R}^n\times\mathbb{R}^m\times\mathbb{R}^q\times\mathbb{R}^o\rightarrow\mathbb{R}^p\) represent the system dynamics and output equation, respectively. We assume that the unknown true system trajectory satisfies \((x_t,u_t,w_t,p)\in\mathbb{Z}\), \(t\in\mathbb{I}_{\geq 0}\), where \[\mathbb{Z}:=\{(x,u,w,p)\in\mathbb{X}\times\mathbb{U}\times\mathbb{W}\times\mathbb{P}:f(x,u,w,p)\in\mathbb{X}\}\] and \(\mathbb{X}\subseteq \mathbb{R}^n, \mathbb{U}\subseteq \mathbb{R}^m, \mathbb{W}\subseteq \mathbb{R}^q, \mathbb{P}\subseteq \mathbb{R}^o\) are some known closed sets. Such constraints typically arise from the physical nature of the system, e.g., non-negativity of partial pressures, mechanically imposed limits, or parameter ranges. Using this information can significantly improve the estimation results, cf. [23]. If no such sets are known a priori, they can simply be chosen as \(\mathbb{X}= \mathbb{R}^n, \mathbb{U}= \mathbb{R}^m, \mathbb{W}= \mathbb{R}^q, \mathbb{P}= \mathbb{R}^o\).

The overall goal is to compute at each time \(t\in\mathbb{I}_{\geq0}\) the estimates \(\hat{x}_t\) and \(\hat{p}_t\) of the true values \(x_t\) and \(p\) using some a priori estimates \(\hat{x}_0\) and \(\hat{p}_0\) and the past measured input-output sequence \(\{(u_j,y_j)\}_{j=0}^{t-1}\). To this end, we require suitable detectability and excitation properties.

For the stacked input vector \(\bar{w}^\top=[w^\top,p^\top]\), Assumption [ass:IOSS] is equivalent to exponential i-IOSS, which became a standard detectability condition in the context of MHE (for state estimation) in recent years, see, e.g., [1]–[4], [23]. This property implies that the difference between any two state trajectories is bounded by the differences of their initial states, their disturbance inputs, their parameters, and their outputs. Assumption [ass:IOSS] is not restrictive; in fact, by a straightforward extension of the results from [2], [24], it is necessary and sufficient for the existence of robustly stable state estimators if the true parameter is known, and for practically stable state estimators with respect to the parameter error. Moreover, Assumption [ass:IOSS] can be verified using LMIs, cf., e.g, [4].

The set \(\mathbb{E}_T\) contains all trajectory pairs of length \(T\) that exhibit a \(T\)-step distinguishability property with respect to the parameters. Specifically, for two trajectories that share the same initial state and the same disturbance inputs, if they form a pair contained in the set \(\mathbb{E}_T\), it holds that the sum of their output differences is zero if and only if their parameters are the same. In Section 3.3, we discuss the relation between detectability of the states (Assumption [ass:IOSS]), excited trajectory pairs (Definition [def:obs]), uniform excitation, and a uniform joint detectability condition for both the states and the parameters.

3 Moving horizon estimation without uniformly persistently excited data↩︎

3.1 Design↩︎

At each time \(t\in\mathbb{I}_{\geq 0}\), the proposed MHE scheme considers measured past input-output sequences of the system 1 within a moving horizon of length \(N_t = \min\{t, N\}\) for some \(N\in\mathbb{I}_{\geq0}\). The current state and parameter estimates are obtained by solving the following nonlinear program: \[\tag{6} \begin{align}\tag{7} &\min_{\hat{x}_{t-N_t|t},\hat{p}_{|t},\hat{w}_{\cdot|t}}\; J_t(\hat{x}_{t-N_t|t},\hat{p}_{|t},\hat{w}_{\cdot|t},\hat{y}_{\cdot|t}) \\ &\text{s.t. } \hat{x}_{j+1|t}=f(\hat{x}_{j|t},u_j,\hat{w}_{j|t},\hat{p}_{|t}), \;j\in\mathbb{I}_{[t-N_t,t-1]}, \tag{8} \\ & \hat{y}_{j|t}=h_{}(\hat{x}_{j|t},u_j,\hat{w}_{j|t},\hat{p}_{|t}),\;j\in\mathbb{I}_{[t-N_t,t-1]}, \tag{9} \\ &(\hat{x}_{j|t},u_j,\hat{w}_{j|t},\hat{p}_{|t})\in\mathbb{Z},\;j\in\mathbb{I}_{[t-N_t,t-1]} \tag{10} \end{align}\] for all \(t\in\mathbb{I}_{\geq0}\). The decision variables \(\hat{x}_{t-N_t|t}\), \(\hat{p}_{|t}\), and \(\hat{w}_{\cdot|t} = \{\hat{w}_{j|t}\}_{j=t-N_t}^{t-1}\) denote the current estimates of the state at the beginning of the horizon, the parameter, and the disturbance sequence over the horizon, respectively, estimated at time \(t\). Given the past input sequence \(\{u_j\}_{j=t-N_t}^{t-1}\) applied to system 1 , these decision variables (uniquely) define a sequence of state and output estimates \(\{\hat{x}_{j|t}\}_{j=t-N_t}^{t}\) and \(\{\hat{y}_{j|t}\}_{j=t-N_t}^{t-1}\) under 8 and 9 , respectively. We use the cost function \[\begin{align} &J_t(\hat{x}_{t-N_t|t},\hat{p}_{|t},\hat{w}_{\cdot|t},\hat{y}_{\cdot|t}) \nonumber \\ &{=\,} \gamma(N_t)\|\hat{x}_{t-N_t|t}{\,-\,}\bar{x}_{t-N_t}\|_{W_{t-N_t}}^2{+\,} \eta_1^{N_t}\|\hat{p}_{|t}{\,-\,}\bar{p}_{t-N_t}\|_{V_{t-N_t}}^2 \nonumber \\ & \;\;+\sum_{j=1}^{N_t}\eta_2^{j-1}\big(\|\hat{w}_{t-j|t}\|_{Q_{t-j}}^2+\|\hat{y}_{t-j|t}-y_{t-j}\|_{R_{t-j}}^2\big), \label{eq:MHE95objective} \end{align}\tag{11}\] where \(\{y_j\}_{j=t-N_t}^{t-1}\) is the measured output sequence of system 1 . The prior estimates \(\bar{x}_{t-N_t}\) and \(\bar{p}_{t-N_t}\) as well as the cost function parameters \(\gamma(\cdot)\), \(\eta_1\), \(\eta_2\), \(W_t\), \(V_t\), \(Q_t\), and \(R_t\) are defined below. Note that we consider the prediction form of the estimation problem, i.e., without using the most recent measurement \(y_t\) in 11 , which is commonly done to simplify the notation in the theoretical analysis of MHE schemes, cf., e.g., [1]–[4], and see also [23] for a discussion on this topic. We impose the following assumption on the cost function weights.

Assumption [ass:bounds2] ensures large tuning capabilities of the cost function 11 while satisfying certain relations to the detectability and excitation properties from Assumption [ass:IOSS] and Definition [def:obs]. This is conceptually similar to the recent MHE literature (for state estimation) and typically permits a less conservative stability analysis due to the structural similarities between the detectability condition, MHE, and stability, cf., e.g., [4] and compare also [2], [3]. Potential time dependency of the weighting matrices in 14 –17 can be used to incorporate additional knowledge (e.g., by choosing Kalman filter covariance update laws [25], [26]), which can be beneficial in practice to increase estimation performance. Assumption [ass:bounds2] implies that the cost function 11 is radially unbounded in the decision variables, which together with continuity of \(f\) and \(h\) ensures that the estimation problem 6 and 11 admits a globally optimal solution at any time \(t\in\mathbb{I}_{\geq0}\), cf. [23]. We denote a minimizer by \((\hat{x}^*_{t-N_t|t},\hat{p}^*_{|t},\hat{w}_{\cdot|t}^*)\), and the corresponding optimal state sequence by \(\{\hat{x}^*_{j|t}\}_{j=t-N_t}^t\). The resulting state and parameter estimates at time \(t\in\mathbb{I}_{\geq0}\) are then given by \(\hat{x}_t = \hat{x}^*_{t|t}\) and \(\hat{p}_t = \hat{p}^*_{|t}\), yielding the estimation error \[\label{eq:MHE95error} e_t^\top = \big[e_{\mathrm{x},t}^\top,\;e_{\mathrm{p},t}^\top\big] = \big[(\hat{x}_t-x_t)^\top,\;(\hat{p}_t-p)^\top\big].\tag{18}\] It remains to define suitable update laws for the prior estimates to ensure a proper regularization of the cost function 11 . For the state prior, we select \(\bar{x}_{t}=\hat{x}_{t}\) (which is typically called the filtering prior, cf. [23]). For the parameter prior, we propose the following update law: \[\label{eq:pbar} \bar{p}_{t} = \begin{cases} \hat{p}_t,& \text{if t\in\mathbb{I}_{\geq N} and X_{t}\in\mathbb{E}_{N}},\\ \bar{p}_{t-N_t}, & \text{otherwise} \end{cases}\tag{19}\] with \(\bar{p}_0 = \hat{p}_0\), where \(X_{t}:=\Big(\{(\hat{x}_{j|t}^*,u_j,\hat{w}_{j|t}^*,\hat{p}_{|t}^*)\}_{j=t-N_t}^{t-1}, \{({x}_{j},u_j,{w}_{j},p)\}_{j=t-N_t}^{t-1}\Big)\) is the pair of the (unknown) true and the currently optimal trajectory. The update in 19 depends on the currently present level of excitation; in Section 4, we propose a suitable method to practically check whether \(X_{t}\in\mathbb{E}_{N}\) online. In the following, we show how the horizon length \(N\) must be chosen so that the estimation error 18 exhibits a certain robust stability property that is valid regardless of the parameter excitation.

3.2 Stability Analysis↩︎

We consider the two Lyapunov function candidates: \[\begin{align} &\Gamma_1(t,\hat{x},x,\hat{p},p) = U(\hat{x},x) + \|\hat{p}-p\|_{V_t}^2, \tag{20}\\ &\Gamma_2(\hat{x},x,\hat{p},p) = U(\hat{x},x) + c\|\hat{p}-p\|_{\overline{V}}^2, \;c\geq 1, \tag{21} \end{align}\] where \(U\) is from Assumption [ass:IOSS] and \(V_t\preceq\overline{V}\) is from 11 under Assumption [ass:bounds2]. The following two auxiliary results establish fundamental properties of \(\Gamma_1\) and \(\Gamma_2\) for the two cases where the current level of excitation is too low (\(X_t \notin E_N\), Lemma [lem:nonPE]) or sufficiently high (\(X_t \in \mathbb{E}_N\), Lemma [lem:PE]); the proofs can be found in Appendix 7.

Now, let \[\label{eq:rho} \rho := \max\left\{\eta_1^{-N}c_{1}(N)(\eta_\mathrm{x}^N{\,+\,}\gamma(N))\overline{\lambda}(\overline{W},\underline{U}),\eta_2^N\right\}\tag{27}\] and \(c\) be such that \[c = {2}c_{1}(N)/(1-\rho) + 1 \label{eq:c95def}\tag{28}\] with \(c_{1}\) from 23 . The robustness guarantees for the proposed MHE scheme require satisfaction of the following conditions on the horizon length \(N\): \[\begin{align} &2\overline{\lambda}(\overline{W},\underline{U})c_2(c,N)\gamma(N) <1,\tag{29}\\ &c_2(c,N)\eta_1^N<1,\tag{30}\\ &\eta_1^{-N}c_{1}(N)(\eta_\mathrm{x}^N+\gamma(N))\overline{\lambda}(\overline{W},\underline{U}) < 1,\tag{31} \end{align}\] where \(c_1(s)\) and \(\gamma(s)\) are from 23 and 12 , respectively. The conditions 29 and 30 imply that \(\mu\in[0,1)\) in 25 and hence ensure contraction of the state and parameter estimation error in case the excitation condition used in 19 , i.e., \(X_t\in\mathbb{E}_N\), is met. Condition 31 implies \(\rho\in[0,1)\) in 27 and \(c\geq1\) in 28 and hence ensures boundedness of the estimation error in case the excitation condition is not met. Under Assumption [ass:bounds2], there always exists \(N\) sufficiently large such that the contraction conditions 29 –31 are satisfied (to see this, note that the left-hand side of each of these conditions can be bounded by a function that exponentially decays to zero as \(N\rightarrow\infty\), which follows by invoking 12 –13 and uniform boundedness of \(c_1\) and \(c_2\) from 23 and 26 ).

At each time \(t\in\mathbb{I}_{\geq0}\), we split the interval \([0,t]\) into sub-intervals of length \(N\) and the remainder \(l=t-\lfloor{t/N}\rfloor N\): \[\label{eq:times} t = l + \sum_{m=1}^{k} (i_m+1)N + jN,\tag{32}\] where \(k\in\mathbb{I}_{\geq0}\), \(i_m\in\mathbb{I}_{[1,k]}\) for \(k\in\mathbb{I}_{\geq1}\), and \(j\in\mathbb{I}_{\geq0}\) are defined as follows. To this end, let the set \(\mathcal{T}_t := \left\{\tau\in\mathbb{I}_{[N,t]} : t-\tau-N\left\lfloor\frac{t-\tau}{N}\right\rfloor=0, X_{\tau}\in\mathbb{E}_{N}\right\}\) contain all past time instances from the set \(\{t, t-N, t-2N, ... \}\) at which the excitation condition used in 19 was met (in the following also referred to as PE horizons for simplicity). The variable \(k\) denotes the total number of PE horizons that occurred until the current time \(t\) and is defined as the cardinality of \(\mathcal{T}_t\), i.e., \(k:= |\mathcal{T}_t|\). Suppose that \(k\in\mathbb{I}_{\geq1}\). The sequence \(\{t_m\}_{m=1}^k\) contains time instances corresponding to PE horizons, where \(t_1 := \max\{\tau\in\mathcal{T}_t\}\) and \(t_{m+1} = \max\{\tau\in\mathcal{T}_t:\tau<t_m\}\) for \(m\in\mathbb{I}_{[1,k-1]}\) if \(k\in\mathbb{I}_{\geq 2}\). The sequence \(\{i_m\}_{m=1}^k\) denotes the numbers of non-PE horizons (i.e., past time instances where the excitation condition 19 was not met) between two successive times \(t_m\) and \(t_{m+1}\) with \(i_{m} = ({t_m-N-t_{m+1}})/{N}\) for \(m\in\mathbb{I}_{[1,k-1]}\) if \(k\in\mathbb{I}_{\geq 2}\), and² \(i_k = (t_k-N-l)/N\). Finally, \(j\) stands for the number of non-PE horizons that occurred between time \(t\) and \(t_1\) if \(k\geq1\) (and between time \(t\) and \(l\) if \(k=0\)), i.e., \(j := ({t-\max\{\tau\in\mathcal{T}_t,l\}})/{N}\). Overall, the partitioning in 32 allows us to theoretically cover occasional occurrence of PE horizons, which of course includes the special cases in which PE horizons never occur (\(k=0\)), or in which all horizons are PE (\(j=0\), \(i_m=0\) for all \(m=1,...,k\)).

We are now in a position to state our main result. Our key argument is that \(\Gamma_2\) 21 decreases from \(t_m\) to \(t_{m+1}\).

Before proving Theorem [thm:nonPE], we want to highlight some key properties.

1. Theorem [thm:nonPE] provides a bound on the state and parameter estimation error³ that is valid regardless of the parameter excitation; it also applies if the excitation condition in 19 is never met (which corresponds \(k=0\) in 33 ) and constitutes a bounded-disturbance bounded-estimation-error property.

2. Satisfaction of the excitation condition in 19 for some \(t\in\mathbb{I}_{\geq N}\) (leading to non-zero values of \(k\) in 33 ) always improve the error bound 33 with respect to the initial estimation error.

3. If \(k\rightarrow\infty\) for \(t\rightarrow\infty\), then the estimation error converges to a ball centered at the origin with the radius defined by the true disturbances.

If the time between two consecutive PE horizons can be uniformly bounded for all times, Theorem [thm:nonPE] specializes to robust global exponential stability of the estimation error.

3.3 Special case: uniform persistent excitation↩︎

In the following, we consider the special case in which the excitation condition in 19 is always satisfied by the following uniform PE condition.

Assumption [ass:param] essentially imposes that any two system trajectories of certain length form a persistently excited trajectory pair. Robust stability guarantees for joint state and parameter estimation under Assumptions [ass:IOSS] and [ass:param] are provided by Corollary [cor:kappa] (with \(\kappa=0\)); however, we have the following implications, which are proven below.

Proposition [prop:jointIOSS] essentially implies that state detectability (Assumptions [ass:IOSS]) and uniform PE of the parameters (Assumption [ass:param]) is equivalent to uniform detectability (exponential i-IOSS) of the augmented state \(x_\mathrm{a}^\top=[x^\top,p^\top]\). Consequently, under these assumptions one could simply consider the augmented state \(x_\mathrm{a}\) and apply MHE schemes and theory for state estimation (e.g., [1], [2], [4]). However, Assumption [ass:param] is restrictive, usually not satisfied in practice, and its a priori verification is generally impossible. The proposed method from Section 3.1, on the other hand, provides a strict relaxation, since it is applicable in the practically relevant case where the parameters are only rarely (or never) excited (which violates Assumption [ass:param] and hence implies that the augmented state cannot be uniformly detectable (i.e., exponentially i-IOSS) and no joint i-IOSS Lyapunov function satisfying 55 and 56 can exist).

4 Verifying persistent excitation↩︎

The MHE scheme from Section 3 requires detecting whether the PE condition in the update law 19 at a given time \(t\in\mathbb{I}_{\geq N}\) is satisfied or not. To this end, we propose a sufficient condition that can be verified online.

In the following, we use the definitions \(z := (x,u,w,p)\) and \(\tilde{z} := (\tilde{x},u,\tilde{w},\tilde{p})\) for any \((x,u,w,p),(\tilde{x},u,\tilde{w},\tilde{p})\in\mathbb{Z}\). In the remainder of this section, we assume that \(f\) and \(h\) are at least twice continuously differentiable in all of its arguments. Let \(z_s(s) := z + s(\tilde{z}-z)\) for \(s\in[0,1]\) and \[\begin{align} &\textstyle A(z,\tilde{z}) := \int_0^1\frac{\partial f}{\partial x}(z_s(s))ds,\textstyle C(z,\tilde{z}) := \int_0^1\frac{\partial h}{\partial x}(z_s(s))ds,\nonumber \\ &\textstyle B(z,\tilde{z}) := \int_0^1\frac{\partial f}{\partial w}(z_s(s))ds,\textstyle D(z,\tilde{z}) := \int_0^1\frac{\partial h}{\partial w}(z_s(s))ds,\nonumber \\ &\textstyle E(z,\tilde{z}) := \int_0^1\frac{\partial f}{\partial p}(z_s(s))ds,\textstyle F(z,\tilde{z}) := \int_0^1\frac{\partial h}{\partial p}(z_s(s))ds\label{eq:lin} \end{align}\tag{62}\] for all \(z,\tilde{z}\in\mathbb{Z}\). In the following, we require some boundedness properties of the terms in 62 .

Assumption [ass:bounds95MVT] is naturally satisfied for special classes of systems (e.g., with additive \(w\) and \(h\) linear in \(x\), which renders \(B,C,D\) constant) or generally if \(\mathbb{Z}\) is compact.

Consider the trajectory pair \[\left(\{(x_t,u_t,w_t,p)\}_{t=0}^{T-1},\{(\tilde{x}_t,u_t,\tilde{w}_t,\tilde{p})\}_{t=0}^{T-1}\right) \in \mathbb{Z}^T\times\mathbb{Z}^T \label{eq:traj}\tag{65}\] for some \(T\in\mathbb{I}_{\geq0}\), where \(x_{t+1}=f(x_t,u_t,w_t,p)\) and \(\tilde{x}_{t+1}=f(\tilde{x}_t,u_t,\tilde{w}_t,\tilde{p})\) for all \(t\in\mathbb{I}_{[0,T-1]}\). For the sake of brevity, define \[z_t := (x_t,u_t,w_t,p), \;\tilde{z}_t := (\tilde{x}_t,u_t,\tilde{w}_t,\tilde{p})\] for all \(t\in\mathbb{I}_{[0,T-1]}\), and \[\begin{align} Z &:= \left(x_0,p,\{{u}_t\}_{t=0}^{T-1},\{w_t\}_{t=0}^{T-1}\right), \tag{66}\\ \tilde{Z} &:= \left(\tilde{x}_0,\tilde{p},\{{u}_t\}_{t=0}^{T-1},\{\tilde{w}_t\}_{t=0}^{T-1}\right). \tag{67} \end{align}\]

The following result provides a sufficient condition for the trajectory pair 65 to be an element of the set \(\mathbb{E}_T\) satisfying Definition [def:obs] based on matrix recursions.

Verification of 68 requires the knowledge of both trajectories in 65 . This is not the case when applied to the estimation problem presented in Section 3, since the true trajectory is generally unknown. However, we can make local statements based on data from only one of the trajectories that are valid in a surrounding neighborhood. To this end, we define the closed ball centered at some \(Z\in\mathbb{X}\times\mathbb{P}\times\mathbb{U}^T\times\mathbb{W}^T\) of radius \(r>0\) by \(\mathcal{B}(Z,r):=\{\tilde{Z}\in\mathbb{X}\times\mathbb{P}\times\mathbb{U}^T\times\mathbb{W}^T: \|Z-\tilde{Z}\|\leq r\}\). Then, we can evaluate 68 at \((Z,Z)\) and define \[\label{eq:O95T} \mathcal{O}_T(Z) := \mathcal{C}_T(Z,Z).\tag{83}\]

To summarize, for the trajectory pair 65 and \(Z\) and \(\tilde{Z}\) from 66 and 67 , we can make the following conclusion: if \(\tilde{Z}\in\mathcal{B}(Z,r)\) with \(r\) small enough, then \(\mathcal{O}_T(Z)\geq\alpha'>0\Rightarrow \mathcal{C}_T(Z,\tilde{Z})>\alpha>0\) by Proposition [prop:O], which implies that the trajectory pair 65 is an element of the set \(\mathbb{E}_T\) by Proposition [prop:C95PE].

5 Numerical example↩︎

To illustrate our results, we consider the following system \[\begin{align} x_1^+ &= x_1 + t_{\Delta}b_1(x_2-a_1x_1-a_2x_1^2-a_3x_1^3) + w_1,\\ x_2^+ &= x_2 + t_{\Delta}(x_1-x_2+x_3) + w_2,\\ x_3^+ &= x_3-t_{\Delta}b_2x_2 + w_3,\\ y &= x_1 + w_4, \end{align}\] which corresponds to the Euler-discretized modified Chua’s circuit system from [30] using the sampling interval \(t_{\Delta} = 0.01\) under additional disturbances \(w\in\mathbb{R}^4\) and (noisy) output measurements. The parameters are \(b_1=12.8\), \(b_2 = 19.1\), \(a_1=0.6\), \(a_2 = -1.1\), \(a_3 = 0.45\), which leads to a chaotic behavior of the system. In the following, we treat \(w\) as a uniformly distributed random variable with \(|w_i|\leq 10^{-3}, i=1,2,3\) for the process disturbance and \(|w_4|\leq 0.1\) for the measurement noise. We consider the initial condition \(x_0=[1,0,-1]^\top\) and assume that \(x_t\) evolves in the (known) set \(\mathbb{X} = [-1,3]\times[-1,1]\times[-3,3]\). Furthermore, we consider the case where the exact parameter \(a_3=:p\) is unknown but contained in the set \(\mathbb{P}=[0.2,0.8]\). The objective is to compute the state and parameter estimates \(\hat{x}_t\) and \(\hat{p}_t\) by applying the MHE scheme proposed in Section 3 using the initial estimates \(\hat{x}_0=[-1,0.1,2]^\top\) and \(\hat{p}_0=0.2\).

We construct the i-IOSS Lyapunov function \(U\) (Assumption [ass:IOSS]) and the set \(\mathbb{E}_N\) (Definition [def:obs]) using the methods from [4] and Section 4, respectively. The computations are carried out in MATLAB using YALMIP [31] and the semidefinite programming solver MOSEK [32]; details concerning the verification procedure and MHE simulation (including the code) can be found online [33]. We select constant weighting matrices \(W_t = 2P_\mathrm{p}\), \(V_t = 100S_\mathrm{p}\), \(Q_t = 2(Q_\mathrm{x} + Q_\mathrm{p})\), and \(R_t=R_\mathrm{x}+R_\mathrm{p}\) for all \(t\in\mathbb{I}_{\geq0}\), the discount factors \(\eta_1=0.934\), and \(\eta_2=0.9997\), and the horizon length \(N=150\). These choices satisfy Assumption [ass:bounds2] and the contraction conditions 29 –31 (the minimal horizon length for the selected parameters is \(N_{\min}=137\)). While the relatively high value of the horizon length indicates conservatism in our results, it also seems natural here, as the sampling interval \(t_\Delta\) is relatively small compared to the dynamics of the system. In the following, we consider the MHE scheme presented in Section 3 with two different settings: first, without explicit excitation monitoring by naively assuming uniform PE (Assumption [ass:param]); second, using the proposed excitation-dependent adaptive regularization and selecting the parameter estimate \(\hat{p}_t\) in accordance with Remark [rem:est95p]. For the latter, we use the method proposed in Section 4 and check if \(X_t\in\mathbb{E}_{N_t}\) by evaluating \(\mathcal{O}_{N_t}(Z_t^*)\) defined in 83 at the current optimal solution \(Z_t^* = \big(\hat{x}_{t-N_t|t}^*,\hat{p}^*_{|t},\{\hat{w}^*_{j|t}\}_{j=t-N_t}^{t-1}\big)\). If the lower bound \(\alpha_t:=\underline{\lambda}(\mathcal{O}_T(Z_t^*))\) satisfies \(\alpha_t\geq\alpha\) for the predefined threshold \(\alpha=10^{-3}\), we consider that \(X_t\in\mathbb{E}_{N_t}\), and \(X_t\notin\mathbb{E}_{N_t}\) otherwise.

The estimation results are shown Figure 1. Here, we first note that MHE can lead to significant state estimation errors under parametric uncertainties if incorrect (fixed) parameters are used (compare the yellow and purple curves), which highlights the necessity and effectiveness of estimating states and parameters jointly in a combined MHE framework. In particular, both the naive and adaptive (joint) MHE schemes show rapid convergence of the state and parameter estimation errors, see the blue and red curves. Here, the high accuracy of the parameter estimates at the beginning (see the second and third plot) results from a sufficiently high excitation (compare the bottom plot). During the simulation, however, there are also phases of weak excitation—especially in the time interval \([2000,3000]\), where the excitation level \(\alpha_t\) is close to zero. This essentially renders the parameter unobservable, thereby violating Assumption [ass:param]. As a result, the naive MHE scheme provides poor parameter estimates in this interval, which can be clearly seen from the blue curve in the second and third plots in Figure 1. In contrast, MHE with the proposed adaptive regularization explicitly takes the excitation level into account and is thus able to efficiently compensate for phases of weak excitation, resulting in significantly better parameter estimation results (see the red curves).

While parameter unobservability must obviously be addressed appropriately to ensure reliable parameter estimates, it has almost no impact on the state estimates (the blue and red curves coincide in the top plot in Figure 1). In fact, the state trajectories produced by the naive and the adaptive MHE schemes are almost identical, which can be attributed to the fact that the parameter has a vanishing influence on the cost function when it is unobservable. This is also evident in Table ¿tbl:tab:RMSE?, where we provide the root mean square error (RMSE) for the state and parameter estimates. Specifically, the RMSE of the state estimation errors is equal for both MHE schemes (their difference lies in the order of magnitude of \(10^{-15}\)), whereas the RMSE of the parameter estimation error for MHE using adaptive regularization is more than six times smaller than the naive MHE scheme without excitation monitoring.

We now compare the average time \(\tau_{\mathrm{avrg}}\) required to solve the nonlinear program in 6 at each time step \(t\) for MHE with and without adaptive regularization using a standard laptop. Interestingly, the last column in Table ¿tbl:tab:RMSE? shows that the computation times are very similar. This can be attributed to the fact that the additional computation of \(\alpha_t\) required for the adaptation mechanism (which mainly consists of matrix operations) saves time in solving the actual optimization problem due to a better conditioning of the cost function compared to naive MHE, especially due to the regularization update 19 .

6 Conclusion↩︎

In this paper, we have proposed an MHE scheme to estimate the states and the unknown constant parameters of general uncertain nonlinear discrete-time systems. The cost function involves an adaptive regularization term that is adjusted according to the current parameter excitation. We have derived a bound for the state and parameter estimation error that is valid regardless of excitation of the parameter, and in particular also applies if the parameter is never or only rarely excited during operation. Furthermore, we provided a practical method for online excitation monitoring that is applicable to general nonlinear systems. The numerical example has shown that the proposed MHE scheme is able to efficiently compensate for phases of weak excitation and ensures reliable estimation results for all times. An extension of the MHE scheme for time-varying parameters is provided in [34]. Future work might consider component-wise adaptation mechanisms in MHE depending on component-wise excitation monitoring to further improve the overall estimation performance.

7 Proofs of Lemmas [lem:nonPE] and [lem:PE]↩︎

References↩︎