An Ellipsoidal Predictor–Corrector State Estimation Scheme for Linear Continuous-Time Systems With Bounded Parameters and Bounded Measurement Errors

For linear time-invariant dynamic systems with exactly known coefficients of their system matrices for which measurements with bounded errors are available at discrete time instants, an optimal polygonal state estimation scheme was recently published. This scheme allows for tightly enclosing all possible state trajectories in presence of uncertain, but bounded, system inputs which may be varying arbitrarily within in their bounds. Moreover, this approach is also capable of accounting for uncertainty related to the measurement time instants. However, the drawback of this polygonal technique is its rapidly increasing complexity for larger system dimensions. For that reason, the polygonal state enclosures are replaced by a computationally less expensive, but nearly optimal, ellipsoidal enclosure technique in this paper. Numerical simulations for representative benchmark examples focusing both on applications with precisely known and uncertain parameters conclude this contribution.

1 Introduction

Linear continuous-time system models with uncertain, but temporally constant parameters in their system matrices as well as bounded, arbitrarily varying external inputs can be used for the mathematical representation of a large number of system models in the frame of control design and state estimation. These system models belong to the broad class of differential inclusions for which bounds on the temporal variation rates of state variables are formulated in terms of inequality constraints or interval variables. For an overview of the methodological challenges related to this class of dynamic system models, the reader is referred to the works of Brogliato and Tanwani (2020), Stewart (2011), Filippov (1988), and the references therein.

Possible applications can be found in distributed heating systems (Rauh et al., 2015), where external inputs related to heat conduction and radiation might have temporally varying or unknown state-dependent characteristics, modeling drive trains with elastic shafts (Amann et al., 2004) characterized by uncertain load and disturbance torques, mechanical positioning systems which are included, for example, in piezo servo hydraulic actuation techniques for camless combustion engines (Haus et al., 2014), current and torque control for permanent magnet synchronous machines in a d-q-coordinate system (Mousavi et al., 2020), or oscillation attenuation and trajectory tracking in robotics applications (Rauh et al., 2013). In all of these applications, it is typically desired to reconstruct internal states on the basis of measured data with uncertainty and to make guaranteed statements about the applicability and safety of the resulting state trajectories in terms of a computation of outer enclosures for those states that are reachable at a certain point of time. Examples are soft landing capabilities for controlled valves in the aforementioned combustion engine (Di Bernardo et al., 2012), the guaranteed compliance of controlled electric motors with hard current and torque constraints or the guaranteed collision-free trajectory control of autonomous robots (Zhang et al., 2019).

The uncertain and bounded system inputs then contain (depending on the actual application as already partially mentioned above) either control inputs which are not perfectly known due to actuator imprecision or the influence of external disturbances.

It should be pointed out that state estimation for linear dynamic systems in a point-valued and stochastic context is a well developed area. The classically applied estimation schemes are based on Luenberger observers and Kalman filters (Kalman, 1960; Luenberger, 1964; Stengel, 1994; Davis, 2002; Anderson and Moore, 2005). However, these techniques either provide purely point-valued state estimates in the case of deterministic approaches (Luenberger-like observers) or aim at a reconstruction of an expected value of a probability distribution coming along with confidence bounds in the case of the Kalman filter. Unfortunately, these techniques cannot provide strictly guaranteed outer bounds on the domains of reachable states if models described by differential inclusions are considered, where information about the probability density functions of the uncertain but bounded variables is totally absent.

Moreover, also classical simulation approaches making use of set-valued uncertainty representations are not applicable in the context of differential inclusions. Approaches such as AWA (Lohner, 1987; Lohner, 1988), VNODE and VNODE-LP (Nedialkov, 1999; Nedialkov, 2007; Nedialkov, 2011), VSPODE (Lin and Stadtherr, 2007), COSY VI (Berz and Makino, 1998; Hoefkens, 2001) and RioT (Eble, 2007) are all based on a temporal series expansion of the solution to the dynamic system model under consideration. Due to the fact that the uncertain system inputs are assumed to be bounded in this paper, however, without any knowledge on their variation rates, these techniques are not applicable because they heavily rely on high-order temporal series expansions. The only set-valued approach that is widely known to work in such contexts is the pure computation of bounds on the basis of a Picard iteration (Deville et al., 2002). This approach, however, as a set-valued generalization of the explicit Euler method, inevitably leads to a blow-up of the computed state enclosures which makes it only applicable for very short prediction horizons. This problem can be circumvented partially by defining exponential state enclosures as presented in Rauh et al. (2016). These enclosures, however, suffer from a non-negligible overestimation due to the wrapping effect of interval analysis (Jaulin et al., 2001b; Lohner, 2001) if the state equations under consideration are strongly coupled. This overestimation can be countered by using the polygonal solution representations detailed below or by exploiting specific monotonicity properties of the system model (so-called cooperativity, as it is typically done in the frame of continuous- and discrete-time interval observer design (Raïssi et al., 2012; Efimov et al., 2013; Smith, 1995)¹). The property of cooperativity only holds for a quite restricted class of system models (e.g., thermo-fluidic ones) after a first-principle modeling. So, simulation and state estimation techniques that rely on cooperativity often need to be implemented in combination with an additional change of coordinates (Marouani et al., 2021; Rauh and Kersten, 2021). This change of coordinates is on purpose avoided in the current paper due to the following reasons: Especially for uncertain systems, such changes of coordinates are non-trial to find (if they exist at all) and they typically lead to overestimation that can be avoided by the use of the ellipsoidal enclosure technique derived in this paper.

Note, also a pure replacement of continuous-time system models by using explicit or implicit Euler schemes, Runge–Kutta methods or other alternatives (Hairer et al., 2000) does not solve the issues related with differential inclusions, even if measured data were available at some discrete instants of time in the context of state estimation. The reason for this are the arising temporal discretization errors that may lead to approximations that do not cover the true system dynamics. Thinking of a set-valued context—which represents the reachable state domains by intervals, polygons, ellipsoids, etc.—a naive discretization then provides sets that often do not include all of the actually reachable states during the transition from one measurement instant to the next so that the property of guaranteed state enclosures would be invalidated.

In previous work (Rohou and Jaulin, 2021) of the authors, a polygonal approach was presented that allows for enclosing the domains of reachable states for this class of linear continuous-time system models, representing differential inclusions, in an optimal way. However, the drawback of the technique is the large computational complexity for increasing system dimensions. To circumvent this difficulty, an alternative to the polygonal technique is presented in the current paper which makes use of a computationally less demanding ellipsoidal state enclosure technique. Besides its reduced computational complexity, which makes it applicable to systems of higher dimensions, this approach is also directly applicable to linear time-invariant systems with uncertain parameters and provides enclosures for the domains of reachable states which are still nearly optimal with respect to the widths of the enclosures.

For the purpose of a guaranteed state estimation, a guaranteed evaluation of continuous-time system models as a kind of state prediction algorithm is interfaced with a state correction scheme that is executed at those discrete time instants at which measurements become available.

In contrast to the aforementioned polygonal approach, a small amount of overestimation is introduced by the new ellipsoidal approach due to the fact that the mapping of ellipsoidal domains over continuous-time state equations with uncertain parameters does not result in solution sets that are exactly ellipsoidal. The same also holds for the correction step based on the intersection of ellipsoids with measured data. As shown by means of numerical simulations, the state enclosures remain reasonably tight for sufficiently tight interval parameters included in the system matrices and for sufficiently accurate knowledge of the external system inputs.

This paper is structured as follows. Section 2 provides a review of fundamentals of set-membership state estimation procedures. On this basis, Sections 3, 4 present details of the optimal polygonal and the proposed novel ellipsoidal state estimation scheme. A numerical benchmark example, aiming at a comparison of these enclosure techniques is discussed in Section 5. Besides considering linear time-invariant systems with exactly known parameters, it also visualizes the novel ellipsoidal predictor–corrector state estimator for an example containing bounded uncertainty in the system matrix. Finally, an outlook on future work summarizes this paper in Section 6.

2 Set-Membership State Estimation Approach

The approach provided in Rohou and Jaulin (2021) considers linear time-invariant dynamical systems that are modeled in the form

$\dot{\mathbf{x}}=\mathbf{A}\left(\mathbf{p}\right)\cdot \mathbf{x}+\mathbf{B}\left(\mathbf{p}\right)\cdot \mathbf{u},(1)$

where $\mathbf{A}\left(\mathbf{p}\right)$ is the system matrix and $\mathbf{B}\left(\mathbf{p}\right)$ the corresponding input matrix. Both $\mathbf{A}\left(\mathbf{p}\right)$ and $\mathbf{B}\left(\mathbf{p}\right)$ may depend on uncertain parameters $\mathbf{p}\in {\mathbb{R}}^{n_{\text{p}}}$, which are assumed to be time-invariant and bounded by the interval $\mathbf{p}\in \left[\mathbf{p}\right]$, where the element-wise relation ${\underline{p}}_i\le p_i\le {\bar{p}}_i$ holds for each vector component i ∈ {1, … , n_p}. In addition to the dynamic system model above, we assume that direct measurements of the state vector $\mathbf{x}\in {\mathbb{R}}^{n_{\text{x}}}$ may be available and bounded by a box $[{\mathbf{y}}_k]\in {\mathbb{IR}}^{n_{\text{x}}}$ at the time instant t_k, where ${\mathbb{IR}}^{n_{\text{x}}}$ denotes the set of axis-aligned interval boxes in ${\mathbb{R}}^{n_{\text{x}}}$, cf. Jaulin et al. (2001b); Mayer (2017); Moore et al. (2009).

2.1 State Tube

A tube $\mathbb{X}(\cdot )$² (Kurzhanski and Filippova (1993)) is used to continuously enclose the feasible states of a dynamic system. For that purpose, $\mathbb{X}(\cdot )$ is defined as a set of trajectories $\mathbb{X}(\cdot ):\mathbb{R}\mapsto \mathcal{P}({\mathbb{R}}^{n_{\text{x}}})$ which satisfies

$\mathbf{x}\left(\cdot \right)\in \mathbb{X}\left(\cdot \right)\iff \forall t,\mathbf{x}\left(t\right)\in \mathbb{X}\left(t\right).(2)$

Assume that for $t\in \left[0;\bar{t}\right]$, the state trajectory x(⋅) is known to be inside a prior tube $\check{\mathbb{X}}(\cdot )$. This prior tube is initialized with the direct measurements of the states given by intervals [y_k]. It may also be unbounded at the time instants for which no state information is available a priori. We want to compute recursively the tightest tube $\mathbb{X}(\cdot )$ for x(⋅) that is consistent with both the prior tube $\check{\mathbb{X}}(\cdot )$ (including the measurements) and the state Eq. 1. This sequence can be interpreted as an extension of the state estimator proposed in Bertsekas and Rhodes (1971) or Jaulin et al. (2001a) to continuous-time systems.

For a given initial vector x₁ defined at t₁, the state at time t₂ is expressed by

${\mathbf{x}}_2={\mathbf{\Phi }}_{t_1,t_2}^{\mathbf{u}\left(\cdot \right)}\left({\mathbf{x}}_1\right),(3)$

