Union and Intersection Operators for Thick Ellipsoid State Enclosures: Application to Bounded-Error Discrete-Time State Observer Design

Abstract

Thick ellipsoids were recently introduced by the authors to represent uncertainty in state variables of dynamic systems, not only in terms of guaranteed outer bounds but also in terms of an inner enclosure that belongs to the true solution set with certainty. Because previous work has focused on the definition and computationally efficient implementation of arithmetic operations and extensions of nonlinear standard functions, where all arguments are replaced by thick ellipsoids, this paper introduces novel operators for specifically evaluating quasi-linear system models with bounded parameters as well as for the union and intersection of thick ellipsoids. These techniques are combined in such a way that a discrete-time state observer can be designed in a predictor-corrector framework. Estimation results are presented for a combined observer-based estimation of state variables as well as disturbance forces and torques in the sense of an unknown input estimator for a hovercraft.

1. Introduction

Predictor-corrector approaches for the model-based estimation of state variables and disturbances were presented in several previous research activities [1,2]. These approaches rely on predicting state enclosures for sets of ordinary differential equations either in an interval-based form [3,4], perform Taylor series expansions with respect to time as well as uncertain parameters and initial conditions [5], or employ Taylor model-based approaches [6,7]. However, these techniques are not directly capable of representing inner and outer bounds of the estimated domain in a unified manner as they were determined by means of affine sets [8], combinations of affine arithmetic with automatic differentiation [9], and data-driven overapproximation techniques for a reachability analysis [10] in current references. On the one hand, the predictor-corrector methods mentioned above provide guaranteed outer bounds of all state variables that are consistent with the state equations as well as with measured output information. Both, state equations and output models, the latter representing sensor characteristics with bounded noise, may include interval parameters to reflect limited knowledge and the influence of model simplifications. On the other hand, however, the lack of information concerning domains that belong to the exact sets of reachable states makes it impossible to distinguish between the phenomena of wide state enclosures as a consequence of uncertainty or the result of wide bounds due to pessimism in numerical evaluations. The latter aspect is well-known in the frame of interval analysis and denoted as overestimation that can be traced back to the so-called dependency effect (numerous variables are treated as independent despite underlying physical or mathematical correlations) or to the wrapping effect (complex-shaped solution sets are replaced conservatively by more simple outer bounds which are subsequently propagated further) [11,12,13,14].

From this description, it becomes obvious that it is desirable to express uncertainty on the bounds of the domains of reachable states in the frame of simulation and state estimation techniques. Recently, corresponding set representations were developed which allow to express either a notion of variables certainly belonging to a solution set or certainly not belonging to the set (denoted as a relative distance measure interval arithmetic [15]); alternatively, the notion of thick intervals, thick boxes, thick functions, and thick ellipsoids can be seen as an approach to handle this kind of uncertainty on set boundaries [16]. Moreover, combinations of fuzzy and interval methods as those presented in [17] represent other currently investigated techniques that aim at expressing uncertainty on state boundaries.

Especially the use of ellipsoids as enclosures for reachable sets [18,19] is a promising approach because such domains naturally arise as descriptions of provable stability domains for nonlinear dynamic systems if the corresponding analysis is performed with the help of quadratic Lyapunov functions. Moreover, propagating ellipsoids over a dynamic system model with sufficient smoothness properties typically leads again to a solution set with boundaries that can effectively be described by enclosing ellipsoids if the domains of interest are sufficiently small. For a proof of this property, see [20], where it has been shown that the Hausdorff distance between the true solution and the ellipsoidal enclosure reduces at least linearly for decreasing domain sizes.

Computations of ellipsoidal enclosures inherently have the property of a constant complexity if they are applied recursively. In addition, their enclosure quality enhances automatically for reducing domain sizes which is in contrast to the propagation of polytopic and zonotopic domains [21,22]. There, the number of vertices of the solution sets—especially after intersecting the results of measurement-based innovation steps with the outcome of state prediction stages—tends to increase in each time step. This leads to the necessity to employ specific techniques to reduce the complexity at least after a certain number of evaluation steps [23]. Although these polytopic and zonotopic techniques help to reduce overestimation, a constant complexity over subsequent time steps for the evaluation of dynamic systems can typically only be observed if the solutions are restricted to domains with a fixed number of vertices such as axis-aligned intervals or parallelepipeds as a special kind of preconditioned interval boxes after a linear change of coordinates.

Previous work of the authors has focused on the development of thick ellipsoid function evaluations [20] and recursive simulation approaches [24], where the inner bounds represent the domains of certainly reachable states and the outer bounds reflect the worst-case uncertainty. So far, only extensions of sufficiently smooth function evaluations (the basic arithmetic operations of addition, subtraction, multiplication and division) as well as the evaluation of nonlinear standard functions (such as trigonometric or exponential functions) were investigated. In order to implement state estimation schemes, it is further necessary to define outer and inner solution enclosures in ellipsoidal form for union and intersection operators. These operators are investigated in more detail in this paper so that they form the basis for the implementation of a discrete-time state estimation procedure. The estimator’s structure resembles the one of a discrete-time Kalman filter [25] or of suitable generalizations for nonlinear systems [26,27,28,29], where a state prediction stage is employed to forecast future reachable states and this a-priori knowledge is subsequently refined in an innovation (i.e., correction) stage as soon as measured data become available.

In Section 2 of this paper, a review of thick ellipsoids as a representation of uncertain state domains for dynamic systems is given. Based on their introduction, inner and outer union as well as intersection operators are derived and summarized in Section 3 which form the basis for state estimation in the frame of (quasi-)linear discrete-time system models. The corresponding novel, thick ellipsoid state estimation procedure is derived in Section 4 before a representative benchmark scenario for state and disturbance estimation of a hovercraft model is discussed in Section 5. Finally, this article is concluded by an outlook on future work in Section 6.

2. Thick Ellipsoids

Thick ellipsoids were introduced by the authors in [20,24] to represent conservative outer bounds of the domains of reachable states of dynamic systems in parallel with inner bounds which represent those states that belong to the reachable solutions with absolute certainty. Geometrically, a thick ellipsoid can be understood as an ellipsoid with an uncertain boundary as illustrated in the left-hand side of Figure 1 that allows for enclosing a complexly shaped domain ${\mathcal{A}}_k$ according to

$${\mathcal{E}}_k^{\mathrm{I}}\subseteq {\mathcal{A}}_k\subseteq {\mathcal{E}}_k^{\mathrm{O}}.$$

(1)

The previous publications [20,24] focused on propagating this ellipsoidal domain via a nonlinear map

$${\mathbf{x}}_{k+1}=\mathbf{f}\left({\mathbf{x}}_k\right)$$

(2)

so that the true outputs ${\mathcal{A}}_{k+1}$ after performing a state prediction according to (2) are again enclosed by a thick ellipsoid in a computationally cheap manner. For that reason, linear matrix inequalities are only employed in offline phases for the derivation of the proposed solution technique. Our procedure avoids the search for the globally optimal (i.e., tightest possible) solution of state prediction (and also the subsequent innovation stages) by a replacement of the online use of linear matrix inequalities or complex minmax optimization tasks by low-dimensional optimizations. Classically, the task of finding ellipsoidal state enclosures would be solved by using the aforementioned, often time-consuming techniques in each time step if the approaches published in [19,30,31,32] were applied directly. However, for the goal of a reduction of the computational effort and, hence, for a simplification of real-time implementability, we propose low-dimensional optimization procedures which avoid the search for the typically large number of degrees of freedom if the optimal inner and outer ellipsoids ${\mathcal{E}}_k^{\mathrm{I}}$ and ${\mathcal{E}}_k^{\mathrm{O}}$, respectively, were computed.

In the following, we recapitulate the mathematical definitions of thick ellipsoids together with binary operators and function extensions. Note, similar concepts were also introduced in [16] for so-called thick intervals (which represent scalar interval variables with uncertain lower and upper bounds) as well as thick boxes as the Cartesian product of thick intervals. A closely related concept is the so-called type-2 interval arithmetic discussed in detail in [15].

Definition 1

(Thick ellipsoid). Define a thick ellipsoid $\left(\left(\mathcal{E}\right)\right)=\left(\left(\mathcal{E}\right)\right)\left(\mathit{\mu },\mathbf{\Gamma },\left[\underset{¯}{\rho };\overline{\rho }\right]\right)$ as a subset of the power set $\mathcal{P}\left({\mathbb{R}}^n\right)$ so that

$$\left(\left(\mathcal{E}\right)\right)=\left\{\mathcal{A}\in \mathcal{P}\left({\mathbb{R}}^n\right)|{\mathcal{E}}^{\mathrm{I}}\subseteq \mathcal{A}\subseteq {\mathcal{E}}^{\mathrm{O}}\right\}$$

(3)

with

$$\begin{array}{cc}\hfill {\mathcal{E}}^{\mathrm{I}}& =\left\{\mathbf{x}\in {\mathbb{R}}^n|{\left(\mathbf{x}-\mathit{\mu }\right)}^T{\left(\underset{¯}{\rho }\mathbf{\Gamma }\right)}^{-T}{\left(\underset{¯}{\rho }\mathbf{\Gamma }\right)}^{-1}\left(\mathbf{x}-\mathit{\mu }\right)\le 1\right\},\hfill \\ \hfill {\mathcal{E}}^{\mathrm{O}}& =\left\{\mathbf{x}\in {\mathbb{R}}^n|{\left(\mathbf{x}-\mathit{\mu }\right)}^T{\left(\overline{\rho }\mathbf{\Gamma }\right)}^{-T}{\left(\overline{\rho }\mathbf{\Gamma }\right)}^{-1}\left(\mathbf{x}-\mathit{\mu }\right)\le 1\right\}\hfill \end{array}$$

(4)

and $0\le \underset{¯}{\rho }\le \overline{\rho }$.

Definition 2.

(Thick ellipsoid binary operators and function extensions). A thick ellipsoid extension of the binary operators $\diamond \in \{+,-,·,/,\cup ,\cap \}$ (For the case of division, the value zero is assumed not to belong to the denominator expression in analogy to classical interval arithmetic [11]) satisfies the relation

