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

- Lab-STICC, ROBEX, ENSTA Bretagne, Brest, France

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)^{1}). 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

where *i* ∈ {1, … , *n*_{p}}. In addition to the dynamic system model above, we assume that direct measurements of the state vector *t*_{k}, where

### 2.1 State Tube

A tube ^{2} (Kurzhanski and Filippova (1993)) is used to continuously enclose the feasible states of a dynamic system. For that purpose,

Assume that for **x**(⋅) is known to be inside a *prior tube* **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 **x**(⋅) that is consistent with both the prior tube

For a given initial vector **x**_{1} defined at *t*_{1}, the state at time *t*_{2} is expressed by

where **u**(⋅). The flow can be extended to sets (Aubin and Frankowska (1990)) as follows^{3}:

Furthermore, when the input **u**(⋅) is uncertain but known to be inside a tube

We now define the *posterior state tube* as the tightest tube **x**(⋅) consistent with the prior tube **u**(⋅) in

as illustrated in Figure 1.

**FIGURE 1**. Illustration of one slice of a posterior tube *t*_{2} (red hatched part). This slice is defined as the intersection set *t* to *t*_{2}, according to the feasible inputs *t*_{1} and from *t*_{3}, is depicted in blue at the time instant *t*_{2}.

The state estimation consists in the approximation of

**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

### 2.2 Recursive Procedure

We now define the sets **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

**Theorem 1. **Given the sampling times of measurements

a prior tube **x**(⋅), and a tube

with

### 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

##### 2.3.1 Principle

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

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.,

Moreover, given an interval matrix [**A**] and a set of vectors **A**⋅**x** assuming that **A** ∈ [**A**], i.e.,

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

a prior tube **x**(⋅), and a piecewise constant tube

containing **u**(⋅), where

with *t* is not consistent with the tube discretization, i.e., for

**FIGURE 2**. Illustration of a tube *u*]_{k=2} slice containing all feasible values for *u*(*t*),

### 2.4 Wrappers for Enclosing ${\mathbb{X}}_{t}$

Now that we have a reliable definition for the enclosure **u**](⋅) and guaranteed to enclose **u**(⋅), it remains to reliably implement the sets **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 *polygonal sequence* as a specific formulation of the general approach described in Eq. 15:

• First, forward in time propagation, for

• then, backward in time propagation, for

• finally, polygonal enclosure between different sampling points,

From Theorem 2, we know that for all

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

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

These reasons motivate the study of ellipsoids as an alternative wrapper for

## 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

is given, where

In addition to the dynamic system model above, we assume that measurements **y**_{m,k+1} with an ellipsoidal uncertainty model

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

In this inequality (Eq. 24), the 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 *n*_{x}-dimensional state space.

Using this notation, entries in the vector *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

with the positive definite shape matrix

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

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

For finding outer ellipsoidal enclosures of the subexpression

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:

Step P1: Apply

to the ellipsoid

where *α*_{k+1} ≥ 0 is the smallest value for which the linear matrix inequality (LMI)

is satisfied for all

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 **Q**_{k} are significantly different^{4}. 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

with the identity matrix

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

which accounts for a non-zero ellipsoid midpoint with **x**_{k}, **p** defined according to Eqs. 27, 29, 30. Inflate the outer ellipsoid bound described by the shape matrix (Eq. 32) with

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

The outer ellipsoidal enclosure **x**_{Φ,k+1} at the time instant *k*+1 then becomes

where

Step P4: Determine an ellipsoidal enclosure for the summand

according to

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

with the new midpoint

and the updated shape matrix

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

with

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 **C**.

With the help of this information, the Kalman gain matrix

can be computed to define the updated ellipsoid midpoint

Now, both ellipsoids to be intersected are enclosed during the Step C2 by new ellipsoids centered at the midpoint

and

are determined which represent the maximum distances of the new midpoint computed in Eq. (50) from the original ellipsoid surfaces. Here,

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

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

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

**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

### 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

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 *n*_{x}, interval bounds for the matrix *i*-th column of

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

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

is required for that purpose, where the input signal of the dynamic system model is bounded according to

##### 4.3.1 Approach 1: LMI-Based Enclosures

In the first approach, the bounds are computed after determining a convex polytope **x**_{Ψ,k+1} with the corresponding vertices *j* ∈ {1, … , *L*}. From these vertex points, an ellipsoidal enclosure

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

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

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

Here,

where the regular preconditioning matrix

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

**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

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

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

with the point-valued measured data *t*_{i} listed in Table 1.

**TABLE 1**. Measurement instants *t*_{i} and measured data

##### 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 *ellipsoid (box)* refers to Eqs. 62–65 and *ellipsoid (LMI)* to the solution of Eqs. 60, 61.

**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*_{1} and *x*_{2} 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

**Remark 7. **The polygonal state estimator is implemented in such a way that the product

##### 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

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*_{1} and *x*_{2}, 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**. 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**. 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

### 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

where

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*_{1} = *p*_{2} = 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**. 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*_{1} and *p*_{2}: 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

^{1}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.

^{2}The *dot* notation (⋅) is used in this paper for representing the independent variable.

^{3}For ease of reading, the notation **u**(⋅) is used instead of **u**(⋅) are independent of those outside the interval

^{4}Due to the fact that the matrix **M**, where the matrix determinant equals the product of its eigenvalues (Hall, 2015).

^{5}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

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

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).

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

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

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).

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

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

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.

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

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

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

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

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

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

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.

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.

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

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

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

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

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

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

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.

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).

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.

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

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

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

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

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).

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

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.

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).

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.

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

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

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.

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

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

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.

Keywords: linear time-invariant systems, bounded uncertainty, state estimation, ellipsoidal enclosures, differential inclusions

Citation: Rauh A, Rohou S and Jaulin L (2022) An Ellipsoidal Predictor–Corrector State Estimation Scheme for Linear Continuous-Time Systems With Bounded Parameters and Bounded Measurement Errors. *Front. Control. Eng.* 3:785795. doi: 10.3389/fcteg.2022.785795

Received: 29 September 2021; Accepted: 24 January 2022;

Published: 25 March 2022.

Edited by:

Prashant Mhaskar, McMaster University, CanadaReviewed by:

Soorathep Kheawhom, Chulalongkorn University, ThailandYuanlong Li, Shanghai Jiao Tong University, China

Copyright © 2022 Rauh, Rohou and Jaulin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Andreas Rauh, Andreas.Rauh@interval-methods.de; Simon Rohou, simon.rohou@ensta-bretagne.fr

^{†}**Present address:** Andreas Rauh, Group: Distributed Control in Interconnected Systems, Department of Computing Science, Carl von Ossietzky Universität Oldenburg, Oldenburg, Germany