where ${\mathbf{\Phi }}_{t_1,t_2}$ is the flow of the dynamical system according to its input u(⋅). The flow can be extended to sets (Aubin and Frankowska (1990)) as follows³:

${\mathbf{\Phi }}_{t_1,t_2}^{\mathbf{u}\left(\cdot \right)}\left({\mathbb{X}}_1\right)=\left\{{\mathbf{x}}_2|\exists {\mathbf{x}}_1\in {\mathbb{X}}_1,{\mathbf{x}}_2={\mathbf{\Phi }}_{t_1,t_2}^{\mathbf{u}\left(\cdot \right)}\left({\mathbf{x}}_1\right)\right\}.(4)$

Furthermore, when the input u(⋅) is uncertain but known to be inside a tube $\mathbb{U}(\cdot )$, the set flow given by Eq. 4 becomes

${\mathbf{\Phi }}_{t_1,t_2}^{\mathbb{U}\left(\cdot \right)}\left({\mathbb{X}}_1\right)=\bigcup _{\mathbf{u}\left(\cdot \right)\in \mathbb{U}\left(\cdot \right)}{\mathbf{\Phi }}_{t_1,t_2}^{\mathbf{u}\left(\cdot \right)}\left({\mathbb{X}}_1\right).(5)$

We now define the posterior state tube as the tightest tube $\widehat{\mathbb{X}}(\cdot )$ for x(⋅) consistent with the prior tube $\check{\mathbb{X}}(\cdot )$, the feasible inputs u(⋅) in $\mathbb{U}(\cdot )$, and the dynamic system model. We obtain

${\widehat{\mathbb{X}}}_t=\bigcup _{\mathbf{u}\left(\cdot \right)\in \mathbb{U}\left(\cdot \right)}\bigcap _{\tau \in \left[0;\bar{t}\right]}{\mathbf{\Phi }}_{\tau ,t}^{\mathbf{u}\left(\cdot \right)}\left({\check{\mathbb{X}}}_{\tau }\right)(6)$

as illustrated in Figure 1.

FIGURE 1

FIGURE 1. Illustration of one slice of a posterior tube $\widehat{\mathbb{X}}(\cdot )$ at time t₂ (red hatched part). This slice is defined as the intersection set ${\widehat{\mathbb{X}}}_{t_2}={\mathbf{\Phi }}_{t_1,t_2}^{\mathbb{U}(\cdot )}({\check{\mathbb{X}}}_{t_1})\cap {\mathbf{\Phi }}_{t_3,t_2}^{\mathbb{U}(\cdot )}({\check{\mathbb{X}}}_{t_3})$, with the prior states ${\check{\mathbb{X}}}_{t_1}$ and ${\check{\mathbb{X}}}_{t_3}$. The flow function ${\mathbf{\Phi }}_{t,t_2}^{\mathbb{U}(\cdot )}$ transports the prior sets from t to t₂, according to the feasible inputs $\mathbb{U}(\cdot )$ over $\left[t;t_2\right]$. The result of each transportation, from t₁ and from t₃, is depicted in blue at the time instant t₂.

The state estimation consists in the approximation of $\widehat{\mathbb{X}}(\cdot )$, $\forall t\in \left[0;\bar{t}\right]$. A numerical method can be used to this end that is implemented by means of a recursive algorithm.

Remark 1. For system models (1), in which the parameters p are uncertain and temporally varying, we assume in this paper that the system matrix A only contains the time-invariant parts of the dynamics, while all further time-varying effects are assumed to be included in the additive offset term $\mathbf{B}\left(\mathbf{p}\right)\cdot \mathbf{u}$ by means of a suitable reformulation.

2.2 Recursive Procedure

We now define the sets ${\mathbb{X}\to \overrightarrow{\mathbb{X}}}_t$ and ${\stackrel{⃖}{\mathbb{X}}}_t$ that correspond respectively to the sets of all x(t) consistent with the past (before the point of time t) and the future (after t). The following theorem allows us to implement the exact sequence to compute the tube $\widehat{\mathbb{X}}(\cdot )$. This theorem as well as Theorem 2 are published in Rohou and Jaulin (2021) together with a detailed proof.

Theorem 1. Given the sampling times of measurements

$\mathbb{T}=\left\{0,\delta ,2\delta ,\dots ,\bar{k}\delta \right\}=\left\{t_0,t_1,t_2,\dots ,t_{\bar{k}}\right\},(7)$

a prior tube $\check{\mathbb{X}}(\cdot )$ containing the state trajectory x(⋅), and a tube $\mathbb{U}(\cdot )$ for the inputs. The posterior tube $\widehat{\mathbb{X}}(\cdot )$ can be defined recursively by

$\begin{array}{c}\hfill {\mathbb{X}}_{t_k}={\mathbf{\Phi }}_{t_{k-1},t_k}^{\mathbb{U}\left(\cdot \right)}\left({\mathbb{X}}_{t_{k-1}}\right)\cap \bigcup _{\mathbf{u}\left(\cdot \right)\in \mathbb{U}\left(\cdot \right)}\bigcap _{\tau \in \left[t_{k-1};t_k\right]}{\mathbf{\Phi }}_{\tau ,t_k}^{\mathbf{u}\left(\cdot \right)}\left({\check{\mathbb{X}}}_{\tau }\right)\\ \hfill {\stackrel{⃖}{\mathbb{X}}}_{t_k}={\mathbf{\Phi }}_{t_{k+1},t_k}^{\mathbb{U}\left(\cdot \right)}\left({\stackrel{⃖}{\mathbb{X}}}_{t_{k+1}}\right)\cap \bigcup _{\mathbf{u}\left(\cdot \right)\in \mathbb{U}\left(\cdot \right)}\bigcap _{\tau \in \left[t_k;t_{k+1}\right]}{\mathbf{\Phi }}_{\tau ,t_k}^{\mathbf{u}\left(\cdot \right)}\left({\check{\mathbb{X}}}_{\tau }\right)\\ {\widehat{\mathbb{X}}}_{t_k}={\mathbb{X}}_{t_k}\cap {\stackrel{⃖}{\mathbb{X}}}_{t_k}\hfill \end{array}(8)$

with ${\mathbb{X}}_{t_0}=\check{\mathbb{X}}(t_0)$ and ${\stackrel{⃖}{\mathbb{X}}}_{t_{\bar{k}}}=\check{\mathbb{X}}(t_{\bar{k}})$. Moreover, for $t\in \left[t_k;t_{k+1}\right]\backslash \mathbb{T}$, we have

${\widehat{\mathbb{X}}}_t={\mathbf{\Phi }}_{t_k,t}^{\mathbb{U}\left(\cdot \right)}\left({\widehat{\mathbb{X}}}_{t_k}\right).(9)$

2.3 State Estimator in the Linear Case

The exact sequence suggested by Theorem 1 is valid regardless whether the system is linear or nonlinear. Nevertheless, it can be implemented exactly on a computer only in the linear case. This is due to the fact that we generally do not have an explicit expression for the flow when dealing with nonlinear systems described by $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u})$. This section shows how an accurate representation of the exact sequence for state estimation that is given by Eq. 8 can be implemented in the linear case.

2.3.1 Principle

For a linear dynamic system (1), a closed-form expression for the flow is given by

${\mathbf{\Phi }}_{t_1,t_2}^{\mathbf{u}\left(\cdot \right)}\left(\mathbf{x}\right)=e^{\mathbf{A}\left(\mathbf{p}\right)\cdot \left(t_2-t_1\right)}\mathbf{x}+{\int }_{t_1}^{t_2}{e}^{\mathbf{A}\left(\mathbf{p}\right)\cdot \left(t_2-\tau \right)}\cdot \mathbf{B}\left(\mathbf{p}\right)\cdot \mathbf{u}\left(\tau \right)d\tau .(10)$

To use Eq. 10 in the case where interval uncertainties exist for the point of time t, we need to introduce the concept of the exponential for interval matrices [A] (in our case, A([p])) which has to be understood in terms of a set-theoretic meaning (Goldsztejn and Neumaier (2014)). It is defined as the tightest enclosure in the form of an interval matrix which contains all feasible exponentials of A(p), assuming that p ∈ [p], i.e.,

$e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)}=\left[\left\{\mathbf{M}|\exists \mathbf{p}\in \left[\mathbf{p}\right],\mathbf{M}=e^{\mathbf{A}\left(\mathbf{p}\right)}\right\}\right].(11)$

Moreover, given an interval matrix [A] and a set of vectors $\mathbb{X}$, we define the product $[\mathbf{A}]\cdot \mathbb{X}$ as the tightest box enclosing all feasible products A⋅x assuming that $\mathbf{x}\in \mathbb{X}$ and A ∈ [A], i.e.,

$\left[\mathbf{A}\right]\cdot \mathbb{X}=\left[\left\{\mathbf{y}|\exists \mathbf{A}\in \left[\mathbf{A}\right],\exists \mathbf{x}\in \mathbb{X},\mathbf{y}=\mathbf{A}\cdot \mathbf{x}\right\}\right].(12)$

The following theorem corresponds to an implementation of the exact sequence provided by Theorem 1. For a proof, the reader is referred to Rohou and Jaulin (2021).

Theorem 2. Given the sampling times of measurements

$\mathbb{T}=\left\{0,\delta ,2\delta ,\dots ,\bar{k}\delta \right\}=\left\{t_0,t_1,t_2,\dots ,t_{\bar{k}}\right\},(13)$

a prior tube $\check{\mathbb{X}}(\cdot )$ containing the state trajectory x(⋅), and a piecewise constant tube

$\mathbb{U}\left(\cdot \right)=\left\{\mathbf{u}\left(\cdot \right)|\forall k,\forall t\in \left[k\delta ;\left(k+1\right)\delta \right],\mathbf{u}\left(t\right)\in {\left[\mathbf{u}\right]}_k\right\}(14)$

containing u(⋅), where ${[\mathbf{u}]}_k,k\in \{0,\dots ,\bar{k}-1\}$ is a boxed slice of the tube $\mathbb{U}(\cdot )$ as illustrated by Figure 2. We have

$\begin{array}{c}{\mathbb{X}}_{t_k}\subset {\check{\mathbb{X}}}_{t_k}\cap \left\{e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \delta }\cdot {\mathbb{X}}_{t_{k-1}}+\delta e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]}\cdot \mathbf{B}\left(\left[\mathbf{p}\right]\right)\cdot {\left[\mathbf{u}\right]}_{k-1}\right\}\hfill \\ {\stackrel{⃖}{\mathbb{X}}}_{t_k}\subset {\check{\mathbb{X}}}_{t_k}\cap \left\{e^{-\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \delta }\cdot {\stackrel{⃖}{\mathbb{X}}}_{t_{k+1}}-\delta e^{-\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]}\cdot \mathbf{B}\left(\left[\mathbf{p}\right]\right)\cdot {\left[\mathbf{u}\right]}_k\right\}\hfill \\ {\widehat{\mathbb{X}}}_{t_k}={\mathbb{X}}_{t_k}\cap {\stackrel{⃖}{\mathbb{X}}}_{t_k}\hfill \end{array}(15)$