$$\left\{\begin{array}{c}\mathcal{A}\in \left(\left(\mathcal{A}\right)\right)\hfill \\ \mathcal{B}\in \left(\left(\mathcal{B}\right)\right)\hfill \\ \mathcal{C}=\mathcal{A}\diamond \mathcal{B}\hfill \end{array}\right.⟹\mathcal{C}\in \left(\left(\mathcal{A}\right)\right)\diamond \left(\left(\mathcal{B}\right)\right).$$

(5)

The quantity $\left(\left(\mathcal{C}\right)\right)=\left(\left(\mathcal{A}\right)\right)\diamond \left(\left(\mathcal{B}\right)\right)$ is also a thick ellipsoid, which is typically neither minimal with respect to its width nor uniquely defined. Analogously, $\left(\left(\mathbf{f}\right)\right)$ is a thick ellipsoid function extension of $\mathbf{f}:{\mathbb{R}}^n\mapsto {\mathbb{R}}^m$, if

$$\left\{\begin{array}{c}\mathcal{A}\in \left(\left(\mathcal{A}\right)\right)\hfill \\ \mathcal{B}=\mathbf{f}\left(\mathcal{A}\right)\hfill \end{array}\right.⟹\mathcal{B}\in \left(\left(\mathcal{B}\right)\right)=\left(\left(\mathbf{f}\right)\right)\left(\left(\mathcal{A}\right)\right).$$

(6)

For the implementation of algorithms that allow for evaluating nonlinear mapping functions (2) with the help of thick ellipsoids, the reader is referred to [20,24]. In this paper, we impose a specific, (quasi-)linear structure of the considered system models, so that the thick ellipsoid function evaluation can be simplified by a newly derived simulation scheme according to Section 3.1.

3. Union and Intersection of Thick Ellipsoid Enclosures

Parameter-dependent quasi-linear system models in the form

$${\mathbf{x}}_{k+1}=\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{x}}_k$$

(7)

with ${\mathbf{p}}_k\in \left[{\underset{¯}{\mathbf{p}}}_k;{\overline{\mathbf{p}}}_k\right]$, where ${\underset{¯}{p}}_{j,k}\le p_{j,k}\le {\overline{p}}_{j,k}$, $j\in \{1,\dots ,m\}$, are a special case of the nonlinear system representation introduced in (2).

In many cases, they arise if discrete-time state equations (2) with an equilibrium state at ${\mathbf{x}}_{k+1}={\mathbf{x}}_k=\mathbf{0}$ are re-written by means of symbolic formula manipulation in which the dependence on the state vector ${\mathbf{x}}_k$ is factored out from the system model (2). In such cases, it needs to be ensured that all entries of $\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)$ are well-defined because they do not contain any singularities after factoring out the state vector ${\mathbf{x}}_k$ and are, therefore, finite for the complete operating domain of interest. In addition, these models arise if centered-form representations of nonlinear system models are computed either by determining interval extensions of the system’s Jacobian or by means of techniques for slope arithmetic. For details concerning these latter two options, the reader is referred to [33,34,35].

3.1. Mapping of Thick Ellipsoidal Domains via (Quasi-)Linear System Models

In general, the evaluation of system models in the form (7) is possible with the help of the general solution approach presented in [20,24] if an augmented state vector is introduced that contains not only the actual state variables ${\mathbf{x}}_k$ but also the parameter vector ${\mathbf{p}}_k$ in terms of discretized integrator disturbance models for which ${\mathbf{p}}_{k+1}={\mathbf{p}}_k$ holds. However, the examples investigated in [24] have shown that doing so leads to conservativeness of the solutions. The reason for this observation is that the parameters ${\mathbf{p}}_k$ become correlated with the state vector ${\mathbf{x}}_k$ so that their outer bounds start to expand despite the aforementioned assumption of uncertain but temporally constant values. Hence, the methods summarized in this subsection (together with the intersection operators included in the state estimator according to Section 4) resolve the corresponding issue by introducing uncertain parameters as entries in the system matrix. When doing so, their enclosures may even be reset at each sampling instant k to new interval bounds. As shown subsequently, this option allows for analogously dealing with (outer) state bounds influencing the dynamics matrix $\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)$.

For the sake of compactness, the short-hand notation ${\mathbf{A}}_k=\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)\in \left[{\mathbf{A}}_k\right]$ is used in the following to denote an interval matrix containing the worst-case outer range bounds of all matrix entries at a specific time instant k.

3.1.1. Outer Bounds

Given an ellipsoid

$${\mathcal{E}}_k=\left\{{\mathbf{x}}_k\in {\mathbb{R}}^n|{\mathbf{x}}_k^T·{\mathbf{Q}}_k^{-1}·{\mathbf{x}}_k\le 1\right\}$$

(8)

that is described by a positive definite shape matrix ${\mathbf{Q}}_k={\mathbf{Q}}_k^T>0$ and centered at the origin ${\mathbf{x}}_k=\mathbf{0}$ of the state space, it is desired to enclose the exact (generally non-ellipsoidal) solution set of the mapping (7) in terms of a guaranteed outer ellipsoidal bound

$${\mathcal{E}}_{k+1}^{\mathrm{O}}=\left\{{\mathbf{x}}_{k+1}\in {\mathbb{R}}^n|\frac{1}{{\alpha }_{\mathrm{O},k+1}^2}·{\mathbf{x}}_{k+1}^T·{\mathbf{Q}}_{k+1}^{-1}·{\mathbf{x}}_{k+1}\le 1\right\}$$

(9)

with degrees of freedom concerning the choice of the shape matrix ${\mathbf{Q}}_{k+1}$ and its stretch parameter ${\alpha }_{\mathrm{O},k+1}>0$. The latter will be computed by means of a simple line-search optimization procedure to avoid the necessity for an online solution of linear matrix inequalities.

With the help of the free matrix variable ${\mathcal{R}}_k$, the shape matrix of the predicted ellipsoid is described by

$${\mathbf{Q}}_{k+1}={\mathbf{Q}}_{k+1}^T={\tilde{\mathbf{A}}}_k·{\mathcal{R}}_k·{\tilde{\mathbf{A}}}_k^T>0,{\mathcal{R}}_k={\mathcal{R}}_k^T>0,$$

(10)

where ${\tilde{\mathbf{A}}}_k$ is some invertible point matrix typically satisfying

$${\tilde{\mathbf{A}}}_k\in \left[{\mathbf{A}}_k\right].$$

(11)

This definition is inherited from the propagation of a representative point in the interior of ${\mathcal{E}}_k$ which was chosen to be its midpoint in [20,24].

Remark 1.

Possible choices for the parameterization of the matrices ${\mathcal{R}}_k$ and ${\tilde{\mathbf{A}}}_k$ introduced in (10) and (11) are discussed in the illustrating example in Section 3.1.3.

Assuming regularity of $\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{Q}}_k·{\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)$ for all ${\mathbf{x}}_k\in {\mathcal{E}}_k$, ${\mathbf{p}}_k\in \left[{\mathbf{p}}_k\right]$, the exact solution of the mapping (7) is defined by

$$\begin{array}{cc}\hfill {\mathcal{E}}_{k+1}^{*}=& \{{\mathbf{x}}_{k+1}\in {\mathbb{R}}^n|{\mathbf{x}}_{k+1}^T·{\left({\mathbf{A}}_k·{\mathbf{Q}}_k·{\mathbf{A}}_k^T\right)}^{-1}·{\mathbf{x}}_{k+1}\le 1,\hfill \\ \hfill & \mathrm{where}{\mathbf{A}}_k=\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)\mathrm{with}{\mathbf{x}}_k\in {\mathcal{E}}_k,{\mathbf{p}}_k\in \left[{\mathbf{p}}_k\right]\},\hfill \end{array}$$

(12)

where ${\mathcal{E}}_{k+1}^{*}\subseteq {\mathcal{E}}_{k+1}^{\mathrm{O}}$ needs to be satisfied for (9) to represent a guaranteed outer bound of the solution domain. Note, even if this regularity assumption does not hold, the following Theorem 1 remains true.

Theorem 1

(Guaranteed outer bound). The ellipsoid ${\mathcal{E}}_{k+1}^{\mathrm{O}}$ according to Equation (9) is a guaranteed outer bound for ${\mathcal{E}}_{k+1}^{*}$ if the matrix inequality

$${\mathcal{M}}_{k+1}:=\left[\begin{array}{cc}-{\mathbf{Q}}_k^{-1}& {\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\tilde{\mathbf{A}}}_k^{-T}\\ {\tilde{\mathbf{A}}}_k^{-1}·\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)& -{\alpha }_{\mathrm{O},k+1}^2{\mathcal{R}}_k\end{array}\right]\le 0$$

(13)

holds for all ${\mathbf{x}}_k\in {\mathcal{E}}_k$ and ${\mathbf{p}}_k\in \left[{\mathbf{p}}_k\right]$.

Proof. 

Without any assumptions on the regularity of $\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{Q}}_k·{\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)$, the relation ${\mathcal{E}}_{k+1}^{*}\subseteq {\mathcal{E}}_{k+1}^{\mathrm{O}}$ is satisfied, if mapping the domain ${\mathcal{E}}_{k+1}^{\mathrm{O}}$ one time step back leads to a domain ${\mathcal{S}}_k$ that is a superset of the ellipsoid ${\mathcal{E}}_k$ according to ${\mathcal{S}}_k\supseteq {\mathcal{E}}_k$ with

$$\begin{array}{cc}\hfill {\mathcal{S}}_k=& \{{\mathbf{x}}_k\in {\mathbb{R}}^n|{\mathbf{x}}_k^T·{\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·\frac{1}{{\alpha }_{\mathrm{O},k+1}^2}·{\left({\tilde{\mathbf{A}}}_k·{\mathcal{R}}_k·{\tilde{\mathbf{A}}}_k^T\right)}^{-1}·\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{x}}_k\le 1,\hfill \\ \hfill & \mathrm{where}{\mathbf{x}}_k\in {\mathcal{E}}_k,{\mathbf{p}}_k\in \left[{\mathbf{p}}_k\right]\}.\hfill \end{array}$$

(14)

Geometrically speaking, this means that the level curve of height 1 of the quadratic form in (8), which corresponds to the bound of the ellipsoid ${\mathcal{E}}_k$, must be an inner bound for the level curve of identical height of ${\mathcal{S}}_k$. This condition corresponds to the scalar inequality

$${\mathbf{x}}_k^T·{\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·\frac{1}{{\alpha }_{\mathrm{O},k+1}^2}·{\left({\tilde{\mathbf{A}}}_k·{\mathcal{R}}_k·{\tilde{\mathbf{A}}}_k^T\right)}^{-1}·\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{x}}_k\le {\mathbf{x}}_k^T·{\mathbf{Q}}_k^{-1}·{\mathbf{x}}_k$$

(15)

for all ${\mathbf{x}}_k\in {\mathcal{E}}_k$ and ${\mathbf{p}}_k\in \left[{\mathbf{p}}_k\right]$.

Factoring out the vector ${\mathbf{x}}_k$ to the left and right of (15) turns this scalar inequality into the matrix inequality

$${\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·\frac{1}{{\alpha }_{\mathrm{O},k+1}^2}·{\left({\tilde{\mathbf{A}}}_k·{\mathcal{R}}_k·{\tilde{\mathbf{A}}}_k^T\right)}^{-1}·\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)-{\mathbf{Q}}_k^{-1}\le 0,$$

(16)

which (for the choice of a regular matrix ${\tilde{\mathbf{A}}}_k$) is equivalent to the condition

$${\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\tilde{\mathbf{A}}}_k^{-T}·{\left({\alpha }_{\mathrm{O},k+1}^2{\mathcal{R}}_k\right)}^{-1}·{\tilde{\mathbf{A}}}_k^{-1}·\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)-{\mathbf{Q}}_k^{-1}\le 0.$$

(17)

Removing the difference of the two matrices in (17) with the help of the Schur complement formula [36,37] leads to Equation (13) which completes the proof of Theorem 1. □

Remark 2.

