A tube X ( ⋅ ) ² (Kurzhanski and Filippova (1993)) is used to continuously enclose the feasible states of a dynamic system. For that purpose, X ( ⋅ ) is defined as a set of trajectories X ( ⋅ ) : R ↦ P ( R n x ) which satisfies x ⋅ ∈ X ⋅ ⇔ ∀ t , x t ∈ X t . (2)

Assume that for t ∈ 0 ; t ̄ , the state trajectory x(⋅) is known to be inside a prior tube X ̌ ( ⋅ ) . 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 X ( ⋅ ) for x(⋅) that is consistent with both the prior tube X ̌ ( ⋅ ) (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 x 2 = Φ t 1 , t 2 u ⋅ x 1 , (3) where Φ 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 ³ : Φ t 1 , t 2 u ⋅ X 1 = x 2 | ∃ x 1 ∈ X 1 , x 2 = Φ t 1 , t 2 u ⋅ x 1 . (4)

Furthermore, when the input u(⋅) is uncertain but known to be inside a tube U ( ⋅ ) , the set flow given by Eq. 4 becomes Φ t 1 , t 2 U ⋅ X 1 = ⋃ u ⋅ ∈ U ⋅ Φ t 1 , t 2 u ⋅ X 1 . (5)

We now define the posterior state tube as the tightest tube X ̂ ( ⋅ ) for x(⋅) consistent with the prior tube X ̌ ( ⋅ ) , the feasible inputs u(⋅) in U ( ⋅ ) , and the dynamic system model. We obtain X ̂ t = ⋃ u ⋅ ∈ U ⋅ ⋂ τ ∈ 0 ; t ̄ Φ τ , t u ⋅ X ̌ τ (6) as illustrated in Figure 1.

The state estimation consists in the approximation of X ̂ ( ⋅ ) , ∀ t ∈ 0 ; t ̄ . 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 B p ⋅ u by means of a suitable reformulation.

We now define the sets X → X → t and 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 X ̂ ( ⋅ ) . 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 T = 0 , δ , 2 δ , … , k ̄ δ = t 0 , t 1 , t 2 , … , t k ̄ , (7) a prior tube X ̌ ( ⋅ ) containing the state trajectory x(⋅), and a tube U ( ⋅ ) for the inputs. The posterior tube X ̂ ( ⋅ ) can be defined recursively by X t k = Φ t k − 1 , t k U ⋅ X t k − 1 ∩ ⋃ u ⋅ ∈ U ⋅ ⋂ τ ∈ t k − 1 ; t k Φ τ , t k u ⋅ X ̌ τ X ⃖ t k = Φ t k + 1 , t k U ⋅ X ⃖ t k + 1 ∩ ⋃ u ⋅ ∈ U ⋅ ⋂ τ ∈ t k ; t k + 1 Φ τ , t k u ⋅ X ̌ τ X ̂ t k = X t k ∩ X ⃖ t k (8) with X t 0 = X ̌ ( t 0 ) and X ⃖ t k ̄ = X ̌ ( t k ̄ ) . Moreover, for t ∈ t k ; t k + 1 \ T , we have X ̂ t = Φ t k , t U ⋅ X ̂ t k . (9)

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 x ̇ = f ( x , 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.

Now that we have a reliable definition for the enclosure X ( ⋅ ) from a bounded input U ( ⋅ ) (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 X ( ⋅ ) . Several wrappers can be considered, such as intervals, polygons, or ellipsoids. The use of intervals amounts to representing 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.

Consider a polygon P k which contains 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 ∈ { 1 , … , k ̄ } :

From Theorem 2, we know that for all t k ∈ T , X ̂ t k ⊂ P k , ∀ t ∈ t k ; t k + 1 , X ̂ t ⊂ P k : k + 1 . (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 Φ p = e A p ⋅ δ and Ψ p = δ ⋅ e A p ⋅ δ ⋅ B p (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.

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.

These reasons motivate the study of ellipsoids as an alternative wrapper for X ( ⋅ ) sets in the following section.

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 Φ p ⋅ 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 p 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 Φ p ⋅ x k in Eq. 22, we assume that this part of the system model is reformulated in the form x Φ , k + 1 = Φ p ⋅ x k = Φ p ⋅ x k + Φ ̃ ⋅ μ k + Φ p − Φ ̃ ⋅ μ 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: E k = E k μ k , Γ k denotes the ellipsoidal uncertainty on the non-origin centered states x k , (27) E k = E k 0 , Γ k the uncertainty of x ̌ k after shifting the ellipsoid to the origin, and (28) Φ ̃ = Φ m i d p the midpoint approximation of the quasi - linear system matrix with (29) p ∈ p = p ̲ ; p ̄ , and the interval midpoint m i d p = 1 2 ⋅ p ̲ + p ̄ . (30)

Step P1: Apply

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 Φ ̃ − 1 ⋅ Φ p − 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 Λ = b l k d i a g β I , β − 1 I (35) with the identity matrix I ∈ R n x × n x and the square root β = min λ i Q k (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

The outer ellipsoidal enclosure E Φ , k + 1 of x _Φ,k+1 at the time instant k+1 then becomes x Φ , k + 1 ∈ E Φ , k + 1 μ Φ , k + 1 , Γ Φ , k + 1 , (40) where Γ Φ , k + 1 = α k + 1 ⋅ 1 + ρ O , k + 1 ⋅ Φ ̃ ⋅ Γ k . (41)

Step P4: Determine an ellipsoidal enclosure for the summand

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 E Φ , k + 1 and E Ψ , k + 1 according to E Φ , k + 1 ⊕ E Ψ , k + 1 ⊆ E k + 1 μ k + 1 , Γ k + 1 (44) with the new midpoint μ k + 1 = μ Φ , k + 1 + μ Ψ , k + 1 (45) and the updated shape matrix Γ k + 1 = Q 1 2 , (46) which is given in closed-form by the nearly optimal shape matrix Q k + 1 = 1 + 1 β ⋅ Γ Φ , k + 1 ⋅ Γ Φ , k + 1 T + 1 + β ⋅ Γ Ψ , k + 1 ⋅ Γ Ψ , k + 1 T (47) with β = t r a c e Γ Φ , k + 1 ⋅ Γ Φ , k + 1 T t r a c e Γ Ψ , k + 1 ⋅ Γ Ψ , k + 1 T . (48)

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

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;

Preliminary work in Rauh et al. (2021) has shown that an efficient heuristic approach for the computation of the common center point μ ̃ k + 1 in Step C1 of the two ellipsoids E μ k + 1 , Γ k + 1 , Q k + 1 = Γ k + 1 Γ k + 1 T and E μ m , k + 1 , Γ m , k + 1 , Q m , k + 1 = Γ m , k + 1 Γ 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 E μ m , k + 1 , Γ m , k + 1 characterizes the measurement model (Eq. 23) with the output matrix C.

With the help of this information, the Kalman gain matrix L k + 1 = Q k + 1 ⋅ C T ⋅ C ⋅ Q k + 1 ⋅ C T + Q m − 1 (49) can be computed to define the updated ellipsoid midpoint μ ̃ k + 1 = μ k + 1 + L k + 1 ⋅ μ m , k + 1 − μ k + 1 . (50)

Now, both ellipsoids to be intersected are enclosed during the Step C2 by new ellipsoids centered at the midpoint μ ̃ k + 1 . For that purpose, the scaling factors ζ k + 1 = 1 + Δ μ 1 , k + 1 with Δ μ 1 , k + 1 = Q k + 1 − 1 2 ⋅ μ ̃ k + 1 − μ k + 1 (51) and ζ m , k + 1 = 1 + Δ μ m , k + 1 with Δ μ m , k + 1 = P m ′ 1 2 ⋅ μ ̃ k + 1 − μ m , k + 1 ′ (52) are determined which represent the maximum distances of the new midpoint computed in Eq. (50) from the original ellipsoid surfaces. Here, Δ μ 1 , k + 1 and Δ μ 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 E k + 1 = E k + 1 μ ̃ k + 1 , Γ k + 1 ′ with Γ k + 1 ′ = Q k + 1 ′ 1 2 , (53) where the shape matrix is determined with the closed-form expression Q k + 1 ′ = 4 ⋅ 2 ⋅ ζ k + 1 2 ⋅ Q k + 1 − 1 + ζ m , k + 1 − 2 ⋅ P m ′ − 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 P m ′ 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 μ m , k + 1 ′ . From all three alternatives, the ellipsoid with the smallest volume is selected when continuing the evaluation of the state equations.

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 Φ p ∈ e A p ⋅ δ , (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 p 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 Φ p 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 Φ p can be obtained by solving the auxiliary IVP x ̇ t = A p ⋅ x t (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 Ψ p according to Ψ p ∈ δ ⋅ e A p ⋅ 0 ; δ ⋅ B p (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 x Ψ , k + 1 = Ψ p ⋅ u k (58) is required for that purpose, where the input signal of the dynamic system model is bounded according to u k ∈ u k , again without any assumptions on the temporal variation rates of this system input. To determine this ellipsoid, the following two different options are possible.

As a benchmark example, we consider the system model x ̇ 1 x ̇ 2 = x 2 − x 1 − x 2 + u (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 t ∈ cos t + sin t 3 + t 10 + − 0.1 ; 0.1 , (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 ∈ 0 ; 10 . 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 x t i ∈ y m t i + − 0.01 ; 0.01 − 0.01 ; 0.01 (68) with the point-valued measured data y m t i at points of time t _i listed in Table 1.

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 x ̇ 1 x ̇ 2 = x 2 − 1 + p 1 ⋅ x 1 − 1 + p 1 ⋅ x 2 + u , (69) where p 1 ∈ − 0.1 ; 0.1 and p 2 ∈ − 0.1 ; 0.1 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.

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.