with ${\mathbb{X}}_{t_0}={\check{\mathbb{X}}}_{t_0}$ and ${\stackrel{⃖}{\mathbb{X}}}_{t_{\bar{k}}}={\check{\mathbb{X}}}_{t_{\bar{k}}}$.Moreover, when t is not consistent with the tube discretization, i.e., for $t\in \left[t_k;t_{k+1}\right]\backslash \mathbb{T}$, we have

${\widehat{\mathbb{X}}}_t=\left\{e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left(t-t_k\right)}\cdot {\vec{\mathbb{X}}}_{t_k}+\left(t-t_k\right)e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;t-t_k\right]}\cdot \mathbf{B}\left(\left[\mathbf{p}\right]\right)\cdot {\left[\mathbf{u}\right]}_k\right\}.(16)$

FIGURE 2

FIGURE 2. Illustration of a tube $\mathbb{U}(\cdot )=[\mathrm{u}](\cdot )$ made of boxed slices. This representation can be used in order to reliably enclose sets of trajectories. For instance, the thick blue box depicts the [u]_k=2 slice containing all feasible values for u(t), $t\in \left[t_2;t_3\right]$, cf. Rohou and Jaulin (2021).

2.4 Wrappers for Enclosing ${\mathbb{X}}_t$

Now that we have a reliable definition for the enclosure $\mathbb{X}(\cdot )$ from a bounded input $\mathbb{U}(\cdot )$ (without any information about its temporal variation rates within the given bounds), that can be numerically represented by a boxed tube [u](⋅) and guaranteed to enclose u(⋅), it remains to reliably implement the sets $\mathbb{X}(\cdot )$. Several wrappers can be considered, such as intervals, polygons, or ellipsoids. The use of intervals amounts to representing ${\mathbb{X}}_t$ as a box (which is the implementation of tubes described, for instance, in Rohou et al. (2017), and used for [u](⋅)). The following section provides the polygon-based algorithm provided in Rohou and Jaulin (2021) that has been shown to be optimal. The further sections of this paper then introduce a new method using ellipsoids as wrappers. In the last part of the paper, we will compare all three methods to assess the pros and cons of each of these alternatives.

3 Optimal Polygonal Method

3.1 Polygonal Sequences

Consider a polygon ${\mathcal{P}}_k$ which contains ${\check{\mathbb{X}}}_{t_k}$. We apply the following polygonal sequence as a specific formulation of the general approach described in Eq. 15:

• First, forward in time propagation, for $k\in \{1,\dots ,\bar{k}\}$:

${\mathcal{P}}_k≔{\mathcal{P}}_k\cap \left\{e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \delta }\cdot {\mathcal{P}}_{k-1}+\delta e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]}\cdot \mathbf{B}\left(\left[\mathbf{p}\right]\right)\cdot {\left[\mathbf{u}\right]}_{k-1}\right\},(17)$

• then, backward in time propagation, for $k\in \{\bar{k}-1,\dots ,0\}$:

$\begin{array}{c}\hfill {\mathcal{P}}_k≔{\mathcal{P}}_k\cap \left\{e^{-\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \delta }\cdot {\mathcal{P}}_{k+1}-\delta e^{-\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]}\cdot \mathbf{B}\left(\left[\mathbf{p}\right]\right)\cdot {\left[\mathbf{u}\right]}_k\right\}\hfill \end{array},(18)$

• finally, polygonal enclosure between different sampling points, $k\in \{0,\dots ,\bar{k}-1\}$:

${\mathcal{P}}_{k:k+1}=\left\{e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]}\cdot {\mathcal{P}}_k+\left[0,\delta \right]e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]}\cdot \mathbf{B}\left(\left[\mathbf{p}\right]\right)\cdot {\left[\mathbf{u}\right]}_k\right\}.(19)$

From Theorem 2, we know that for all $t_k\in \mathbb{T}$,

$\left\{\begin{array}{c}{\widehat{\mathbb{X}}}_{t_k}\subset {\mathcal{P}}_k,\hfill \\ \forall t\in \left[t_k;t_{k+1}\right],{\widehat{\mathbb{X}}}_t\subset {\mathcal{P}}_{k:k+1}.\hfill \end{array}\right.(20)$

The implementation of this procedure resembles the exact sequence given by Theorem 2. As a consequence, it is not necessary in practice to apply the polygonal sequence several times in a recursive manner to obtain an accurate enclosure.

In the following section, the abbreviations

$\mathbf{\Phi }\left(\mathbf{p}\right)=e^{\mathbf{A}\left(\mathbf{p}\right)\cdot \delta }\text{and}\mathbf{\Psi }\left(\mathbf{p}\right)=\delta \cdot e^{\mathbf{A}\left(\mathbf{p}\right)\cdot \delta }\cdot \mathbf{B}\left(\mathbf{p}\right)(21)$

are used to denote the exact temporal discretization of a linear time-invariant system used repeatedly in expressions such as Eq. 17, where suitable extensions will be defined for the case of interval parameters and uncertainty in the sampling (respectively, measurement) time instants.

Remark 2. In contrast to the numerical, approximating discretization schemes mentioned in the introduction of this paper, this kind of discretization is exact and does not suffer from temporal truncation errors if all matrix exponentials involved are represented in terms of interval bounds, cf. Eq. 11.

Remark 3. Due to the fact that the proposed estimation scheme involves temporal forward and backward evaluations of the solution tubes, one of the evaluation directions commonly deals with stable and the other one with unstable dynamics. Hence, the examples included in this paper directly show the applicability of the proposed method to both stable and unstable systems, where the evaluation in the stable direction may (depending on the bounds of uncertain parameters and inputs) show a contraction of the computed state enclosures and the opposite direction of evaluation a corresponding widening caused by the unstable dynamics.

3.2 Pros and Cons of the Polygonal Approach

Set-membership methods are known to be guaranteed because the propagation of uncertainties over time is made rigorously by using reliable operations such as interval analysis (Jaulin et al., 2001b; Mayer, 2017). In the case of the polygonal approach proposed in the previous work of Rohou and Jaulin (2021), the method is also optimal since it is a direct extension of Theorem 2 without wrapping effects. Such wrapping effects would result from enclosing complex shaped sets by more simple outer bounds (such as boxes) which are subsequently propagated further (in time) and then again enclosed by simplified outer bounds. Indeed, each operation during the analysis of linear time-invariant dynamic systems involving convex polygons as state enclosures results in another convex polygon, unfortunately, commonly with an increasing number of vertices required to represent the exact solution domains. This mapping of bounds to the same class of set representation, however, is not the case with boxes, which inevitably leads to a strong pessimism in recursive algorithms (Krasnochtanova et al., 2010).

As a consequence, if we assume that the sampling time δ is infinitely small and that the computer works exactly with real numbers instead of floating point numbers, the polygonal approach can be qualified as exact since it does not introduce any pessimism and does not lose any feasible values.

Despite these advantages, the polygonal approach also comes with two drawbacks:

• Set operations on polygons lead to an increasing number of vertices, because the output of a vertex by an inclusion function is a box. The computational complexity can be reduced by removing vertices in some reliable way. This is done in the work of Rohou and Jaulin (2021) but at the expense of significant computation times. The optimality of the results is also lost during the simplification procedure of these polygons.

• The reliable implementation of polygons and their operations are complex and limited to problems of low dimensions. This is the reason why the implementation proposed in the previous work is limited to sets $\mathbb{X}(\cdot )$ in ${\mathbb{R}}^2$, which can be seen as a strong limitation for a variety of real-life applications, where set-valued state estimation algorithms need to be implemented on embedded systems in a real-time capable manner.

These reasons motivate the study of ellipsoids as an alternative wrapper for $\mathbb{X}(\cdot )$ sets in the following section.

4 Ellipsoidal State Estimation Procedure

To derive the alternative ellipsoidal wrapper and corresponding state estimation procedure, assume first that a discrete-time system model in the form

${\mathbf{x}}_{k+1}=\mathbf{\Phi }\left(\mathbf{p}\right)\cdot {\mathbf{x}}_k+\mathbf{\Psi }\left(\mathbf{p}\right)\cdot {\mathbf{u}}_k(22)$

is given, where $\mathbf{\Phi }\left(\mathbf{p}\right)$ is the system matrix and $\mathbf{\Psi }\left(\mathbf{p}\right)$ the corresponding input matrix, see also Eq. 21 and Remark 1. Both $\mathbf{\Phi }\left(\mathbf{p}\right)$ and $\mathbf{\Psi }\left(\mathbf{p}\right)$ may depend on uncertain parameters $\mathbf{p}\in {\mathbb{R}}^{n_{\text{p}}}$ as it was already discussed for the continuous-time system model (Eq. 1). The following state estimation procedure is a generalization of the thick ellipsoid state estimation presented in Rauh et al. (2021) in the sense that it allows for taking into consideration the additive term $\mathbf{\Psi }\left(\mathbf{p}\right)\cdot {\mathbf{u}}_k$ in the state Eq. 22. This term is either described in terms of an interval enclosure that is subsequently bounded by a guaranteed outer ellipsoid or by the vertices of a polygon as motivated by the previous section. Also in the latter case, this term will be enclosed by a tight (Löwner-John-type, cf. John (1948)) ellipsoid bound to obtain a simple to evaluate state estimation procedure.

In addition to the dynamic system model above, we assume that measurements y_m,k+1 with an ellipsoidal uncertainty model

${\left(\mathbf{C}\cdot {\mathbf{x}}_{k+1}-{\mathbf{y}}_{\text{m},k+1}\right)}^T\cdot {\mathbf{Q}}_{\text{m}}^{-1}\cdot \left(\mathbf{C}\cdot {\mathbf{x}}_{k+1}-{\mathbf{y}}_{\text{m},k+1}\right)\le 1(23)$

are available. To make the ellipsoidal set intersection operator introduced in Rauh et al. (2021) applicable to perform the state estimation tasks, we assume that the inequality (Eq. 23) can be reformulated equivalently in terms of the constraint

${\left({\mathbf{x}}_{k+1}-{\mathbf{y}}_{\text{m},k+1}^{\prime }\right)}^T\cdot {\mathbf{P}}_{\text{m}}^{\prime }\cdot \left({\mathbf{x}}_{k+1}-{\mathbf{y}}_{\text{m},k+1}^{\prime }\right)\le 1.(24)$

In this inequality (Eq. 24), the matrix ${\mathbf{P}}_{\text{m}}^{\prime }$ is a purely diagonal matrix if the system’s output matrix C in Eq. 23 represents a direct measurement of selected components of the state vector x_k+1, i.e., being an all-zero matrix except for a single entry with the value one per row. If less outputs than state variables are available as measured data, the matrix ${\mathbf{P}}_{\text{m}}^{\prime }$ is not invertible and represents a degenerate ellipsoid in the n_x-dimensional state space.

Using this notation, entries in the vector ${\mathbf{y}}_{\text{m},k+1}^{\prime }$ that correspond to non-measured components of the state vector at the time instant k + 1 are set to the associated ellipsoid midpoint μ_k+1 obtained from the state prediction, while all other components correspond to the point-valued measured data y_m,k+1.

In general, ellipsoidal state bounds at the time instant k are denoted as

${\mathcal{E}}_k\left({\mathit{\mu }}_k,{\mathbf{\Gamma }}_k\right)≔\left\{{\mathbf{x}}_k\in {\mathbb{R}}^{n_{\text{x}}}|{\left({\mathbf{x}}_k-{\mathit{\mu }}_k\right)}^T{\mathbf{\Gamma }}_k^{-T}{\mathbf{\Gamma }}_k^{-1}\left({\mathbf{x}}_k-{\mathit{\mu }}_k\right)\le 1\right\}(25)$

with the positive definite shape matrix ${\mathbf{Q}}_k={\mathbf{\Gamma }}_k{\mathbf{\Gamma }}_k^T\succ 0$ and the midpoint vector ${\mathit{\mu }}_k\in {\mathbb{R}}^{n_{\text{x}}}$.

In the following two subsections, a predictor–corrector type state estimation scheme is described. It consists of firstly propagating ellipsoidal state enclosures from the time step k to the next sampling instant k + 1 according to the system dynamics (Eq. 22). Subsequently, this a-priori information is updated in the set-valued correction step on the basis of the measurement models (Eqs. 23, 24), respectively, describing bounded uncertainty in the sensor readings.

4.1 Ellipsoidal State Prediction Step: Propagation of Outer Ellipsoidal State Bounds

The first three steps of the following state prediction procedure for discrete-time systems (Eq. 22) were basically published and proven in Rauh et al. (2021). The major difference is that the quoted publication makes use of a thick ellipsoid state enclosure approach in which inner and outer ellipsoid bounds for the first term $\mathbf{\Phi }\left(\mathbf{p}\right)\cdot {\mathbf{x}}_k$ in Eq. 22 are determined simultaneously.

However, for the sake of a pure state estimation in the form of computing guaranteed outer enclosures, as it has already been done by using a polygonal approach in the previous section, only the outer bounds of those thick ellipsoids are of interest. Therefore, the state prediction technique described in this subsection is inspired by the approach from Rauh et al. (2021) but leaves out the computation of inner state bounds. However, if inner bounds of the solution set were determined additionally, they could be used to detect scenarios reliably in which the bounded parameters p are the source for significant overestimation. In such cases, the outer and inner bounds of the thick ellipsoids from Rauh et al. (2021), Rauh and Jaulin (2021a), and Rauh and Jaulin (2021b) have a large distance. Then, it is often helpful to subdivide the parameter intervals $\left[\mathbf{p}\right]$ into individual subboxes, to evaluate the prediction individually for each of these boxes, and to compute the ellipsoidal union over them as demonstrated in Section 3.3 of Rauh et al. (2021).

For finding outer ellipsoidal enclosures of the subexpression $\mathbf{\Phi }\left(\mathbf{p}\right)\cdot {\mathbf{x}}_k$ in Eq. 22, we assume that this part of the system model is reformulated in the form

${\mathbf{x}}_{\mathbf{\Phi },k+1}=\mathbf{\Phi }\left(\mathbf{p}\right)\cdot {\mathbf{x}}_k=\mathbf{\Phi }\left(\mathbf{p}\right)\cdot {\mathbf{x}}_k+\tilde{\mathbf{\Phi }}\cdot {\mathit{\mu }}_k+\left(\mathbf{\Phi }\left(\mathbf{p}\right)-\tilde{\mathbf{\Phi }}\right)\cdot {\mathit{\mu }}_k,(26)$

so that an origin-centered ellipsoid can be propagated with the help of the first term in Eq. 26, while the remaining terms account for the influence of the generally non-zero ellipsoid midpoint. For that purpose, introduce the following notation already used in Eq. 26. It is employed during the state prediction phase consisting of the Steps P1–P5:

${\mathcal{E}}_k={\mathcal{E}}_k\left({\mathit{\mu }}_k,{\mathbf{\Gamma }}_k\right)\mathrm{denotes\; the\; ellipsoidal\; uncertainty\; on\; the}\mathrm{non-origin\; centered\; states}{\mathbf{x}}_k,(27)$

${\mathcal{E}}_k={\mathcal{E}}_k\left(\mathbf{0},{\mathbf{\Gamma }}_k\right)\mathrm{the\; uncertainty\; of}{\check{\mathbf{x}}}_k\mathrm{after\; shifting\; the}\mathrm{ellipsoid\; to\; the\; origin,\; and}(28)$

$\tilde{\mathbf{\Phi }}=\mathbf{\Phi }\left(\mathrm{m}\mathrm{i}\mathrm{d}\left(\left[\mathbf{p}\right]\right)\right)\text{the\hspace{0.17em}midpoint\hspace{0.17em}approximation\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}quasi}-\text{linear\hspace{0.17em}system\hspace{0.17em}matrix\hspace{0.17em}with}(29)$