As stated before, the quasi-linear system matrix $\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)$ is bounded by an element-wise defined interval matrix $\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)\in \left[{\mathbf{A}}_k\right]=\left[{\underset{¯}{\mathbf{A}}}_k;{\overline{\mathbf{A}}}_k\right]$. Therefore, a conservative solution to the matrix inequality in Theorem 1 can be determined by searching for a state- and parameter-independent, strictly positive definite matrix ${\alpha }_{\mathrm{O},k+1}^2{\mathcal{R}}_k$ that satisfies the interval-valued matrix inequality

$${\mathcal{M}}_{k+1}\in \left[{\mathcal{M}}_{k+1}\right]=\left[\begin{array}{cc}-{\mathbf{Q}}_k^{-1}& {\left[{\mathbf{A}}_k\right]}^T·{\tilde{\mathbf{A}}}_k^{-T}\\ {\tilde{\mathbf{A}}}_k^{-1}·\left[{\mathbf{A}}_k\right]& -{\alpha }_{\mathrm{O},k+1}^2{\mathcal{R}}_k\end{array}\right]\le 0.$$

(18)

This matrix inequality can be cast into a linear one by the substitution

$${\mathcal{P}}_{k+1}^{\mathrm{O}}={\left({\mathcal{P}}_{k+1}^{\mathrm{O}}\right)}^T={\alpha }_{\mathrm{O},k+1}^2{\mathcal{R}}_k>0$$

(19)

to make standard routines for linear matrix inequalities such as theMatlabtoolboxSeDuMiin combination withYalmip[38,39] applicable if (18) is further replaced by the convex combination of a finite number of extremal realizations with the help of the ideas given in [40,41]. The extremal realizations can be obtained by computing an interval extension of the state- and parameter-dependent system matrix $\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)$ with the help of toolboxes such asIntLab[42] forMatlaband determining all matrix entries that have a non-zero diameter. Then, it is possible to represent (18) by the convex polytope

$$\mathcal{D}=\left\{{\mathcal{M}}_{k+1}|{\mathcal{M}}_{k+1}\left(\xi \right)=\sum _{v=1}^{n_v}{\xi }_v·{\mathcal{M}}_{v,k+1};\sum _{v=1}^{n_v}{\xi }_v=1;{\xi }_v\ge 0\right\},$$

(20)

which is only used in the following example for the sake of comparison with the proposed method, where for the list of all $n_v=2^{n^2}$ interval vertices formed of ${\underset{¯}{\mathbf{A}}}_k$ and ${\overline{\mathbf{A}}}_k$ the corresponding extremal realizations

$${\mathcal{M}}_{v,k+1}=\left[\begin{array}{cc}-{\mathbf{Q}}_k^{-1}& {\mathbf{A}}_{v,k}^T·{\tilde{\mathbf{A}}}_k^{-T}\\ {\tilde{\mathbf{A}}}_k^{-1}·{\mathbf{A}}_{v,k}& -{\mathcal{P}}_{k+1}^{\mathrm{O}}\end{array}\right],v\in \{1,\dots ,n_v\},$$

(21)

are defined and need to be negative semi-definite with ${\mathcal{M}}_{v,k+1}\le 0$. Note, this formulation of all (21) directly accounts for the structural symmetry of (18). A (nearly) optimal solution can be obtained by minimizing $\mathrm{trace}\left\{{\mathcal{P}}_{k+1}^{\mathrm{O}}\right\}$, cf. [43,44], which can be employed as a substitute for the optimization task involving the logarithm of the shape matrix determinant according to ([45], Appendix C).

To avoid the necessity to solve linear matrix inequalities in each temporal discretization step of a dynamic system, the following Theorem 2 describes a relaxed version of the solution procedure after a-priori fixing the matrix ${\mathcal{R}}_k$, for example, by setting ${\mathcal{R}}_k={\mathbf{Q}}_k$.

Theorem 2

(Eigenvalue criterion for outer ellipsoidal bounds). For a fixed combination of the shape and stretch parameters ${\mathcal{R}}_k$ and ${\alpha }_{\mathrm{O},k+1}$, respectively, ${\mathcal{E}}_{k+1}^{\mathrm{O}}$ is a guaranteed outer enclosure of the solution set ${\mathcal{E}}_{k+1}^{*}$, if

$${\lambda }_i\left(\mathrm{mid}\left(\left[{\mathcal{M}}_{k+1}\right]\right)\right)+\rho \left(\mathrm{rad}\left(\left[{\mathcal{M}}_{k+1}\right]\right)\right)\le 0$$

(22)

holds for all $i\in \{1,\dots ,2n\}$, where ${\lambda }_i\left(\mathrm{mid}\left(\left[{\mathcal{M}}_{k+1}\right]\right)\right)$ is the i-th eigenvalue of $\frac{1}{2}·\left({\underset{¯}{\mathcal{M}}}_{k+1}+{\overline{\mathcal{M}}}_{k+1}\right)$ and $\rho \left(\mathrm{rad}\left(\left[{\mathcal{M}}_{k+1}\right]\right)\right)$ the spectral radius of $\frac{1}{2}·\left({\overline{\mathcal{M}}}_{k+1}-{\underset{¯}{\mathcal{M}}}_{k+1}\right)$.

Proof. 

Due to the structural symmetry of $\left[{\mathcal{M}}_{k+1}\right]$, all possible eigenvalues ${\lambda }_i\left({\mathcal{M}}_{k+1}\right)$ of ${\mathcal{M}}_{k+1}\in \left[{\mathcal{M}}_{k+1}\right]$ are guaranteed to be real-valued. According to Rohn [46,47], interval bounds for the range of each eigenvalue are given by the enclosures

$${\lambda }_i\left({\mathcal{M}}_{k+1}\right)\in \left[{\lambda }_{i,\mathcal{M}}\right]={\lambda }_i\left(\mathrm{mid}\left(\left[{\mathcal{M}}_{k+1}\right]\right)\right)+\left[-1;1\right]·\rho \left(\mathrm{rad}\left(\left[{\mathcal{M}}_{k+1}\right]\right)\right)$$

(23)

located on the real axis of the complex plane. If $sup\left(\left[{\lambda }_{i,\mathcal{M}}\right]\right)\le 0$ holds for all $i\in \{1,\dots ,2n\}$, Theorem 1 is satisfied due to the negative semi-definiteness of the interval matrix (18). □

Remark 3.

To enhance the efficiency of the enclosure technique, the eigenvalue bounds from (23) can be combined with the upper eigenvalue bound ${\lambda }_{max,\mathcal{M}}$ obtained from the Gershgorin circle criterion [48]

$${\lambda }_{max,\mathcal{M}}\le \underset{i\in \{1,\dots ,2n\}}{max}\left\{sup\left\{\left[{\mathcal{M}}_{i,i}\right]+\sum _{j=1,j\ne i}^{2n}\left|\left[{\mathcal{M}}_{i,j}\right]\right|\right\}\right\}.$$

(24)

If this is not efficient, all vertices of the interval matrix (18) could be tested for negative semi-definiteness.

Remark 4.

The optimal, i.e., smallest volume ellipsoid after fixing the matrix ${\mathcal{R}}_k$, for example, by the choice ${\mathcal{R}}_k={\mathbf{Q}}_k$, is obtained with the help of Theorem 2 if the parameter ${\alpha }_{\mathrm{O},k+1}>0$ is minimized by a line search procedure so that (22) remains satisfied. For ${\alpha }_{\mathrm{O},k+1}\to \infty$, the outer bound becomes infinitely large.

3.1.2. Inner Bounds

Inner ellipsoidal bounds

$${\mathcal{E}}_{k+1}^{\mathrm{I}}=\left\{{\mathbf{x}}_{k+1}\in {\mathbb{R}}^n|\frac{1}{{\alpha }_{\mathrm{I},k+1}^2}·{\mathbf{x}}_{k+1}^T·{\mathbf{Q}}_{k+1}^{-1}·{\mathbf{x}}_{k+1}\le 1\right\}$$

(25)

with the shape matrix defined in (10) describe domains that are guaranteed to belong to the solution set according to ${\mathcal{E}}_{k+1}^{\mathrm{I}}\subseteq {\mathcal{E}}_{k+1}^{*}$ are characterized by the following Theorem 3.

Theorem 3

(Guaranteed inner bound). The ellipsoid ${\mathcal{E}}_{k+1}^{\mathrm{I}}$ according to Equation (25) is a guaranteed outer bound for ${\mathcal{E}}_{k+1}^{*}$ if the matrix inequality

$${\mathcal{N}}_{k+1}:=\left[\begin{array}{cc}{\alpha }_{\mathrm{I},k+1}^{-2}·{\mathcal{R}}_k^{-1}& {\left({\tilde{\mathbf{A}}}_k^{-1}·\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)\right)}^{-T}\\ {\left({\tilde{\mathbf{A}}}_k^{-1}·\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)\right)}^{-1}& {\mathbf{Q}}_k\end{array}\right]\ge 0$$

(26)

holds for all ${\mathbf{x}}_k\in {\mathcal{E}}_k$ and ${\mathbf{p}}_k\in \left[{\mathbf{p}}_k\right]$.

Proof. 

The proof of Theorem 3 follows a similar reasoning as the proof of Theorem 1. If the matrix $\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{Q}}_k·{\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)$ has full rank for all ${\mathbf{x}}_k\in {\mathcal{E}}_k$ and ${\mathbf{p}}_k\in \left[{\mathbf{p}}_k\right]$, ${\mathcal{E}}_{k+1}^{\mathrm{I}}\subseteq {\mathcal{E}}_{k+1}^{*}$ holds, if the level curve of height 1 of the quadratic form in (25), which corresponds to the bound of the ellipsoid ${\mathcal{E}}_{k+1}^{\mathrm{I}}$, is fully contained in the interior of ${\mathcal{E}}_{k+1}^{*}$ according to the scalar inequality

$${\mathbf{x}}_{k+1}^T·\frac{1}{{\alpha }_{\mathrm{I},k+1}^2}·{\left({\tilde{\mathbf{A}}}_k·{\mathcal{R}}_k·{\tilde{\mathbf{A}}}_k^T\right)}^{-1}·{\mathbf{x}}_{k+1}\ge {\mathbf{x}}_{k+1}^T·{\left(\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{Q}}_k·{\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)\right)}^{-1}·{\mathbf{x}}_{k+1}$$

(27)

This condition corresponds to the matrix inequality

$$\frac{1}{{\alpha }_{\mathrm{I},k+1}^2}·{\left({\tilde{\mathbf{A}}}_k·{\mathcal{R}}_k·{\tilde{\mathbf{A}}}_k^T\right)}^{-1}-{\left(\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{Q}}_k·{\mathbf{A}}^T\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)\right)}^{-1}\ge 0,$$

(28)

which is equivalent to

$$\frac{1}{{\alpha }_{\mathrm{I},k+1}^2}·{\mathcal{R}}_k^{-1}-{\tilde{\mathbf{A}}}_k^T·{\mathbf{A}}^{-T}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\mathbf{Q}}_k^{-1}·{\mathbf{A}}^{-1}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\tilde{\mathbf{A}}}_k\ge 0$$

(29)

under the above-mentioned regularity assumption.

Removing the difference of matrices in (29) by the application of the Schur complement formula completes the proof of Theorem 3. □

Remark 5.

In analogy to Equation (18), a conservative formulation of Theorem 3 can be given by the interval-valued matrix inequality

$${\mathcal{N}}_{k+1}\in \left[{\mathcal{N}}_{k+1}\right]=\left[\begin{array}{cc}{\alpha }_{\mathrm{I},k+1}^{-2}·{\mathcal{R}}_k^{-1}& {\left({\tilde{\mathbf{A}}}_k^{-1}·\left[{\mathbf{A}}_k\right]\right)}^{-T}\\ {\left({\tilde{\mathbf{A}}}_k^{-1}·\left[{\mathbf{A}}_k\right]\right)}^{-1}& {\mathbf{Q}}_k\end{array}\right]\ge 0$$

(30)

for which a common solution after the linearizing variable substitution

$${\mathcal{P}}_{k+1}^{\mathrm{I}}={\left({\mathcal{P}}_{k+1}^{\mathrm{I}}\right)}^T={\alpha }_{\mathrm{I},k+1}^{-2}·{\mathcal{R}}_k^{-1}>0$$

(31)

needs to be found. An optimization of this solution is possible by a minimization of $\mathrm{trace}\left\{{\mathcal{P}}_{k+1}^{\mathrm{I}}\right\}$, (replacing the optimization task involving the logarithm of the shape matrix determinant according to [45] Appendix C), where for a polytopic representation of the domain $\left[{\mathcal{N}}_{k+1}\right]$ according to (21) the entry-wise defined vertices of the matrix inverse ${\left({\tilde{\mathbf{A}}}_k^{-1}·\left[{\mathbf{A}}_k\right]\right)}^{-1}$ need to be found. A relaxation into a computationally less expensive eigenvalue test is given by Theorem 4.

Theorem 4

(Eigenvalue criterion for inner ellipsoidal bounds). For a fixed combination of the shape and stretch parameters ${\mathcal{R}}_k$ and ${\alpha }_{\mathrm{I},k+1}$, respectively, ${\mathcal{E}}_{k+1}^{\mathrm{I}}$ is a guaranteed inner enclosure of the solution set ${\mathcal{E}}_{k+1}^{*}$, if

$${\lambda }_i\left(\mathrm{mid}\left(\left[{\mathcal{N}}_{\mathbf{k}+\mathbf{1}}\right]\right)\right)+\rho \left(\mathrm{rad}\left(\left[{\mathcal{N}}_{\mathbf{k}+\mathbf{1}}\right]\right)\right)\ge 0$$

(32)

holds for all $i\in \{1,\dots ,2n\}$, where the $\mathrm{mid}$ and $\mathrm{rad}$ operators are defined as in Theorem 2.

Proof. 

The proof of Theorem 4 results from a direct application of the proof of Theorem 2 after exchanging the sign condition of the respective inequalities, where $inf\left(\left[{\lambda }_{i,\mathcal{N}}\right]\right)\ge 0$ needs to hold for all $i\in \{1,\dots ,2n\}$. □

Remark 6.

To enhance the efficiency of the enclosure technique, the eigenvalue bounds from (32) can be combined with the lower eigenvalue bound ${\lambda }_{min,\mathcal{N}}$ obtained from the Gershgorin circle criterion

$${\lambda }_{min,\mathcal{N}}\ge \underset{i\in \{1,\dots ,2n\}}{min}\left\{inf\left\{\left[{\mathcal{N}}_{i,i}\right]-\sum _{j=1,j\ne i}^{2n}\left|\left[{\mathcal{N}}_{i,j}\right]\right|\right\}\right\}.$$

(33)

If this is not efficient, all vertices of the interval matrix (30) could be tested for positive semi-definiteness. However, this is only recommended for sufficiently small values of n.

Remark 7.

The optimal, i.e., largest volume ellipsoid after fixing the matrix ${\mathcal{R}}_k$, for example, by the choice ${\mathcal{R}}_k={\mathbf{Q}}_k$, is obtained with the help of Theorem 4 if the parameter ${\alpha }_{\mathrm{O},k+1}>0$ is maximized by a line search procedure so that (32) remains satisfied. For ${\alpha }_{\mathrm{I},k+1}=0$ or non-regular matrices ${\tilde{\mathbf{A}}}_k^{-1}·\left[{\mathbf{A}}_k\right]$, the inner ellipsoidal bound is the empty set.

Using identical matrices ${\mathcal{R}}_k$ for the outer and inner ellipsoids, the thick ellipsoid sets is defined by

$${\left(\left(\mathcal{E}\right)\right)}_{k+1}={\left(\left(\mathcal{E}\right)\right)}_{k+1}\left(\mathbf{0},{\mathcal{R}}_k^{\frac{1}{2}}·{\tilde{\mathbf{A}}}_k,\left[{\alpha }_{\mathrm{I},k+1};{\alpha }_{\mathrm{O},k+1}\right]\right).$$

(34)

3.1.3. Illustrating Example

As an illustrating example, consider the system model

$${\mathbf{x}}_{k+1}=\left[\begin{array}{cc}0.5& p_1\\ p_2& 0.6\end{array}\right]·{\mathbf{x}}_k$$

(35)

with the independent parameters $p_1$ and $p_2$, where

$${\mathbf{x}}_k\in {\mathcal{E}}_k\mathrm{with}{\mathbf{Q}}_k=\left[\begin{array}{cc}1& 0\\ 0& 2\end{array}\right].$$

(36)

To show the influence of parameter uncertainty on the computed ellipsoidal enclosures, the following cases are distinguished in Figure 2:

case 1

$p_1\in \left[-1;2\right]$, $p_2\in \left[-1;2\right]$;

case 2

$p_1\in \left[-0.25;0.5\right]$, $p_2\in \left[-0.25;0.5\right]$;

case 3

$p_1\in \left[-0.1;0.2\right]$, $p_2\in \left[-0.1;0.2\right]$.

In the left column of Figure 2, the parameter matrix ${\mathcal{R}}_k={\mathbf{Q}}_k$ is chosen, while it is set to ${\mathcal{R}}_k={\tilde{\mathbf{A}}}_k^{-T}·{\mathbf{Q}}_k·{\tilde{\mathbf{A}}}_k^{-1}$ in the right column, where ${\tilde{\mathbf{A}}}_k$ is the interval midpoint of the uncertain system matrix.

Obviously, for the smaller uncertainty in the cases 2 and 3, the first setting is superior. This is not astonishing because it corresponds to that choice of matrix that would be employed for the covariance prediction in an Extended Kalman Filter (EKF). As it is well known, strong deviations from a point-valued system matrix ${\mathbf{A}}_k$ may turn an EKF quite inaccurate. Then, other choices of the forecasted covariance are often superior, as they could be determined by the Unscented Kalman Filter technique [27]. This approach was not investigated here. However, the almost symmetric interval bounds of $\left[{\mathbf{A}}_k\right]$ motivate the second choice of ${\mathcal{R}}_k$, which corresponds to ${\mathbf{Q}}_{k+1}={\mathbf{Q}}_k$.

Moreover, it should be pointed out that the case 1 in Figure 2a,b corresponds to a scenario in which the inner bound is empty, the case 2 in Figure 2c,d to a scenario in which the inner bound just appears due to the matrix $\left[{\mathbf{A}}_k\right]$ containing purely regular realizations; finally, Figure 2e,f shows the case of small uncertainty with a clearly visible inner solution set, where the proposed one-parameter optimization task from the previous subsections provides results that are close to the solution of the robust linear matrix inequality design, where ${\mathcal{R}}_k$ is not predefined but rather optimized with the help of the matrices ${\mathcal{P}}_{k+1}^{\mathrm{O}}$ and ${\mathcal{P}}_{k+1}^{\mathrm{I}}$ to obtain independent outer and inner enclosures (dashed and dotted lines, respectively).

3.2. Dikin Ellipsoids for the Intersection of Ellipsoids

The intersection of two ellipsoidal (state) domains is a fundamental operation for the estimation of state variables in a bounded-error framework. For that reason, this section introduces corresponding procedures that are readily applicable to thick ellipsoids, starting with the well-known case of two ellipsoids with identical midpoints. These results are further generalized to the case, where the ellipsoid midpoints may be different from each other.

3.2.1. Intersection of Ellipsoids with Identical Midpoints

Outer as well as inner ellipsoidal representations for the intersection of various ellipsoidal domains were investigated in numerous publications. Especially, the reference [49] gives a thorough overview of a large number of applicable techniques that can be cast into computationally feasible optimization problems. Basically, all of these routines are based on the use of linear matrix inequalities. At each intersection stage, such linear matrix inequalities need to be solved when computing so-called inner and outer bounds by means of a technique that is based on findings of Löwner and John [50]. Analogously, this holds for the Nerimovski approximation [51] that is compared with the first option in [49].

To avoid the necessity for solving linear matrix inequalities explicitly during online applications (note, these matrix inequalities are exploited for the offline proof of the validity of the following enclosure relations), a conservative approximation based on the so-called Dikin ellipsoids [49] is considered in this paper. These ellipsoids can be applied to describe inner and outer bounds for the intersection ${\mathcal{I}}_k$ of m ellipsoids

$${\mathcal{E}}_{i,k}=\left\{{\mathbf{x}}_k\in {\mathbb{R}}^n|{\mathbf{x}}_k^T·{\mathbf{Q}}_{i,k}^{-1}·{\mathbf{x}}_k\le 1,i\in \{1,\dots ,m\}\right\}$$

(37)

with an identical center located in the origin of the state space. The following relations for the shape matrices hold equally for identical non-zero center points.

Using the results of Lemma 5 in [49], the inner and outer bounds are given by

$${\mathcal{E}}_{\mathrm{D},\mathrm{I},k}\subseteq {\mathcal{I}}_k\subseteq {\mathcal{E}}_{\mathrm{D},\mathrm{O},k}$$

(38)

with ${\mathcal{E}}_{\mathrm{D},\mathrm{I},k}={\mathcal{E}}_{\mathrm{D},k}^2$ and ${\mathcal{E}}_{\mathrm{D},\mathrm{O},k}={\mathcal{E}}_{\mathrm{D},k}^{2m}$, respectively, where

$${\mathcal{E}}_{\mathrm{D},k}^r=\left\{{\mathbf{x}}_k\in {\mathbb{R}}^n|{\mathbf{x}}_k^T·r^{-1}·{\mathbf{Q}}_{\mathrm{D},k}^{-1}·{\mathbf{x}}_k\le 1\right\}$$

(39)

is defined for the common shape matrix

$${\mathbf{Q}}_{\mathrm{D},k}^{-1}=2\sum _{i=1}^m{\mathbf{Q}}_{i,k}^{-1}.$$

(40)

This definition of inner and outer ellipsoidal enclosures corresponds to a thick ellipsoid introduced in the Definitions 1 and 2 in the form

$${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}={\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}\left(\mathbf{0},{\mathbf{\Gamma }}_{\mathrm{D},k},\left[\sqrt{2};\sqrt{2m}\right]\right)\mathrm{with}{\mathbf{\Gamma }}_{\mathrm{D},k}={\mathbf{Q}}_{\mathrm{D},k}^{\frac{1}{2}}.$$