$\mathbf{p}\in \left[\mathbf{p}\right]=\left[\underline{\mathbf{p}};\bar{\mathbf{p}}\right],\text{and\hspace{0.17em}the\hspace{0.17em}interval\hspace{0.17em}midpoint}\mathrm{m}\mathrm{i}\mathrm{d}\left(\left[\mathbf{p}\right]\right)=\frac{1}{2}\cdot \left(\underline{\mathbf{p}}+\bar{\mathbf{p}}\right).(30)$

Step P1: Apply

${\check{\mathbf{x}}}_{k+1}=\mathbf{\Phi }\left(\mathbf{p}\right)\cdot {\check{\mathbf{x}}}_k(31)$

to the ellipsoid ${\mathcal{E}}_k$ in Eq. 28. The shape matrix of the outer ellipsoid enclosure of the image set is described by an ellipsoid with the shape matrix

${\check{\mathbf{Q}}}_{k+1}={\alpha }_{k+1}^2\cdot {\mathbf{\Gamma }}_{k+1}\cdot {\mathbf{\Gamma }}_{k+1}^T,(32)$

where α_k+1 ≥ 0 is the smallest value for which the linear matrix inequality (LMI)

${\mathcal{M}}_{k+1}≔\mathbf{\Lambda }\left[\begin{array}{cc}\hfill -{\mathbf{Q}}_k^{-1}\hfill & \hfill {\mathbf{\Phi }}^T\left(\mathbf{p}\right)\cdot {\tilde{\mathbf{\Phi }}}^{-T}\hfill \\ \hfill {\tilde{\mathbf{\Phi }}}^{-1}\cdot \mathbf{\Phi }\left(\mathbf{p}\right)\hfill & \hfill -{\alpha }_{k+1}^2{\mathcal{R}}_k\hfill \end{array}\right]\mathbf{\Lambda }⪯0,{\mathbf{Q}}_k={\mathbf{\Gamma }}_k\cdot {\mathbf{\Gamma }}_k^T(33)$

is satisfied for all $\mathbf{p}\in \left[\mathbf{p}\right]$ with

${\mathcal{R}}_k≔{\mathbf{\Gamma }}_k\cdot {\mathbf{\Gamma }}_k^T.(34)$

As an addition to the original algorithm in Rauh et al. (2021), the symmetric preconditioning matrix Λ = Λ^T ≻ 0 has been introduced in Eq. 33. It helps to optimize the ellipsoidal enclosures especially in scenarios in which the norms of ${\left({\tilde{\mathbf{\Phi }}}^{-1}\cdot \mathbf{\Phi }\left(\mathbf{p}\right)\right)}^{-T}$ and Q_k are significantly different⁴. Then, the non-rescaled equation with Λ = I may provide unnecessarily wide outer bounds. Numerical investigations have shown that a reasonable choice for this scaling is the block diagonal matrix

$\mathbf{\Lambda }=\mathrm{b}\mathrm{l}\mathrm{k}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\beta \mathbf{I},{\beta }^{-1}\mathbf{I}\right)(35)$

with the identity matrix $\mathbf{I}\in {\mathbb{R}}^{n_{\text{x}}\times n_{\text{x}}}$ and the square root

$\beta =\sqrt{\min \left\{{\lambda }_i\left({\mathbf{Q}}_k\right)\right\}}(36)$

of the smallest eigenvalue of Q_k. According to Rauh et al. (2021), the LMI (Eq. 33) characterizes an ellipsoid that encloses all results of the mapping (Eq. 31). To satisfy this inclusion property, the stretch parameter α_k+1 is increased up to the point, where the property of negative (semi-)definiteness holds in Eq. 33.

Step P2: Compute interval bounds for the term

${\mathbf{b}}_k=\left(\mathbf{\Phi }\left(\mathbf{p}\right)-\tilde{\mathbf{\Phi }}\right)\cdot {\mathit{\mu }}_k\in \left[{\mathbf{b}}_k\right](37)$

which accounts for a non-zero ellipsoid midpoint with x_k, $\tilde{\mathbf{\Phi }}$, and p defined according to Eqs. 27, 29, 30. Inflate the outer ellipsoid bound described by the shape matrix (Eq. 32) with

${\mathbf{Q}}_{k+1}={\left(1+{\rho }_{\text{O},k+1}\right)}^2\cdot {\check{\mathbf{Q}}}_{k+1},{\rho }_{\text{O},k+1}=\sup \left\{‖{\alpha }_{k+1}^{-1}\cdot {\mathbf{\Gamma }}_k^{-1}\cdot \left[{\mathbf{b}}_k\right]‖\right\},(38)$

where the interval-valued generalization of the Euclidean norm operator is defined in Rauh and Jaulin (2021a).

Step P3: Compute the updated ellipsoid midpoint as

${\mathit{\mu }}_{\mathbf{\Phi },k+1}=\tilde{\mathbf{\Phi }}\cdot {\mathit{\mu }}_k.(39)$

The outer ellipsoidal enclosure ${\mathcal{E}}_{\mathbf{\Phi },k+1}$ of x_Φ,k+1 at the time instant k+1 then becomes

${\mathbf{x}}_{\mathbf{\Phi },k+1}\in {\mathcal{E}}_{\mathbf{\Phi },k+1}\left({\mathit{\mu }}_{\mathbf{\Phi },k+1},{\mathbf{\Gamma }}_{\mathbf{\Phi },k+1}\right),(40)$

where

${\mathbf{\Gamma }}_{\mathbf{\Phi },k+1}={\alpha }_{k+1}\cdot \left(1+{\rho }_{\text{O},k+1}\right)\cdot \tilde{\mathbf{\Phi }}\cdot {\mathbf{\Gamma }}_k.(41)$

Step P4: Determine an ellipsoidal enclosure for the summand

${\mathbf{x}}_{\mathbf{\Psi },k+1}=\mathbf{\Psi }\left(\mathbf{p}\right)\cdot {\mathbf{u}}_k(42)$

according to

${\mathbf{x}}_{\mathbf{\Psi },k+1}\in {\mathcal{E}}_{\mathbf{\Psi },k+1}\left({\mathit{\mu }}_{\mathbf{\Psi },k+1},{\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}\right).(43)$