(41)

Although these enclosures are generally suboptimal, we restrict ourselves to these Dikin ellipsoids due to their computational inexpensive evaluation which is advantageous if real-time state estimation procedures are desired. A generalization to the case of ellipsoids with different midpoints is presented in the following subsection.

Efficient techniques that allow for checking whether the intersection of two ellipsoids is empty were published in [52,53]. This check is especially interesting if the thick ellipsoid state observer technique derived in the following section is applied to tasks in the area of fault detection. There, the case that a predicted state domain does not intersect with those domains that are compatible with sensor data typically indicates failures of system components, actuators, or sensors. This corresponds to the fact that the investigated (nominal) system model is proven to be infeasible.

A further related technique, however, not only relying on ellipsoidal state enclosures was very recently derived in [54]. There, outer bounds for the domains of reachable states are characterized by ellipsoids, while disturbances are bounded by zonotopes that may become degenerate if one or more coordinates are unbounded. This property of degeneracy will be picked up again in Section 3.2.3, where an example is presented in which only partial state measurements are available in the thick ellipsoid framework.

3.2.2. Generalization to the Intersection of Ellipsoids with Different Midpoints

An extension of the technique for computing Dikin ellipsoids to a scenario in which two ellipsoids with different centers are intersected, is based on the following three-stage procedure:

Step 1

Determine the common center point for the desired inner and outer bounds of the intersection that must be included in all ellipsoids to be intersected;

Step 2

Determine initial approximations of the shape matrices for the inner and outer bounds according to Section 3.2.1;

Step 3

For non-empty inner bounds, correct the outer enclosure so that the inner and outer ellipsoid surfaces become parallel to each other and, thus, form a thick ellipsoid $\left(\left(\mathcal{E}\right)\right)$.

Step 3*

As an alternative to Step 3, the initial outer enclosure remains fixed and the inner ellipsoid surface is adapted to become parallel to each other to form a thick ellipsoid $\left(\left(\mathcal{E}\right)\right)$.

As a heuristic approach for the computation of the common center point ${\tilde{\mathit{\mu }}}_k$ in Step 1 of two ellipsoids that are described by the midpoints ${\mathit{\mu }}_{1,k}$ and ${\mathit{\mu }}_{2,k}$ as well as the shape matrices ${\mathbf{Q}}_{1,k}$ and ${\mathbf{Q}}_{2,k}$, respectively, the first one is interpreted as the prior knowledge in the innovation step of a Kalman filter [25,55] and the second one as the corresponding information of the state measurement.

Under this assumption, one obtains the Kalman gain matrix

$${\mathbf{L}}_k={\mathbf{Q}}_{1,k}·{\left({\mathbf{Q}}_{1,k}+{\mathbf{Q}}_{2,k}\right)}^{-1}$$

(42)

with the updated ellipsoid midpoint

$${\tilde{\mathit{\mu }}}_k={\mathit{\mu }}_{1,k}+{\mathbf{L}}_k·\left({\mathit{\mu }}_{2,k}-{\mathit{\mu }}_{1,k}\right)$$

(43)

that is now used for the solution of the Steps 2 and 3 listed above.

For both ellipsoids $i\in \{1,2\}$ to be intersected, scaling factors ${\xi }_i>$ and ${\zeta }_i>$ are determined which represent the minimum and maximum distances of the new midpoint from the ellipsoid surface according to

$${\xi }_{i,k}=1-∥\Delta {\mathit{\mu }}_{i,k}∥\mathrm{and}{\zeta }_i=1+∥\Delta {\mathit{\mu }}_{i,k}∥,$$

(44)

with

$$\Delta {\mathit{\mu }}_{i,k}={\mathbf{Q}}_{i,k}^{-\frac{1}{2}}·\left({\tilde{\mathit{\mu }}}_k-{\mathit{\mu }}_{i,k}\right)$$

(45)

and $∥\Delta {\mathit{\mu }}_{i,k}∥$ as the Euclidean norm of the vector-valued argument.

In such a way, both ellipsoids to be intersected can be enclosed by thick ellipsoids

$${\left(\left(\mathcal{E}\right)\right)}_{i,k}={\left(\left(\mathcal{E}\right)\right)}_{i,k}\left({\tilde{\mathit{\mu }}}_k,{\tilde{\mathbf{\Gamma }}}_{i,k},\left[{\xi }_{i,k};{\zeta }_{i,k}\right]\right)\mathrm{with}{\tilde{\mathbf{\Gamma }}}_{i,k}={\mathbf{Q}}_{i,k}^{\frac{1}{2}}.$$

(46)

The initial inner approximation of the shape matrix of the resulting intersection is given by

$${\mathbf{Q}}_k^{\mathrm{I}}=2·{\left(2·\left({\left({\xi }_{1,k}^2·{\mathbf{Q}}_{1,k}\right)}^{-1}+{\left({\xi }_{2,k}^2·{\mathbf{Q}}_{2,k}\right)}^{-1}\right)\right)}^{-1}$$

(47)

according to the previous subsection, while its outer counterpart results in

$${\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}=4·{\left(2·\left({\left({\zeta }_{1,k}^2·{\mathbf{Q}}_{1,k}\right)}^{-1}+{\left({\zeta }_{2,k}^2·{\mathbf{Q}}_{2,k}\right)}^{-1}\right)\right)}^{-1}.$$

(48)

Usually, these two shape matrices are not proportional to each other in an element-wise sense. Hence, the inner shape matrix ${\mathbf{Q}}_k^{\mathrm{I}}$ is inflated in Step 3 according to the following theorem so that it contains the outer ellipsoid with certainty.

Theorem 5

(Thick interval representation of the Dikin intersection of ellipsoids with different midpoints). The thick Dikin ellipsoid

$${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}={\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}\left({\tilde{\mu }}_k,{\mathbf{\Gamma }}_{\mathrm{D},k},\left[1;{\left(\underset{i\in \{1,\dots ,n\}}{min}\left\{{\lambda }_i\left({\left({\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}\right)}^{-1}·{\mathbf{Q}}_k^{\mathrm{I}}\right)\right\}\right)}^{-\frac{1}{2}}\right]\right)with{\mathbf{\Gamma }}_{\mathrm{D},k}={\left({\mathbf{Q}}_k^{\mathrm{I}}\right)}^{\frac{1}{2}}.$$

(49)

represents guaranteed inner and outer bounds of the intersection of two ellipsoids with different midpoints if all parameters are chosen according to Equations (42)–(48).

Proof. 

The validity of the inner bound is a direct consequence of Lemma 5 in [49], and the construction of the matrix (45) as a guaranteed underapproximation of the domains to be intersected. The outer bound of (49) is represented by the shape matrix

$${\mathbf{Q}}_k^{\mathrm{O}}={\alpha }_{\mathrm{O},k+1}^2·{\mathbf{Q}}_k^{\mathrm{I}}.$$

(50)

It results from Equation (16) after setting ${\mathbf{A}}_k=\mathbf{I}$, ${\tilde{\mathbf{A}}}_k=\mathbf{I}$, ${\mathbf{Q}}_k={\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}$, and ${\mathcal{R}}_k={\mathbf{Q}}_k^{\mathrm{I}}$ which yields

$${\alpha }_{\mathrm{O},k+1}^{-2}·\mathbf{I}\le {\left({\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}\right)}^{-1}·{\mathbf{Q}}_k^{\mathrm{I}}.$$

(51)

This inequality is satisfied if

$${\alpha }_{\mathrm{O},k+1}^{-2}\le \underset{i\in \{1,\dots ,n\}}{min}\left\{{\lambda }_i\left({\left({\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}\right)}^{-1}·{\mathbf{Q}}_k^{\mathrm{I}}\right)\right\},$$

(52)

where the expression in curly brackets denotes the minimum eigenvalue of the corresponding matrix product. Taking the square root of the inverse of this expression completes the proof of Theorem 5. □

The alternative solution according to Step 3* is given by the thick ellipsoid enclosure

$${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}^{*}={\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}^{*}\left({\tilde{\mathit{\mu }}}_k,{\mathbf{\Gamma }}_{\mathrm{D},k}^{*},\left[{\left(\underset{i\in \{1,\dots ,n\}}{max}\left\{{\lambda }_i\left({\left({\mathbf{Q}}_k^{\mathrm{I}}\right)}^{-1}·{\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}\right)\right\}\right)}^{-\frac{1}{2}};1\right]\right)\mathrm{with}{\mathbf{\Gamma }}_{\mathrm{D},k}^{*}={\left({\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}\right)}^{\frac{1}{2}},$$

(53)

which is a direct consequence of Theorem 5 by searching for an ellipsoid parallel to the outer bound that is guaranteed to belong to the interior of the one with the shape matrix ${\mathbf{Q}}_k^{\mathrm{I}}$.

Remark 8.

Due to the fact that the relation

$$\underset{i\in \{1,\dots ,n\}}{max}\left\{{\lambda }_i\left({\left({\mathbf{Q}}_k^{\mathrm{I}}\right)}^{-1}·{\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}\right)\right\}=\underset{i\in \{1,\dots ,n\}}{min}{\left\{{\lambda }_i\left({\left({\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}\right)}^{-1}·{\mathbf{Q}}_k^{\mathrm{I}}\right)\right\}}^{-1}>0$$

(54)

holds for arbitrary positive definite matrices ${\mathbf{Q}}_k^{\mathrm{I}}$ and ${\mathbf{Q}}_k^{{\mathrm{O}}^{\prime }}$, the ratio between the volumes of the outer and inner ellipsoids in Equations (49) and (53) is identical. Hence, the variant (49) is typically preferred if thick ellipsoids with maximum volume inner enclosures are desired, while (53) provides tighter outer bounds.

Remark 9.

For ellipsoids with identical midpoints ${\tilde{\mathit{\mu }}}_k={\mathit{\mu }}_{1,k}={\mathit{\mu }}_{2,k}$, the ellipsoid of Theorem 5 is identical to the one in (41) if the midpoint parameter ${\tilde{\mathit{\mu }}}_k$ is used instead of the value $\mathbf{0}$ in (41).

3.2.3. Illustrating Example

In this subsection, various examples for the intersection of ellipsoids and the result interpretation by means of the thick ellipsoids ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}$ and ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}^{*}$, respectively, are summarized. Figure 3a deals with the application of the intersection procedure of Section 3.2.1 for two ellipsoids identically centered at the origin with the shape matrices

$${\mathbf{Q}}_{1,k}=\left[\begin{array}{cc}1& 0\\ 0& 8\end{array}\right]\mathrm{and}{\mathbf{Q}}_{2,k}=\left[\begin{array}{cc}4& -5\\ -5& 8\end{array}\right].$$

(55)

Here, both enclosures ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}$ and ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}^{*}$ are identical. Moreover, the fact that the ellipsoid midpoints are identical results in inner and outer enclosures ${\mathcal{E}}_{\mathrm{D},\mathrm{k}}^{\mathrm{I}}$ and ${\mathcal{E}}_{\mathrm{D},\mathrm{k}}^{\mathrm{O}}$ that are quite close to the optimal ellipsoidal result representation.

To visualize the procedure of Section 3.2.2 for the case of two different midpoints

$${\mathit{\mu }}_{1,k}={\left[\begin{array}{cc}0& 0\end{array}\right]}^T\mathrm{and}{\mathit{\mu }}_{2,k}={\left[\begin{array}{cc}1& 2\end{array}\right]}^T$$

(56)

either with the shape matrices (55) or with the new shape matrices

$${\mathbf{Q}}_{1,k}=\left[\begin{array}{cc}1& 0\\ 0& 8\end{array}\right]\mathrm{and}{\mathbf{Q}}_{2,k}=4·\left[\begin{array}{cc}4& -5\\ -5& 8\end{array}\right],$$

(57)

the results in Figure 3c,d are presented. These results confirm the observation described in Remark 8 that the use of Equation (49) leads to tighter inner bounds, while Equation (53) possesses the better outer bounds. Depending on the application to detect states belonging to the solution set with certainty or proving in a guaranteed way that specific state domains are not reachable in the frame of a safety or feasibility test, either the first or the second option should be preferred.

However, the observed larger distance between the resulting inner and outer bounds in comparison with Figure 3a gives rise to the assumption that further improvements of the solution could be possible. This could either be done by means of the procedures in [49] or by intersecting the newly obtained outer bounds of both ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}$ and ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}^{*}$. Note, also the choice of the common ellipsoid midpoint represents a degree of freedom that can be investigated in future work together with possibilities to intersect the enlarged dashed-dotted bounds (which correspond to the outer ellipsoid surfaces of (46)) in a recursive manner with both ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}$ and ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{D},k}^{*}$.

Interestingly, the proposed intersection procedures can also be applied to the case where one of the original ellipsoids becomes degenerate. In practice, this is the case for a pure measurement of a single state variable. This case is considered in Figure 4, where in all subplots the matrices ${\mathbf{Q}}_{1,k}$ from Equations (55) and (57) were replaced with ${\mathbf{Q}}_{1,k}^{-1}=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]$, corresponding to $\left|x_{1,k}\right|\le 1$.

Due to the fact that the basic Formula (40) for computing Dikin ellipsoids only makes use of the inverted shape matrices, the solution procedure remains applicable for the degenerate case. This is a crucial requirement to use this intersection approach in cases in which only a subset to the state variables appears in the system’s output equation when considering state estimation procedures according to the following two sections. In Figure 4b,c, the new ellipsoid midpoint ${\tilde{\mathit{\mu }}}_k$ results from the following adapted Kalman gain

$${\mathbf{L}}_k={\mathbf{Q}}_{2,k}·\left[\begin{array}{c}1\\ 0\end{array}\right]·{\left(\left[\begin{array}{cc}1& 0\end{array}\right]·{\mathbf{Q}}_{2,k}·\left[\begin{array}{c}1\\ 0\end{array}\right]+{\sigma }_{1,k}^2\right)}^{-1},\mathrm{where}{\mathbf{Q}}_{1,k}^{-1}=\left[\begin{array}{cc}{\sigma }_{1,k}^{-2}& 0\\ 0& 0\end{array}\right]$$

(58)

with the updated ellipsoid midpoint

$${\tilde{\mathit{\mu }}}_k={\mathit{\mu }}_{2,k}+{\mathbf{L}}_k·\left(\left[\begin{array}{cc}1& 0\end{array}\right]·\left({\mathit{\mu }}_{1,k}-{\mathit{\mu }}_{2,k}\right)\right),$$

(59)

for which all remaining parameter values are defined in (56) and (57).

3.3. Thick Ellipsoid Union of Two Ellipsoids with Different Midpoints

3.3.1. General Solution Procedure

The thick ellipsoid union of two ellipsoids with different midpoints is a direct consequence of the procedure in Section 3.2.2. Firstly, the shape matrix ${\mathbf{Q}}_k^{\mathrm{I}}$ of the thick ellipsoid

$${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{U},k}={\left(\left(\mathcal{E}\right)\right)}_{\mathrm{U},k}\left({\tilde{\mathit{\mu }}}_k,{\mathbf{\Gamma }}_{\mathrm{U},k},\left[1;{\eta }_k\right]\right)\mathrm{with}{\eta }_k={\eta }_k^{*}\ge 1\mathrm{and}{\mathbf{\Gamma }}_{\mathrm{U},k}={\left({\mathbf{Q}}_k^{\mathrm{I}}\right)}^{\frac{1}{2}}$$

(60)

is determined as in Equation (47). The remaining degree of freedom ${\eta }_k^{*}$ can be found according to the work of John [50]. For that purpose, specific points on the ellipsoids to be merged were extracted in [24] with subsequently inflating the inner bound by the factor ${\eta }_k^{*}\ge 1$ until all points are included in its interior.

To avoid the extraction of individual points, which has to be performed with sufficient care so that the new outer bound is guaranteed to be conservative, Equation (50) with ${\zeta }_{j,k}$ given in (44) is generalized to find the maximum scaling

$${\eta }_k^{*}=\underset{j\in \{1,2\}}{max}\left\{{\left(\underset{i\in \{1,\dots ,n\}}{min}\left\{{\lambda }_i\left({\left({\zeta }_{j,k}^2·{\mathbf{Q}}_{j,k}\right)}^{-1}·{\mathbf{Q}}_k^{\mathrm{I}}\right)\right\}\right)}^{-\frac{1}{2}}\right\}$$

(61)

so that all outer bounds of the thick ellipsoids in (46) are guaranteed to be contained in the union.

Remark 10.

This inflation operation generalizes naturally to an arbitrary integer number j of ellipsoids to be merged after they have been converted into representations with a common midpoint.

As an alternative parameterization (similar to the computation of the intersection of ellipsoids), both outer and inner bounds can be tuned by an alternative option. For that purpose, firstly define an approximation to the resulting shape matrix of the union according to the covariance prediction of a Kalman filter, where the state equation in the Kalman filter corresponds to computing the mean value of the data from each of the ellipsoids, with

$${\tilde{\mathbf{Q}}}_k=\frac{1}{4}·\left({\mathbf{Q}}_{1,k}+{\mathbf{Q}}_{2,k}\right).$$

(62)

Then, the scaling factors according to (49) and (53) yield the shape matrices

$${\mathbf{Q}}_k^{\mathrm{O}},*={\tilde{\mathbf{Q}}}_k·{\left({\eta }_k^{\mathrm{O}}\right)}^2,{\eta }_k^{\mathrm{O}}=\underset{j\in \{1,2\}}{max}\left\{{\left(\underset{i\in \{1,\dots ,n\}}{min}\left\{{\lambda }_i\left({\left({\zeta }_{j,k}^2·{\mathbf{Q}}_{j,k}\right)}^{-1}·{\tilde{\mathbf{Q}}}_k\right)\right\}\right)}^{-\frac{1}{2}}\right\}$$

(63)

for the outer bound as well as

$${\mathbf{Q}}_k^{\mathrm{I}},*={\tilde{\mathbf{Q}}}_k·{\left({\eta }_k^{\mathrm{I}}\right)}^2,{\eta }_k^{\mathrm{I}}={\left(\underset{i\in \{1,\dots ,n\}}{max}\left\{{\lambda }_i\left({\left({\mathbf{Q}}_k^{\mathrm{I}}\right)}^{-1}·{\tilde{\mathbf{Q}}}_k\right)\right\}\right)}^{-\frac{1}{2}}$$

(64)

for the inner bound.

Together, they result in the thick ellipsoid

$${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{U},k}^{*}={\left(\left(\mathcal{E}\right)\right)}_{\mathrm{U},k}^{*}\left({\tilde{\mathit{\mu }}}_k,{\mathbf{\Gamma }}_{\mathrm{U},k}^{*},\left[{\eta }_k^{\mathrm{I}};{\eta }_k^{\mathrm{O}}\right]\right)\mathrm{with}{\eta }_k={\eta }_k^{*}\ge 1\mathrm{and}{\mathbf{\Gamma }}_{\mathrm{U},k}^{*}={\left({\tilde{\mathbf{Q}}}_k\right)}^{\frac{1}{2}}.$$

(65)

3.3.2. Illustrating Example

Using the same examples as in (55)–(57), the procedure from Section 3.3.1 is visualized in Figure 5. Due to the fact that the outer hull in Figure 5c is more conservative than necessary (computed by means of Equation (60)), the alternative procedure from Equation (65) is preferable in the considered example. More complex options for computing optimal outer ellipsoidal enclosures typically need to employ the solution of matrix inequalities and/or minmax optimization problems in each evaluation step. Due to the resulting computational effort, such options are not further considered for the desired online application in this paper. However, the reader is referred to [19,30,56] for related work.

4. Thick Ellipsoid State Estimation Algorithm

As shown in Figure 6, the thick ellipsoid state estimation procedure consists of two stages. Prediction steps, which are based on the (uncertain) discrete-time state equation, are executed up to the point of time at which a new measurement becomes available. There, a measurement-based correction of the a-priori knowledge resulting from the previous state prediction is performed. This corrected state information then serves in a recursive manner as the input for a subsequent state prediction.

4.1. Thick Ellipsoid Prediction Step

For state estimation purposes, typically thick ellipsoids have to be considered in the prediction step according to Section 3.1 which have a non-zero midpoint. Their propagation then consists of separating the state equation as illustrated in Figure 7 into a part that depends on an ellipsoid centered at the origin and an offset term that accounts for the non-zero center according to

$${\mathbf{x}}_{k+1}=\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\check{\mathbf{x}}}_k+{\tilde{\mathbf{A}}}_k·{\mathit{\mu }}_k+\left(\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)-{\tilde{\mathbf{A}}}_k\right)·{\mathit{\mu }}_k,$$

(66)

where

$$\begin{array}{cc}\hfill {\mathbf{x}}_k& \in {\left(\left(\mathcal{E}\right)\right)}_k\left({\mathit{\mu }}_k,{\mathbf{\Gamma }}_k,\left[{\underset{¯}{\rho }}_k;{\overline{\rho }}_k\right]\right),\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill {\check{\mathbf{x}}}_k& \in {\check{\left(\left(\mathcal{E}\right)\right)}}_k\left(\mathbf{0},{\mathbf{\Gamma }}_k,\left[{\underset{¯}{\rho }}_k;{\overline{\rho }}_k\right]\right),\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{A}}}_k& =\mathbf{A}\left({\mathit{\mu }}_k,\mathrm{mid}\left(\left[{\mathbf{p}}_k\right]\right)\right),\mathrm{and}\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill {\mathbf{p}}_k& \in \left[{\mathbf{p}}_k\right]=\left[{\underset{¯}{\mathbf{p}}}_k;{\overline{\mathbf{p}}}_k\right]\mathrm{with}\mathrm{mid}\left(\left[{\mathbf{p}}_k\right]\right)=\frac{1}{2}·\left({\underset{¯}{\mathbf{p}}}_k+{\overline{\mathbf{p}}}_k\right).\hfill \end{array}$$