Details of this step are further described in Section 4.3 for the case that the system model (Eq. 22) results from the exact temporal discretization of a linear time-invariant continuous-time model with uncertain but bounded inputs according to Eq. 21.

Step P5: Compute an ellipsoidal enclosure of the Minkowski sum of the two intermediate results ${\mathcal{E}}_{\mathbf{\Phi },k+1}$ and ${\mathcal{E}}_{\mathbf{\Psi },k+1}$ according to

${\mathcal{E}}_{\mathbf{\Phi },k+1}\oplus {\mathcal{E}}_{\mathbf{\Psi },k+1}\subseteq {\mathcal{E}}_{k+1}\left({\mathit{\mu }}_{k+1},{\mathbf{\Gamma }}_{k+1}\right)(44)$

with the new midpoint

${\mathit{\mu }}_{k+1}={\mathit{\mu }}_{\mathbf{\Phi },k+1}+{\mathit{\mu }}_{\mathbf{\Psi },k+1}(45)$

and the updated shape matrix

${\mathbf{\Gamma }}_{k+1}={\mathbf{Q}}^{\frac{1}{2}},(46)$

which is given in closed-form by the nearly optimal shape matrix

${\mathbf{Q}}_{k+1}=\left(1+\frac{1}{\beta }\right)\cdot \left({\mathbf{\Gamma }}_{\mathbf{\Phi },k+1}\cdot {\mathbf{\Gamma }}_{\mathbf{\Phi },k+1}^T\right)+\left(1+\beta \right)\cdot \left({\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}\cdot {\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}^T\right)(47)$

with

$\beta =\sqrt{\frac{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left\{{\mathbf{\Gamma }}_{\mathbf{\Phi },k+1}\cdot {\mathbf{\Gamma }}_{\mathbf{\Phi },k+1}^T\right\}}{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left\{{\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}\cdot {\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}^T\right\}}}.(48)$

For a derivation of this expression, the reader is referred to Kurzhanskii and Vályi (1997), Halder (2018), and Noack et al. (2009).

4.2 Ellipsoidal Correction Step: Intersection of Ellipsoids With Different Midpoints

To perform the correction step, either for the case of intersecting bounds for the state variables which are compatible with measured data, or for intersecting different guaranteed state enclosures during the ellipsoidal enclosure approach at the same time instant (resulting from a temporal forward and backward evaluation of the system model), we employ the intersection technique for ellipsoids with different midpoints already derived in Rauh et al. (2021). This approach is an extension of the technique for computing Dikin ellipsoids which is discussed in detail in Henrion et al. (2001).

This extension consists of the following two steps:

Step C1: 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 C2: Determine the shape matrices for the outer ellipsoid bound according to the computation of Dikin ellipsoids according to Henrion et al. (2001).

Preliminary work in Rauh et al. (2021) has shown that an efficient heuristic approach for the computation of the common center point ${\tilde{\mathit{\mu }}}_{k+1}$ in Step C1 of the two ellipsoids $\mathcal{E}\left({\mathit{\mu }}_{k+1},{\mathbf{\Gamma }}_{k+1}\right)$, ${\mathbf{Q}}_{k+1}={\mathbf{\Gamma }}_{k+1}{\mathbf{\Gamma }}_{k+1}^T$ and $\mathcal{E}\left({\mathit{\mu }}_{\text{m},k+1},{\mathbf{\Gamma }}_{\text{m},k+1}\right)$, ${\mathbf{Q}}_{\text{m},k+1}={\mathbf{\Gamma }}_{\text{m},k+1}{\mathbf{\Gamma }}_{\text{m},k+1}^T$ is given by an approach motivated by the innovation step of a Kalman filter (Kalman, 1960; Stengel, 1994). Here, the ellipsoid $\mathcal{E}\left({\mathit{\mu }}_{\text{m},k+1},{\mathbf{\Gamma }}_{\text{m},k+1}\right)$ characterizes the measurement model (Eq. 23) with the output matrix C.

With the help of this information, the Kalman gain matrix

${\mathbf{L}}_{k+1}={\mathbf{Q}}_{k+1}\cdot {\mathbf{C}}^T\cdot {\left(\mathbf{C}\cdot {\mathbf{Q}}_{k+1}\cdot {\mathbf{C}}^T+{\mathbf{Q}}_{\text{m}}\right)}^{-1}(49)$

can be computed to define the updated ellipsoid midpoint

${\tilde{\mathit{\mu }}}_{k+1}={\mathit{\mu }}_{k+1}+{\mathbf{L}}_{k+1}\cdot \left({\mathit{\mu }}_{\text{m},k+1}-{\mathit{\mu }}_{k+1}\right).(50)$

Now, both ellipsoids to be intersected are enclosed during the Step C2 by new ellipsoids centered at the midpoint ${\tilde{\mathit{\mu }}}_{k+1}$. For that purpose, the scaling factors

${\zeta }_{k+1}=1+‖{\mathbf{\Delta }\mathit{\mu }}_{1,k+1}‖\text{with}{\mathbf{\Delta }\mathit{\mu }}_{1,k+1}={\mathbf{Q}}_{k+1}^{-\frac{1}{2}}\cdot \left({\tilde{\mathit{\mu }}}_{k+1}-{\mathit{\mu }}_{k+1}\right)(51)$

and

${\zeta }_{\text{m},k+1}=1+‖{\mathbf{\Delta }\mathit{\mu }}_{\text{m},k+1}‖\text{with}{\mathbf{\Delta }\mathit{\mu }}_{\text{m},k+1}={\left({\mathbf{P}}_{\text{m}}^{\prime }\right)}^{\frac{1}{2}}\cdot \left({\tilde{\mathit{\mu }}}_{k+1}-{\mathit{\mu }}_{\text{m},k+1}^{\prime }\right)(52)$

are determined which represent the maximum distances of the new midpoint computed in Eq. (50) from the original ellipsoid surfaces. Here, $‖{\mathbf{\Delta }\mathbf{\mu }}_{1,k+1}‖$ and $‖{\mathbf{\Delta }\mathbf{\mu }}_{\text{m},k+1}‖$ are the Euclidean norms of the corresponding vector-valued arguments.

The outer bound of the intersection of these two rescaled ellipsoids is given by Eq. 53 in Rauh et al. (2021) by

${\mathcal{E}}_{k+1}={\mathcal{E}}_{k+1}\left({\tilde{\mathit{\mu }}}_{k+1},{\mathbf{\Gamma }}_{k+1}^{\prime }\right)\text{with}{\mathbf{\Gamma }}_{k+1}^{\prime }={\left({\mathbf{Q}}_{k+1}^{\prime }\right)}^{\frac{1}{2}},(53)$

where the shape matrix is determined with the closed-form expression

${\mathbf{Q}}_{k+1}^{\prime }=4\cdot {\left(2\cdot \left({\left({\zeta }_{k+1}^2\cdot {\mathbf{Q}}_{k+1}\right)}^{-1}+{\zeta }_{\text{m},k+1}^{-2}\cdot {\mathbf{P}}_{\text{m}}^{\prime }\right)\right)}^{-1}(54)$

with Q_k+1 from Eq. 38. As it was shown in Rauh et al. (2021), this approach is applicable also in cases in which not all components of the state vector are measured, i.e., if the matrix ${\mathbf{P}}_{\text{m}}^{\prime }$ introduced in Eq. 24 represents a degenerate ellipsoid with bounds that may be infinitely wide in some components of the state space. Note that this procedure may be less accurate than the technique presented in Becis-Aubry (2020). However, it has the advantage of a simple closed-form representation that does not require the solution of a specific optimization task at each sampling instant at which measurements are available.

Remark 4. To avoid pessimism that may result from the ellipsoid midpoint choice in Eq. 50 when the sizes of both ellipsoids to be intersected are significantly different, the example in the following section makes use of a three-fold evaluation of the intersection step, where besides the choice in Eq. 50 also intersection results are determined which are centered either at μ_k+1 or at ${\mathit{\mu }}_{\text{m},k+1}^{\prime }$. From all three alternatives, the ellipsoid with the smallest volume is selected when continuing the evaluation of the state equations.

4.3 Application to Linear Time-Invariant Continuous-Time Models With Bounded Inputs

To make the predictor–corrector state estimation approach presented in the previous two subsections applicable to the continuous-time models introduced in Sections 2, 3, define the fundamental matrix

$\mathbf{\Phi }\left(\mathbf{p}\right)\in e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \delta },(55)$

where δ is the known discretization step size. This fundamental matrix is either evaluated in symbolic form with a subsequent replacement of all occurrences of the parameters p by their corresponding point values (or by their interval bounds $\left[\mathbf{p}\right]$ to obtain a suitable interval extension of this matrix). If symbolic formula manipulation for the computation of this matrix is inefficient in the case of large dimensions n_x, interval bounds for the matrix $\mathbf{\Phi }\left(\mathbf{p}\right)$ can be determined alternatively by solving a suitable initial value problem (IVP) in a verified way. In detail, an interval enclosures of the i-th column of $\mathbf{\Phi }\left(\mathbf{p}\right)$ can be obtained by solving the auxiliary IVP

$\dot{\mathbf{x}}\left(t\right)=\mathbf{A}\left(\mathbf{p}\right)\cdot \mathbf{x}\left(t\right)(56)$

with the initial conditions x(0) = e_i, where e_i is the i-th unit vector. Possible options for solvers applicable to this task are VNODE-LP (Nedialkov, 2011), VSPODE (Lin and Stadtherr, 2007) or CAPD (Kapela et al., 2020).

In addition, a set-valued enclosure of the matrix $\mathbf{\Psi }\left(\mathbf{p}\right)$ according to

$\mathbf{\Psi }\left(\mathbf{p}\right)\in \delta \cdot e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]}\cdot \mathbf{B}\left(\mathbf{p}\right)(57)$

is necessary to evaluate the Steps P4 and P5 of the ellipsoidal state prediction approach. As described in Step P4, cf. Eq. 42, an ellipsoidal enclosure of the expression

${\mathbf{x}}_{\mathbf{\Psi },k+1}=\mathbf{\Psi }\left(\mathbf{p}\right)\cdot {\mathbf{u}}_k(58)$

is required for that purpose, where the input signal of the dynamic system model is bounded according to ${\mathbf{u}}_k\in \left[{\mathbf{u}}_k\right]$, again without any assumptions on the temporal variation rates of this system input. To determine this ellipsoid, the following two different options are possible.

4.3.1 Approach 1: LMI-Based Enclosures

In the first approach, the bounds are computed after determining a convex polytope ${\mathbb{X}}_{\mathbf{\Psi },k+1}$ representing all possible values for x_Ψ,k+1 with the corresponding vertices ${\mathbf{x}}_{\mathbf{\Psi },k+1}^{⟨j⟩}$, j ∈ {1, … , L}. From these vertex points, an ellipsoidal enclosure