(70)

Then, the following steps are executed during the state prediction:

• Propagate the inner bound of the thick ellipsoid (68) and extract the inner hull of the image set that is obtained by applying the mapping

  $${\check{\mathbf{x}}}_{k+1}=\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)·{\check{\mathbf{x}}}_k.$$

  (71)

  Denote the corresponding shape matrix of the inner ellipsoid by ${\check{\mathbf{Q}}}_{k+1}^{\mathrm{I}}$. This matrix is related to the shape matrix ${\alpha }_{\mathrm{I},k+1}^2·{\mathbf{Q}}_{k+1}$ in Equation (25) by

  $${\check{\mathbf{Q}}}_{k+1}^{\mathrm{I}}={\alpha }_{\mathrm{I},k+1}^2·{\rho }_k^2·{\mathbf{\Gamma }}_k·{\mathbf{\Gamma }}_k^T.$$

  (72)
• Propagate the outer bound of the thick ellipsoid (68) and extract the outer hull of the image set that is obtained by applying the mapping (71). Denote the corresponding shape matrix of the outer ellipsoid by ${\check{\mathbf{Q}}}_{k+1}^{\mathrm{O}}$. This matrix is related to the shape matrix ${\alpha }_{\mathrm{O},k+1}^2·{\mathbf{Q}}_{k+1}$ in Equation (9) by

  $${\check{\mathbf{Q}}}_{k+1}^{\mathrm{O}}={\alpha }_{\mathrm{O},k+1}^2·{\overline{\rho }}_k^2·{\mathbf{\Gamma }}_k·{\mathbf{\Gamma }}_k^T.$$

  (73)
• Compute interval bounds for the term

  $${\mathbf{b}}_k=\left(\mathbf{A}\left({\mathbf{x}}_k,{\mathbf{p}}_k\right)-{\tilde{\mathbf{A}}}_k\right)·{\mathit{\mu }}_k\in \left[{\mathbf{b}}_k\right]$$

  (74)

  with ${\mathbf{x}}_k$, ${\tilde{\mathbf{A}}}_k$, and ${\mathbf{p}}_k$ defined according to (67), (69), and (70) and deflate the inner ellipsoid bound as well as inflate the outer bound according to the procedure described in [20,24] with

  $$\begin{array}{cc}\hfill {\mathbf{Q}}_{k+1}^{\mathrm{I}}& ={\left(1-{\rho }_{\mathrm{I},k+1}\right)}^2·{\check{\mathbf{Q}}}_{k+1}^{\mathrm{I}},{\rho }_{\mathrm{I},k+1}=sup\left\{∥{\alpha }_{\mathrm{I},k+1}^{-1}·{\underset{¯}{\rho }}_k^{-1}·{\mathbf{\Gamma }}_k^{-1}·\left[{\mathbf{b}}_k\right]∥\right\}\hfill \end{array}$$

  (75)

  and

  $$\begin{array}{cc}\hfill 1-1{\mathbf{Q}}_{k+1}^{\mathrm{O}}& ={\left(1+{\rho }_{\mathrm{O},k+1}\right)}^2·{\check{\mathbf{Q}}}_{k+1}^{\mathrm{O}},{\rho }_{\mathrm{O},k+1}=sup\left\{∥{\alpha }_{\mathrm{O},k+1}^{-1}·{\overline{\rho }}_k^{-1}·{\mathbf{\Gamma }}_k^{-1}·\left[{\mathbf{b}}_k\right]∥\right\}.\hfill \end{array}$$

  (76)

  Note, for ${\rho }_{\mathrm{I}}\ge 1$, the inner bound becomes the empty set.
• Compute the updated ellipsoid midpoint as

  $${\mathit{\mu }}_{k+1}={\tilde{\mathbf{A}}}_k·{\mathit{\mu }}_k.$$

  (77)
• Finally, determine the predicted thick ellipsoid set

  $${\mathbf{x}}_{k+1}\in {\left(\left(\mathcal{E}\right)\right)}_{k+1}\left({\mathit{\mu }}_{k+1},{\mathbf{\Gamma }}_{k+1},\left[{\rho }_{k+1};{\overline{\rho }}_{k+1}\right]\right),$$

  (78)

  where

  $$\begin{array}{cc}\hfill {\rho }_{k+1}& ={\rho }_k·{\alpha }_{\mathrm{I},k+1}·\left(1-{\rho }_{\mathrm{I},k+1}\right),\hfill \\ \hfill {\overline{\rho }}_{k+1}& ={\overline{\rho }}_k·{\alpha }_{\mathrm{O},k+1}·\left(1+{\rho }_{\mathrm{O},k+1}\right),\mathrm{and}\hfill \\ \hfill {\mathbf{\Gamma }}_{k+1}& ={\tilde{\mathbf{A}}}_k·{\mathbf{\Gamma }}_k.\hfill \end{array}$$

  (79)

Remark 11.

In cases where the inner ellipsoid becomes the empty set, purely the outer ellipsoid hull is further propagated to represent worst-case outer state enclosures.

4.2. Thick Ellipsoid Correction Step

The correction step of the thick ellipsoid state estimator at the point $k+1$ is given by the application of Theorem 5 with the following step-by-step procedure (in this description, we prefer the version of maximized inner bounds to reduce the likelihood of empty inner enclosures; however, the task of minimizing the outer volume can be stated in full analogy):

• Determine the inner shape matrix on the basis of Equation (47) according to

  $${\mathbf{Q}}_{k+1}^{{\mathrm{I}}^{\prime }}=2·{\left(2·\left({\left({\xi }_{1,k+1}^2·{\mathbf{Q}}_{k+1}^{\mathrm{I}}\right)}^{-1}+{\left({\xi }_{2,k+1}^2·{\mathbf{Q}}_{\mathrm{m},k+1}\right)}^{-1}\right)\right)}^{-1},$$

  (80)

  where (based on the change of the ellipsoids’ midpoint positions according to (43)) ${\mathbf{Q}}_{k+1}^{\mathrm{I}}$ from (75) is the inner bound of the previous prediction; ${\mathbf{Q}}_{\mathrm{m},k+1}$ characterizes the measurement uncertainty (possibly given in terms of a degenerate ellipsoid).
• Use Equation (48) to determine the intermediate outer bound

  $${\mathbf{Q}}_{k+1}^{{\mathrm{O}}^{\prime }}=4·{\left(2·\left({\left({\zeta }_{1,k+1}^2·{\mathbf{Q}}_{k+1}^{\mathrm{O}}\right)}^{-1}+{\left({\zeta }_{2,k+1}^2·{\mathbf{Q}}_{\mathrm{m},k+1}\right)}^{-1}\right)\right)}^{-1}$$

  (81)

  with ${\mathbf{Q}}_{k+1}^{\mathrm{O}}$ from (76).
• Finally, Equation (49) yields the thick Dikin ellipsoid as the result of the measurement-based correction step.

4.3. Visualization of the Thick Ellipsoid State Estimation Procedure

To visualize the individual steps of the thick ellipsoid state estimation procedure, the test example

$$\left[\begin{array}{c}{x}_{1,k+1}\\ x_{2,k+1}\end{array}\right]=\left[\begin{array}{cc}1& 0.05\\ -0.01x_{1,k}& 0.95\end{array}\right]·\left[\begin{array}{c}{x}_{1,k}\\ x_{2,k}\end{array}\right]$$

(82)

is investigated. At the initial point of time $k=0$, all possible system states are included in the set

$${\left(\left(\mathcal{E}\right)\right)}_0={\left(\left(\mathcal{E}\right)\right)}_0\left({\mathit{\mu }}_0,{\mathbf{\Gamma }}_0,\left[1;1\right]\right)\mathrm{with}{\mathit{\mu }}_0=\left[\begin{array}{c}1\\ 1\end{array}\right]\mathrm{and}{\mathbf{\Gamma }}_0={\left(\left[\begin{array}{cc}1& 0.9\\ 0.9& 1\end{array}\right]\right)}^{\frac{1}{2}}.$$

(83)

The corresponding ellipsoid is shown as the dash-dotted line in Figure 8a. Performing a single prediction step leads to a shift in the ellipsoid midpoint according to Equation (77). Moreover, the nonlinearity in the system model (82) causes the resulting distance between the inner and outer bounds of the predicted ellipsoid ${\left(\left(\mathcal{E}\right)\right)}_1$. Now, this ellipsoid is intersected with the measurement information that is highlighted by the vertical strip in Figure 8a,b. According to Figure 8b, the correction step leads again to a shift of the midpoint of the new thick ellipsoid ${\left(\left(\mathcal{E}\right)\right)}_1^{\prime }$ according to Equations (58) and (59). The inner bound of ${\left(\left(\mathcal{E}\right)\right)}_1^{\prime }$ is always fully contained in the band of possible state measurements ${\left(\left(\mathcal{E}\right)\right)}_{\mathrm{m},1}$, while its outer bound becomes tighter than the result of the preceding state prediction.

Now, the result ${\left(\left(\mathcal{E}\right)\right)}_1^{\prime }$ serves as the input for the subsequent state prediction. For the sake of a compact notation and for compatibility with Figure 6 and Figure 7, the prime symbol is suppressed in Figure 8c, where the prediction from $k=1$ to $k=2$ is shown. After this prediction, the next correction step can be executed according to Figure 8d. Note, tight measurement tolerances or strongly nonlinear dependencies in the system model (re-written into the quasi-linear form) may lead to empty inner bounds after either the prediction or the correction step. Then, a slightly modified algorithm could be used. Instead of interpreting the inner bound as the state domain that is guaranteed to belong to the reachable solution set after the execution of all computations until the current instant k, it can be used to determine a measure for the precision of the sequence of prediction and correction steps at a single point of time k. For that purpose, it is necessary to initialize the inner bound of ${\left(\left(\mathcal{E}\right)\right)}_k$ with the same scaling parameter as the outer bound (cf. the interval parameter with identical bounds in (83)) and to use it in analogy to the initialization of the ellipsoid in Figure 8a.

5. Application Scenario: State and Disturbance Estimation for an Underactuated Hovercraft

As a more complex benchmark application, consider the motion of an underactuated hovercraft in its body-fixed coordinate frame [57,58]. In these coordinates, it is desired to employ the ellipsoid techniques summarized in the previous sections to compute guaranteed outer state enclosures in the presence of bounded noise and disturbances. Moreover, the thick ellipsoid state estimation algorithm of Section 4 is employed to reconstruct bounds on external disturbances.

5.1. Modeling

According to Figure 9, denote the surge and sway velocities by $x_1$ and $x_2$, respectively, and the hovercraft’s yaw rate by $x_3$ in a body-fixed coordinate frame. For the following derivation of the equations of motion, assume that the translational velocities are summarized in the vector $\mathbf{V}={\left[\begin{array}{ccc}{x}_1& x_2& 0\end{array}\right]}^T$, while the angular velocity vector $\mathbf{W}$ with respect to the ship’s center of gravity (COG) is given by $\mathbf{W}={\left[\begin{array}{ccc}0& 0& x_3\end{array}\right]}^T$. In this formulation, the dynamics associated with the motion in heave, roll, and pitch are neglected by setting the corresponding vector entries to zero.