${\mathcal{E}}_{\mathbf{\Psi },k+1}\left({\mathit{\mu }}_{\mathbf{\Psi },k+1},{\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}\right)\text{with}{\mathbf{Q}}_{\mathbf{\Psi },k+1}={\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}\cdot {\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}^T(59)$

can be determined which is close to the minimum-volume Löwner-John ellipsoid can be determined by solving the LMI constrained optimization problem.⁵

$\begin{array}{l}\hfill \underset{{\mathbf{Q}}_{\mathbf{\Psi },k+1},{\mathit{\mu }}_{\mathbf{\Psi },k+1}}{\min }\left(\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left\{{\mathbf{Q}}_{\mathbf{\Psi },k+1}\right\}\right)\\ \hfill \left[\begin{array}{cc}\hfill 1\hfill & \hfill {\left({\mathbf{x}}_{\mathbf{\Psi },k+1}^{⟨j⟩}-{\mathit{\mu }}_{\mathbf{\Psi },k+1}\right)}^T\hfill \\ \hfill \left({\mathbf{x}}_{\mathbf{\Psi },k+1}^{⟨j⟩}-{\mathit{\mu }}_{\mathbf{\Psi },k+1}\right)\hfill & \hfill {\mathbf{Q}}_{\mathbf{\Psi },k+1}\hfill \end{array}\right]⪰0\text{for\hspace{0.17em}all}j\in \left\{1,\dots ,L\right\}\end{array}(60)$

${\mathbf{Q}}_{\mathbf{\Psi },k+1}⪰0.$

The LMIs above with the decision variables μ_Ψ,k+1 and Q_Ψ,k+1 are equivalent according to the Schur complement formula to the inequality constraints

${\left({\mathbf{x}}_{\mathbf{\Psi },k+1}^{⟨j⟩}-{\mathit{\mu }}_{\mathbf{\Psi },k+1}\right)}^T\cdot {\mathbf{Q}}_{\mathbf{\Psi },k+1}^{-1}\cdot \left({\mathbf{x}}_{\mathbf{\Psi },k+1}^{⟨j⟩}-{\mathit{\mu }}_{\mathbf{\Psi },k+1}\right)\le 1(61)$

that need to be satisfied for each polygon vertex j ∈ {1, … , L} to have a guaranteed outer enclosure of the polygon.

4.3.2 Approach 2: Simplified Conversion of Box-Type Enclosures Into Ellipsoids

As a computationally less demanding formulation, the shape matrix Q_Ψ,k+1 (which in the limit case may be degenerate due to zero eigenvalues leading to principal axes of vanishing length in some of the dimensions of the state space) can be approximated in terms of

${\mathbf{Q}}_{\mathbf{\Psi },k+1}=n_{k+1}^{\prime }\cdot \mathbf{T}\cdot \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{\left[\begin{array}{ccc}\hfill r_{1,k+1}^2\hfill & \hfill \dots \hfill & \hfill r_{n,k+1}^2\hfill \end{array}\right]\right\}\cdot {\mathbf{T}}^T.(62)$

Here, $n_{k+1}^{\prime }$ is the number of non-zero entries in the vector

${\mathbf{r}}_{k+1}=\left[\begin{array}{c}\hfill r_{1,k+1}\hfill \\ \hfill ⋮\hfill \\ \hfill r_{n,k+1}\hfill \end{array}\right]=\frac{1}{2}\cdot \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{m}\left\{\delta \cdot e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]-\mathrm{m}\mathrm{i}\mathrm{d}\left(\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]\right)}\cdot \mathbf{B}\left(\mathbf{p}\right)\cdot \left[{\mathbf{u}}_k\right]\right\},(63)$

where the regular preconditioning matrix

$\mathbf{T}=e^{\text{mid}\left(\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]\right)}(64)$

has been introduced to reduce the influence of overestimation due to the wrapping effect of interval analysis (Jaulin et al., 2001b).

Using the shape matrix in Eq. 62, the associated ellipsoid midpoint is given by

${\mathit{\mu }}_{\mathbf{\Psi },k+1}=\mathbf{T}\cdot \mathrm{m}\mathrm{i}\mathrm{d}\left(\delta \cdot e^{\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]-\mathrm{m}\mathrm{i}\mathrm{d}\left(\mathbf{A}\left(\left[\mathbf{p}\right]\right)\cdot \left[0;\delta \right]\right)}\cdot \mathbf{B}\left(\mathbf{p}\right)\cdot \left[{\mathbf{u}}_k\right]\right).(65)$

Remark 5. Uncertainty in the discretization step size δ can be accounted for in all equations in this subsection by treating its value as an interval parameter in analogy to the procedure described in (Rohou and Jaulin, 2021), Section 3.4.

5 Numerical Benchmark Example

5.1 Benchmark With a Point-Valued System Matrix

As a benchmark example, we consider the system model

$\left[\begin{array}{c}\hfill {\dot{x}}_1\hfill \\ \hfill {\dot{x}}_2\hfill \end{array}\right]=\left[\begin{array}{c}\hfill x_2\hfill \\ \hfill -x_1-x_2+u\hfill \end{array}\right](66)$

that has been studied originally in Rohou and Jaulin (2021) for the demonstration of the polygonal state estimation procedure. For the state Eq. 66, it is assumed that the initial conditions are unknown to the state estimator so that they have to be reconstructed by the proposed procedure during a temporal backward evaluation of the system model, starting from the first point of time at which measured data are available. For the following simulation case study, we assume further that the control input of the system is described by

$u\left(t\right)\in \cos \left(t\right)+\sin \left(\frac{t}{3}\right)+\frac{t}{10}+\left[-0.1;0.1\right],(67)$

for which an enclosing tube is determined according to Theorem 2 and Figure 2 for the polygonal estimator and by Eqs. 42, 43 for the ellipsoidal counterpart. Note that the exact temporal dependence of this system input, and therefore also its first and all higher-order time derivatives, are unknown to the state estimator. Only outer bounds for the system input, described by a tube parameterized in terms of interval slices according to Figure 2, are assumed to be available.

For the application of the state estimator, measurements are available at selected, exactly known points of time in the window $t\in \left[0;10\right]$. Those measurement instants do not necessarily coincide with the sampling interval with which the state Eq. 66 are discretized in Eq. 22. Hence, prediction steps (in temporal forward and backward direction) are required to synchronize both the prediction and measurement-based correction steps of the state estimators. As far as the measured data are concerned, interval bounds for the corresponding uncertainty are represented by

$\mathbf{x}\left(t_i\right)\in {\mathbf{y}}_{\text{m}}\left(t_i\right)+\left[\begin{array}{c}\hfill \left[-0.01;0.01\right]\hfill \\ \hfill \left[-0.01;0.01\right]\hfill \end{array}\right](68)$

with the point-valued measured data ${\mathbf{y}}_{\text{m}}\left(t_i\right)$ at points of time t_i listed in Table 1.

TABLE 1

TABLE 1. Measurement instants t_i and measured data ${\mathbf{y}}_{\text{m}}\left(t_i\right)$.

5.1.1 Short Time Intervals Between Measurements

For the scenario that the (non-equidistant) time intervals between subsequent measurements are relatively short, we consider all eight measurements given in Table 1. In this case, Figure 3 gives a comparison between the state enclosures obtained by a simple interval-based, the polygonal, and the ellipsoidal wrappers. For the ellipsoidal approach, the two different variants for the computation of ${\mathcal{E}}_{\mathbf{\Psi },k+1}\left({\mathit{\mu }}_{\mathbf{\Psi },k+1},{\mathbf{\Gamma }}_{\mathbf{\Psi },k+1}\right)$ according to Section 4.3 are compared, where the label ellipsoid (box) refers to Eqs. 62–65 and ellipsoid (LMI) to the solution of Eqs. 60, 61.

FIGURE 3

FIGURE 3. Comparison of interval-based, polygonal, and ellipsoidal wrappers for the case of eight measurement instants with different discretization step sizes δ. (A) State enclosures for δ = 0.1. (B) State enclosures for δ = 0.01. (C) Enlarged view for δ = 0.1. (D) Enlarged view for δ = 0.01. (E) Volume of the state enclosures for δ = 0.1. (F) Volume of the state enclosures for δ = 0.01.

The comparison of the enclosures of all feasible states on the equidistant time grid k⋅δ as well as the comparison of the resulting volumes of the state enclosures illustrates that all results are close to each other in this scenario due to the short temporal distances between the individual measurements. As illustrated in Rohou and Jaulin (2021), the setting δ = 0.01 represents a solution that produces almost the optimal bounds for the domains of reachable states, while δ = 0.1 can be used as a practical compromise between the tightness of the state enclosures and the computational effort.

Correlations between the state variables x₁ and x₂ are equally well represented by the polygonal and by the ellipsoidal approach, see the enlarged view in Figures 3E,F. The corresponding enlarged domain is indicated by blue boxes in the first two graphs of Figure 3. It should be pointed out that axis-aligned interval boxes are not capable of representing such dependencies. As motivated earlier in this paper, this phenomenon coincides with the wrapping effect of interval analysis when considering larger time spans for the evaluation of the state equations without a measurement-based reduction of the size of the interval boxes. From a practical point of view, this may lead to excessively conservative interval enclosures. This effect is investigated further in the following subsection.

Remark 6. The larger volumes of the ellipsoidal state enclosures in comparison with the polygonal ones (cf. Figures 3E,F) are mainly caused by the fact that the measured state information according to Eq. 68 are originally given in the form of interval boxes which need to be enclosed by ellipsoids that — in the considered two-dimensional scenario — have a volume that is larger by a factor of $\frac{\pi }{2}\approx 1.57$.

Remark 7. The polygonal state estimator is implemented in such a way that the product $\mathbf{\Psi }\left(\mathbf{p}\right)\cdot {\mathbf{u}}_k$ is enclosed by an interval box in an intermediate evaluation step. Therefore, it is possible that the LMI-based ellipsoidal enclosure approach leads to tighter bounds of this term if strong correlations exist between the individual vector components, for example, if the true set representation of this product resembles a kind of non-axis-aligned parallelepiped.

5.1.2 Long Time Intervals Between Measurements

To create a simulation scenario with longer time intervals between the measurement points, the simulation of the previous subsection is modified in such a way that only the time instants $t_i\in \{\frac{2}{3},9.0\}$ of Table 1 are taken into consideration.

In Figure 4, all solution approaches are depicted for the discretization step size δ = 0.1, while Figure 5 shows the results for δ = 0.01. As in the previous subsection, the polygonal and ellipsoidal solutions are close to each other, where the simplified box-type representation of the input term produces only small overestimation in comparison with the optimized LMI solution. Both polygonal and ellipsoidal enclosures are successfully applicable to capture correlations between the state variables x₁ and x₂, where the larger volumes of the ellipsoidal approach are again explained by the more conservative representation of the interval boxes describing the measured data.

FIGURE 4

FIGURE 4. Comparison of interval-based, polygonal, and ellipsoidal wrappers for the case of two measurement instants with the discretization step size δ = 0.1. (A) Comparison of the state enclosures with an interval-based solution. (B) Polygonal and ellipsoidal state enclosures. (C) State enclosures (enlarged view). (D) Volume of the state enclosures.

FIGURE 5

FIGURE 5. Comparison of interval-based, polygonal, and ellipsoidal wrappers for the case of two measurement instants with the discretization step size δ = 0.01. (A) Comparison of the state enclosures with an interval-based solution. (B) Polygonal and ellipsoidal state enclosures. (C) State enclosures (enlarged view). (D) Volume of the state enclosures.

However, due to the long time spans between measurements, the naive interval wrapper is much more pessimistic. Even though the box widths reduce slightly for smaller values of δ, the shape of the actual state trajectory cannot be traced by the interval wrapper which produces boxes that are overlapping over the whole time interval $t\in \left[0;10\right]$. Note that the visualization of the computed volumes also confirms the exponential growth of the state enclosures in pure interval-based simulations (Figures 4D, 5D), that are well known in the literature, cf. Lohner (1987) and Lohner (2001).

5.2 Benchmark With an Uncertain System Matrix

To illustrate that the ellipsoidal state estimation approach remains applicable in settings with time-invariant interval parameters in the state equations without any algorithmic changes, consider the modified system model

$\left[\begin{array}{c}\hfill {\dot{x}}_1\hfill \\ \hfill {\dot{x}}_2\hfill \end{array}\right]=\left[\begin{array}{c}\hfill x_2\hfill \\ \hfill -\left(1+p_1\right)\cdot x_1-\left(1+p_1\right)\cdot x_2+u\hfill \end{array}\right],(69)$

where $p_1\in \left[-0.1;0.1\right]$ and $p_2\in \left[-0.1;0.1\right]$ are two independent time-invariant uncertain parameters.

For the case of the discretization step size δ = 0.1, Figure 6 gives a comparison of the ellipsoidal state enclosures in the cases of the uncertain parameters and for p₁ = p₂ = 0. Both for short and long time spans between the measurements (chosen as in the previous subsections), the ellipsoidal bounds for the uncertain state equations remain much tighter than the box-type solution enclosures described above (even though they were determined for a model with exactly known parameters). This confirms the statement of very limited overestimation introduced by the ellipsoidal enclosure technique. Moreover, it needs to be pointed out that also the computing times remain practically identical as before. The only step that is now slightly more costly is the computation of the matrix exponentials involved in Eq. 11. However, these matrices can be pre-computed so that they do not impact the potential real-time capability of the state estimation algorithm.

FIGURE 6

FIGURE 6. Comparison of the ellipsoidal wrappers for the case of two and eight measurement instants with the discretization step size δ = 0.1, without and with interval uncertainty in the parameters p₁ and p₂: dark gray – ellipsoidal bounds for the system model (69) with interval parameters; light gray – ellipsoidal bounds for the system model (66) with point-valued parameters. (A) Scenario with two measurement instants. (B) Scenario with eight measurement instants.

Remark 8. In this simulation case study, the simplified box-type representation of the input term produces practically the same bounds as the optimized LMI solution. Therefore, only the latter is depicted in Figure 6.

6 Conclusion and Outlook on Future Work

In this paper, a novel ellipsoidal state estimation procedure has been presented which can be employed as a computationally efficient alternative to an optimal polygonal approach published in previous work. Pessimism in the ellipsoidal approach mostly results from the representation of box-type uncertainty by a common outer ellipsoidal hull. Our contribution is especially relevant for lots of applications that require computational efficiency. For instance, in embedded systems — such as autonomous robots — performing online state estimation is an omnipresent task. For them, it is important to efficiently compute guaranteed outer enclosures of the sets of reachable states even if they are not represented by minimum volume bounds. This “trade-off” is actually crucial and the test cases provided in this paper reveal that our approach is relevant, providing near-optimal enclosures with a fast algorithm.

In future work, this pessimism can be reduced by treating the interval-valued measurement uncertainty individually for each scalar output in terms of degenerate ellipsoids as already demonstrated in Rauh et al. (2021). However, in this case, it will be necessary to derive systematic criteria that allow the user to select the best uncertainty representation (wrt. computational effort and tightness of the enclosures) prior to the simulation so that tedious trial-and-error studies can be avoided which aim at figuring out whether the intersection with a joint ellipsoidal uncertainty representation or multiple intersections with degenerate ellipsoids lead to tighter bounds.

Moreover, combinations of the proposed state estimation procedure with online adaptations of control strategies will be a subject of future work. These include trajectory re-planning procedures for autonomous robots if collisions with obstacles cannot be ruled out with certainty or — more generally — the combination of decentralized adaptive control procedures for distributed systems with set-valued state and disturbance estimation procedures.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding authors.

Author Contributions

The algorithm was designed and implemented by AR and SR. The paper was jointly written by all authors. All authors have read and agreed to the published version of the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹Interval observers provide a kind of framer for the actual system behavior by computing lower and upper bounding trajectories for all possible states in a decoupled manner. Due to the fact that they mimic the structure of a Luenberger observer for both lower and upper bounding systems, they should not be confused with predictor–corrector state estimation schemes in which a prediction based on the open-loop dynamics is performed between the discrete time instants at which measured data are available. Typically, the prediction in predictor–corrector state estimation schemes is carried out without any assumption regarding potential cooperativity of the state equations.

²The dot notation (⋅) is used in this paper for representing the independent variable.

³For ease of reading, the notation u(⋅) is used instead of ${\mathbf{u}}_{\left[t_1;t_2\right]}(\cdot )$, as the values of u(⋅) are independent of those outside the interval $\left[t_1;t_2\right]$.

⁴Due to the fact that the matrix $\mathbf{\Phi }\left(\mathbf{p}\right)$ is defined as the exponential of a continuous-time model’s system matrix according to Eq. 21, the eigenvalues of $\tilde{\mathbf{\Phi }}=\mathbf{\Phi }\left(\mathrm{m}\mathrm{i}\mathrm{d}\left(\left[\mathbf{p}\right]\right)\right)$ are always non-zero and, thus, $\tilde{\mathbf{\Phi }}$ is invertible. This is a direct consequence of the series expansion definition of the matrix exponential and can be shown alternatively by Jacobi’s formula to prove the identity $\det \left(e^{\mathbf{M}}\right)=e^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left\{\mathbf{M}\right\}}$ for arbitrary square complex matrices M, where the matrix determinant equals the product of its eigenvalues (Hall, 2015).

⁵The exact volume minimization task would require the solution of an optimization task in which the minimization of a matrix trace is replaced by the more complex determinant minimization task described in Appendix C of Tarbouriech et al. (2011).

References

Amann, N., Bocker, J., and Prenner, F. (2004). Active Damping of Drive Train Oscillations for an Electrically Driven Vehicle. IEEE/ASME Trans. Mechatron. 9, 697–700. doi:10.1109/tmech.2004.839036

CrossRef Full Text | Google Scholar

Anderson, B. D. O., and Moore, J. B. (2005). Optimal Filtering. Mineola, USA: Dover Publications, Inc.

Google Scholar

Aubin, J., and Frankowska, H. (1990). Set-Valued Analysis. Boston: Birkhäuser.

Google Scholar

Becis-Aubry, Y. (2020). Ellipsoidal Constrained State Estimation in Presence of Bounded Disturbances. Preprint Available at: arxiv.org/pdf/2012.03267v1.pdf accessed: May 25, 2021).

Google Scholar

Bertsekas, D., and Rhodes, I. (1971). Recursive State Estimation for a Set-Membership Description of Uncertainty. IEEE Trans. Automat. Contr. 16, 117–128. doi:10.1109/tac.1971.1099674

CrossRef Full Text | Google Scholar

Berz, M., and Makino, K. (1998). Verified Integration of ODEs and Flows Using Differential Algebraic Methods on High-Order Taylor Models. Reliable Comput. 4, 361–369. doi:10.1023/a:1024467732637

CrossRef Full Text | Google Scholar

Brogliato, B., and Tanwani, A. (2020). Dynamical Systems Coupled with Monotone Set-Valued Operators: Formalisms, Applications, Well-Posedness, and Stability. SIAM Rev. 62, 3–129. doi:10.1137/18m1234795

CrossRef Full Text | Google Scholar

Davis, J. (2002). “Luenberger Observers,” in Foundations of Deterministic and Stochastic Control. Systems & Control: Foundations & Applications (Boston, MA, USA: Birkhäuser).

CrossRef Full Text | Google Scholar

Deville, Y., Janssen, M., and van Hentenryck, P. (2002). Consistency Techniques for Ordinary Differential Equations. Constraint 7, 289–315. doi:10.1023/a:1020573518783

CrossRef Full Text | Google Scholar

Di Bernardo, M., di Gaeta, A., Hoyos Velasco, C., Montanaro, U., and Santini, S. (2012). “Model-Based Soft Landing Control of an Electromechanical Engine Valve Actuator,” in Volume 2: Legged Locomotion; Mechatronic Systems; Mechatronics; Mechatronics for Aquatic Environments; MEMS Control; Model Predictive Control; Modeling and Model-Based Control of Advanced IC Engines, 87–94. doi:10.1115/dscc2012-movic2012-8526

CrossRef Full Text | Google Scholar

Eble, I. (2007). Über Taylor-Modelle. Ph.D. thesis. German: Universität Karlsruhe (TH), Fakultät für Mathematik. In German.

Google Scholar

Efimov, D., Raïssi, T., Chebotarev, S., and Zolghadri, A. (2013). Interval State Observer for Nonlinear Time Varying Systems. Automatica 49, 200–205. doi:10.1016/j.automatica.2012.07.004

CrossRef Full Text | Google Scholar

Filippov, A. F. (1988). Differential Equations with Discontinuous Righthand Sides. Dordrecht, Netherlands: Springer-Verlag.

Google Scholar

Goldsztejn, A., and Neumaier, A. (2014). On the Exponentiation of Interval Matrices. Reliable Comput. 20, 52–72.

Google Scholar

Hairer, E., Nørsett, S., and Wanner, G. (2000). Solving Ordinary Differential Equations I. 2nd edn. Berlin Heidelberg, Germany: Springer-Verlag.

Google Scholar

Halder, A. (2018). “On the Parameterized Computation of Minimum Volume Outer Ellipsoid of Minkowski Sum of Ellipsoids,” in Proc. of IEEE Conference on Decision and Control (CDC), Miami Beach, USA, Dec. 2018 (IEEE), 4040–4045. doi:10.1109/CDC.2018.8619508

CrossRef Full Text | Google Scholar

Hall, B. (2015). “Lie Groups, Lie Algebras, and Representations: An Elementary Introduction,” in Graduate Texts in Mathematics (Cham, Switzerland: Springer International Publishing).

Google Scholar

Haus, B., Mercorelli, P., and Werner, N. (2014). “A Piezo Servo Hydraulic Actuator for Use in Camless Combustion Engines and its Control with MPC,” in Proceeding of the 2014 International Conference on Control, Decision and Information Technologies (CoDIT), Metz, France, Nov. 2014 (IEEE), 471–476.

CrossRef Full Text | Google Scholar