As also shown in [59,60], a force and torque balance in all degrees of freedom yields the equations of motion

$$\begin{array}{cc}\hfill & \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial E}{\partial \mathbf{V}}\right)+\mathbf{W}\times \left(\frac{\partial E}{\partial \mathbf{V}}\right)=\mathbf{F}\mathrm{and}\hfill \\ \hfill & \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial E}{\partial \mathbf{W}}\right)+\mathbf{V}\times \left(\frac{\partial E}{\partial \mathbf{V}}\right)+\mathbf{W}\times \left(\frac{\partial E}{\partial \mathbf{W}}\right)=\mathbf{T},\hfill \end{array}$$

(84)

where

$$E=\frac{1}{2}·\left({\mathbf{V}}^T·\left[\begin{array}{ccc}{m}_{11}& 0& 0\\ 0& m_{22}& 0\\ 0& 0& m_{33}\end{array}\right]·\mathbf{V}+{\mathbf{W}}^T·\left[\begin{array}{ccc}{J}_{11}& 0& 0\\ 0& J_{22}& 0\\ 0& 0& J_{33}\end{array}\right]·\mathbf{W}\right)$$

(85)

is the total kinetic energy of the ship (under the assumption of purely diagonal mass and inertia matrices [60]) with the strictly non-negative entries $m_{ii}\ge 0$ and $J_{ii}\ge 0$, $i\in \{1,2,3\}$. For a hovercraft, it can be assumed that $m_{11}=m_{22}=m$ corresponds to the mass of the ship, without additional hydrodynamic effects, while $J_{33}=J_{\mathrm{r}}$ is the mass moment of inertia of the vessel in the direction of $x_3$.

Defining further the vectors

$$\mathbf{F}=\left[\begin{array}{c}{u}_1-c{x}_1-d_1\\ -c{x}_2-d_2\\ 0\end{array}\right]\mathrm{and}\mathbf{T}=\left[\begin{array}{c}0\\ 0\\ u_2-c_{\mathrm{r}}{x}_3-d_3\end{array}\right]$$

(86)

of non-conservative, external generalized forces and torques with the control inputs $u_1$ and $u_2$ and independent, velocity-proportional damping effects in each coordinate with the coefficient $c>0$ for both translational degrees of freedom as well as $c_{\mathrm{r}}>0$ for the rotary motion, the equations of motion (84) simplify to

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\\ {\dot{d}}_1\\ {\dot{d}}_2\\ {\dot{d}}_3\end{array}\right]=\left[\begin{array}{cccccc}-\frac{c}{m}& x_3& 0& \frac{1}{m}& 0& 0\\ -x_3& -\frac{c}{m}& 0& 0& \frac{1}{m}& 0\\ 0& 0& -\frac{c_{\mathrm{r}}}{J}& 0& 0& \frac{1}{J}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]·\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ d_1\\ d_2\\ d_3\end{array}\right]+\left[\begin{array}{cc}\frac{1}{m}& 0\\ 0& 0\\ 0& \frac{1}{J}\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]·\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right],$$

(87)

where $d_i$, $i\in \{1,2,3\}$, represent external disturbances according to Figure 9 due to wind and currents. These disturbances are included in the state equations (87) by means of integrator disturbance models.

For the following simulation, a linear state feedback controller

$$\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]=-\left[\begin{array}{ccc}{k}_{11}& k_{12}& k_{13}\\ k_{21}& k_{22}& k_{23}\end{array}\right]·\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right],$$

(88)

corresponding to a feedback of the actual velocity, is parameterized by means of pole assignment to decelerate the hovercraft up to standstill at the desired operating point $x_1=x_2=x_3=0$.

For the use of the proposed ellipsoidal state estimation scheme, that will be applied to the closed-loop system, the state equations are discretized in time with the step size T according to

$$\left[\begin{array}{c}{x}_{1,k+1}\\ x_{2,k+1}\\ x_{3,k+1}\\ d_{1,k+1}\\ d_{2,k+1}\\ d_{3,k+1}\end{array}\right]=\left(\mathbf{I}+T·\left[\begin{array}{cccccc}-\frac{c+k_{11}}{m}& x_{3,k}-\frac{k_{12}}{m}& -\frac{k_{13}}{m}& -\frac{1}{m}& 0& 0\\ -x_{3,k}& -\frac{c}{m}& 0& 0& -\frac{1}{m}& 0\\ -\frac{k_{21}}{J}& -\frac{k_{22}}{J}& -\frac{c_{\mathrm{r}}+k_{23}}{J}& 0& 0& \frac{1}{J}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\right)·\left[\begin{array}{c}{x}_{1,k}\\ x_{2,k}\\ x_{3,k}\\ d_{1,k}\\ d_{2,k}\\ d_{3,k}\end{array}\right].$$

(89)

The (normalized) parameters used for the following simulations are: $c=0.01$, $m=5$, $J=0.8$, $c_{\mathrm{r}}=0.01$, $k_{11}=0.5$, $k_{12}=1.5$, $k_{13}=0$, $k_{21}=0$, $k_{22}=0$, $k_{23}=0.1$, and $T=0.1$.

As measured system outputs, all three velocities (surge, sway, and yaw) are assumed to be available. For the simulation, they are described by the three independent sensor models

$$\frac{1}{{\delta }_i^2}·{\left(x_{i,k}-y_{\mathrm{m}i,k}\right)}^2\le 1,$$

(90)

where $y_{\mathrm{m}i,k}$, $i\in \{1,2,3\}$, are the data provided by the velocity sensors and ${\delta }_i>0$ their corresponding maximum uncertainty. Besides filtering the velocity signals, the proposed observer is used to reconstruct the ranges of the disturbance inputs $d_i$ that are compatible with the system model and the sensor data. Due to the assumption of constant but uncertain disturbances, see Equations (87) and (89), it is permitted to intersect disturbance estimates at a specific time step k with those of previous steps and to employ those intersected results as virtual measurements in analogy to the sensor models (90).

In the system model (89), the influence of the yaw rate $x_3$ on the system matrix during the state prediction is accounted for by an interval parameter that corresponds to its estimated range from the preceding correction step.

Without any modification, the proposed method is capable of handling imprecisely known damping coefficients $c\in \left[c;\overline{c}\right]$ and $c_{\mathrm{r}}\in \left[c_{\mathrm{r}};{\overline{c}}_{\mathrm{r}}\right]$ or uncertain masses and mass moments of inertia.

5.2. Simulation Results

To illustrate the performance of the thick ellipsoid state estimation algorithm of Section 4, two different levels of uncertainty are compared for the sensor models (90). These are

$${\delta }_i^2=0.005\phantom{1.}$$

(91)

and

$${\delta }_i^2=0.0001.$$

(92)

with the help of these parameters, bounded measurement noise (in the sense of scaled uniformly distributed random numbers) has been generated for the simulations presented in this section.

Moreover, we assume true initial system states ${\mathbf{x}}_0$ randomly chosen from the interior of an ellipsoid centered at

$${\mathit{\mu }}_0^{*}={\left[\begin{array}{cccccc}2& 0.1& 0.01& 0.1& 0& 0\end{array}\right]}^T$$

(93)

with the diagonal shape matrix

$${\mathbf{Q}}_0^{*}=\mathrm{diag}\left\{\left[\begin{array}{cccccc}1& 1& 0.1& 1& 1& 1\end{array}\right]\right\}.$$

(94)

The hovercraft’s external disturbances originate from setting $d_1=-0.25$ and $d_3=0.05$.

The state estimation procedure is initialized with an ellipsoid (containing the true initial state) with the same shape matrix (94) and the shifted center ${\mathit{\mu }}_0={\mathit{\mu }}_0^{*}-{\left[\begin{array}{cccccc}0& 0& 0& 0.1& 0& 0\end{array}\right]}^T$.

Figure 10 and Figure 11 display the evolution of the outer state enclosures for $k=1000$ steps of the estimation algorithm from Figure 6 in the left columns. In the right columns of these figures, an enlarged view is given which compares outer and inner bounds of the thick ellipsoid after the first prediction step, the propagation of the prediction results in terms of outer bounds and the respective corrected ellipsoids after accounting for the measured data. In addition, the measured and true state values are depicted.

Throughout this section, the intersection of the predicted ellipsoids with the degenerate information according to (90) is performed by the application of formula (53). The gap between the outer and inner enclosures for $k=1$ indicates the non-negligible influence of nonlinearities for sufficiently uncertain values of $x_{3,k}$. For the sake of forecasting the domains of reachable states, only the outer bounds are investigated further in Figure 12 and Figure 13. It can be seen that the smaller measurement uncertainty (92) leads to reduced uncertainty in the state reconstruction. Most noticeable, however, is the fact that the observer becomes also applicable to reconstruct bounds for the disturbance $d_{3,k}$ in this case which are tighter than their initialization and consistent with the system model as well as the uncertain state observations, cf. Figure 13d.

The disturbance estimation was performed as described at the end of the previous subsection by intersecting the current estimate for the disturbance (after the correction step) with the previously obtained best disturbance bounds.

For the case of temporally varying disturbances, this procedure has to be modified. Instead of intersecting the bounds from previous time steps, it would be necessary to account for worst-case variations of the disturbance during the prediction step by inflating the corresponding domains similarly to the influence of process noise in a Kalman filter implementation.

For the time-invariant case, future work can also aim at disturbance reconstruction by accounting for the estimates over a window of multiple time steps k and to find a common solution for the resulting set of algebraic constraints which are obtained when evaluating the state Equations (89) in which both $x_{i,k}$ and $x_{i,k+1}$ are replaced by the outcomes of the respective correction steps.

6. Conclusions and Outlook on Future Work

In this paper, a thick ellipsoid state estimation procedure was presented which allows for a computation of outer and (if the uncertainties are not too large) inner ellipsoidal enclosures. The advantage of the presented technique over other ellipsoidal bounding approaches such as those in [19,56] is the applicability to general (differentiable) nonlinear state equations (in this paper, in quasi-linear form) as well as its reduced complexity by using a solution technique in which complex minmax optimization procedures are replaced by simple one-parameter optimizations which only require the online computation of eigenvalues.

In future work, a combination of stochastic and set-valued uncertainty models can be investigated using the novel thick ellipsoid definitions. Other existing techniques such as those presented in [61] also make use of a Kalman filter-like structure. However, they are mainly focused on linear measurement equations. This can be circumvented by the proposed technique if the system’s output equations are reformulated in a quasi-linear form and the influence of nonlinearities is mapped onto the uncertainty that represents the measurement noise.

In addition, the simulation results in Section 5 have shown that the observer approach is capable of reconstructing bounded disturbances. This property can be exploited in future work also for an observer-based parameter identification. Finally, the thick ellipsoid-based state prediction can be employed directly for a numerical evaluation of classical Luenberger-like observers and Extended Kalman Filters. With these results, it will become possible to perform a stability analysis of these estimation schemes when applied to nonlinear dynamic systems. The procedure will be similar to the one proposed in [33], where a centered-form representation of discrete-time dynamic systems was employed for the derivation of stability contractors.