Henrion, D., Tarbouriech, S., and Arzelier, D. (2001). LMI Approximations for the Radius of the Intersection of Ellipsoids: Survey. J. Optimization Theor. Appl. 108, 1–28. doi:10.1023/a:1026454804250

CrossRef Full Text | Google Scholar

Hoefkens, J. (2001). Rigorous Numerical Analysis with High-Order Taylor Models. Ph.D. thesis. Michigan State University. Available from: http://www.bt.pa.msu.edu/cgi-bin/display.pl?name=hoefkensphd (Accessed May 25, 2021).

Google Scholar

Jaulin, L., Kieffer, M., Braems, I., and Walter, E. (2001a). Guaranteed Non-linear Estimation Using Constraint Propagation on Sets. Int. J. Control. 74, 1772–1782. doi:10.1080/00207170110090642

CrossRef Full Text | Google Scholar

Jaulin, L., Kieffer, M., Didrit, O., and Walter, É. (2001b). Applied Interval Analysis. London: Springer-Verlag.

Google Scholar

John, F. (1948). “Extremum Problems with Inequalities as Subsidiary Conditions,” in Studies and Essays Presented to R. Courant on His 60th Birthday (New York, USA: Interscience Publishers, Inc.), 187–204.

Google Scholar

Kalman, R. E. (1960). A New Approach to Linear Filtering and Prediction Problems. Transaction ASME – J. Basic Eng. 82, 35–45. doi:10.1115/1.3662552

CrossRef Full Text | Google Scholar

Kapela, T., Mrozek, M., Wilczak, D., and Zgliczynski, P. (2020). CAPD:DynSys: A Flexible C++ Toolbox for Rigorous Numerical Analysis of Dynamical Systems. Commun. Nonlinear Sci. Numer. Simulation 101, 105578. doi:10.1016/j.cnsns.2020.105578

CrossRef Full Text | Google Scholar

Krasnochtanova, I., Rauh, A., Kletting, M., Aschemann, H., Hofer, E. P., and Schoop, K.-M. (2010). Interval Methods as a Simulation Tool for the Dynamics of Biological Wastewater Treatment Processes with Parameter Uncertainties. Appl. Math. Model. 34, 744–762. doi:10.1016/j.apm.2009.06.019

CrossRef Full Text | Google Scholar

Kurzhanski, A. B., and Filippova, T. F. (1993). “On the Theory of Trajectory Tubes - A Mathematical Formalism for Uncertain Dynamics, Viability and Control,” in Advances in Nonlinear Dynamics and Control: A Report from Russia. Editor A. B. Kurzhanski (Boston, MA, USA: Birkhäuser), 122–188. doi:10.1007/978-1-4612-0349-0_4

CrossRef Full Text | Google Scholar

Kurzhanskii, A. B., and Vályi, I. (1997). “Ellipsoidal Calculus for Estimation and Control,” in Systems & Control: Foundations & Applications (Boston, MA, USA: Birkhäuser). doi:10.1007/978-1-4612-0277-6

CrossRef Full Text | Google Scholar

Lin, Y., and Stadtherr, M. A. (2007). Validated Solutions of Initial Value Problems for Parametric ODEs. Appl. Numer. Math. 57, 1145–1162. doi:10.1016/j.apnum.2006.10.006

CrossRef Full Text | Google Scholar

Lohner, R. (1987). “Enclosing the Solutions of Ordinary Initial and Boundary Value Problems,” in Computer Arithmetic: Scientific Computation and Programming Languages. Editors E. W. Kaucher, U. W. Kulisch, and C. Ullrich (Stuttgart: Wiley-Teubner Series in Computer Science), 255–286.

Google Scholar

Lohner, R. (1988). Einschließung der Lösung gewöhnlicher Anfangs- und Randwertaufgaben und Anwendungen. Ph.D. thesis. German: Universität Karlsruhe (TH), Fakultät für Mathematik. In German.

Google Scholar

Lohner, R. J. (2001). “On the Ubiquity of the Wrapping Effect in the Computation of Error Bounds,” in Perspectives on Enclosure Methods. Editors U. Kulisch, R. Lohner, and A. Facius (Wien, New York: Springer-Verlag), 201–216. doi:10.1007/978-3-7091-6282-8_12

CrossRef Full Text | Google Scholar

Luenberger, D. G. (1964). Observing the State of a Linear System. IEEE Trans. Mil. Electron. 8, 74–80. doi:10.1109/tme.1964.4323124

CrossRef Full Text | Google Scholar

Marouani, G., Dinh, T. N., Raïssi, T., Wang, X., and Messaoud, H. (2021). Unknown Input Interval Observers for Discrete-Time Linear Switched Systems. Eur. J. Control. 59, 165–174. doi:10.1016/j.ejcon.2020.09.004

CrossRef Full Text | Google Scholar

Mayer, G. (2017). “Interval Analysis and Automatic Result Verification,” in De Gruyter Studies in Mathematics (Berlin/Boston: De Gruyter), 516.

Google Scholar

Moore, R. E., Kearfott, R. B., and Cloud, M. J. (2009). Introduction to Interval Analysis. Philadelphia: SIAM.

Google Scholar

Mousavi, M. H., Karami, M. E., Ahmadi, M., Sharafi, P., and Veysi, F. (2020). Robust Speed Controller Design for Permanent Magnet Synchronous Motor Based on Gain-Scheduled Control Method via LMI Approach. SN Appl. Sci. 2, 1699. doi:10.1007/s42452-020-03453-z

CrossRef Full Text | Google Scholar

Nedialkov, N. (1999). Computing Rigorous Bounds on the Solution of an Initial Value Problem for an Ordinary Differential Equation. Ph.D. thesis. Graduate Department of Computer Science, University of Toronto.

Google Scholar

Nedialkov, N. (2007). “Interval Tools for ODEs and DAEs,” in CD-Proceedings of the 12th GAMM-IMACS International Symposium on Scientific Computing, Computer Arithmetic, and Validated Numerics SCAN 2006, Duisburg, Germany, Sept. 2006 (IEEE Computer Society).

Google Scholar

Noack, B., Klumpp, V., and Hanebeck, U. (2009). “State Estimation with Sets of Densities Considering Stochastic and Systematic Errors,” in Proceedings of the 2009 12th International Conference on Information Fusion, FUSION 2009, Seattle, WA, USA, July 2009 (IEEE), 1751–1758.

Google Scholar

Nedialkov, N. (2011). “Implementing a Rigorous ODE Solver through Literate Programming,” in Modeling, Design, and Simulation of Systems with Uncertainties. Editors A. Rauh, and E. Auer (Berlin, Heidelberg/Germany: Springer-Verlag). Mathematical Engineering. doi:10.1007/978-3-642-15956-5_1

CrossRef Full Text | Google Scholar

Raïssi, T., Efimov, D., and Zolghadri, A. (2012). Interval State Estimation for a Class of Nonlinear Systems. IEEE Trans. Automat. Contr. 57, 260–265. doi:10.1109/tac.2011.2164820

CrossRef Full Text | Google Scholar

Rauh, A., Bourgois, A., and Jaulin, L. (2021). Union and Intersection Operators for Thick Ellipsoid State Enclosures: Application to Bounded-Error Discrete-Time State Observer Design. Algorithms 14, 88. doi:10.3390/a14030088

CrossRef Full Text | Google Scholar

Rauh, A., and Jaulin, L. (2021a). A Computationally Inexpensive Algorithm for Determining Outer and Inner Enclosures of Nonlinear Mappings of Ellipsoidal Domains. Appl. Math. Comp. Sci. AMCS 31, 399–415. doi:10.34768/amcs-2021-0027

CrossRef Full Text | Google Scholar

Rauh, A., and Jaulin, L. (2021b). “A Novel Thick Ellipsoid Approach for Verified Outer and Inner State Enclosures of Discrete-Time Dynamic Systems,” in Proc. of 19th IFAC Symposium System Identification: Learning Models for Decision and Control, Padova, Italy. (online).

CrossRef Full Text | Google Scholar

Rauh, A., and Kersten, J. (2021). Transformation of Uncertain Linear Systems with Real Eigenvalues into Cooperative Form: The Case of Constant and Time-Varying Bounded Parameters. Algorithms 14, 85. doi:10.3390/a14030085

CrossRef Full Text | Google Scholar

Rauh, A., Kersten, J., and Aschemann, H. (2015). “Robust Control for a Spatially Three-Dimensional Heat Transfer Process,” in Proceedings of 8th IFAC Symposium on Robust Control Design ROCOND’15, Bratislava, Slovakia, July 2015, 8–11.

CrossRef Full Text | Google Scholar

Rauh, A., Westphal, R., and Aschemann, H. (2013). “Verified Simulation of Control Systems with Interval Parameters Using an Exponential State Enclosure Technique,” in CD-Proc. of IEEE Intl. Conference on Methods and Models in Automation and Robotics MMAR, Miedzyzdroje, Poland, Aug. 2013 (IEEE).

CrossRef Full Text | Google Scholar

Rauh, A., Westphal, R., Aschemann, H., and Auer, E. (2016). “Exponential Enclosure Techniques for Initial Value Problems with Multiple Conjugate Complex Eigenvalues,” in Proceedings of 16th GAMM-IMACS International Symposium on Scientific Computing, Computer Arithmetic, and Validated Numerics SCAN2014. vol. 9553 of Lecture Notes in Computer Science, April 2016, 87–122.

CrossRef Full Text | Google Scholar

Rohou, S., and Jaulin, L. (2021). Exact Bounded-Error Continuous-Time Linear State Estimator. Syst. Control. Lett. 153, 104951. doi:10.1016/j.sysconle.2021.104951

CrossRef Full Text | Google Scholar

Rohou, S., Jaulin, L., Mihaylova, L., Le Bars, F., and Veres, S. M. (2017). Guaranteed Computation of Robot Trajectories. Robotics Autonomous Syst. 93, 76–84. doi:10.1016/j.robot.2017.03.020

CrossRef Full Text | Google Scholar

Smith, H. (1995). Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, Vol 41. Renton, Washington DC, USA: Providence Mathematical Surveys and Monographs, American Mathematical Soc.

Google Scholar

Stengel, R. (1994). Optimal Control and Estimation. New York, USA: Dover Publications, Inc.

Google Scholar

Stewart, D. E. (2011). Dynamics with Inequalities. Philadelphia, USA: Society for Industrial and Applied Mathematics.

Google Scholar

Tarbouriech, S., Garcia, G., Gomes da Silva, J., and Queinnec, I. (2011). Stability and Stabilization of Linear Systems with Saturating Actuators. London, UK: Springer-Verlag. doi:10.1007/978-0-85729-941-3

CrossRef Full Text | Google Scholar

Zhang, X., Ma, J., Huang, S., Cheng, Z., and Lee, T. H. (2019). “Integrated Planning and Control for Collision-free Trajectory Generation in 3D Environment with Obstacles,” in Proc. of the 19th Intl. Conference on Control, Automation and Systems (ICCAS), Jeju, Korea (South), Oct. 2019 (IEEE), 974–979.

CrossRef Full Text | Google Scholar