# Offline and Online Use of Interval and Set-Based Approaches for Control and State Estimation: A Selection of Methodological Approaches and Their Application

Andreas Rauh, Marit Lahme, Simon Rohou, Luc Jaulin, Thach Ngoc Dinh, Tarek Raissi, Mohamed Fnadi

► To cite this version:

Andreas Rauh, Marit Lahme, Simon Rohou, Luc Jaulin, Thach Ngoc Dinh, et al. Offline and Online Use of Interval and Set-Based Approaches for Control and State Estimation: A Selection of Methodological Approaches and Their Application. Logical Methods in Computer Science, 2025, 21 (4), 10.46298/lmcs-21(4:11)2025. [https://arxiv.org/abs/2004.04117](https://arxiv.org/abs/2004.04117)

HAL Id: hal-05329422

https://cnam.hal.science/hal-05329422v1

Submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# OFFLINE AND ONLINE USE OF INTERVAL AND SET-BASED APPROACHES FOR CONTROL AND STATE ESTIMATION: A SELECTION OF METHODOLOGICAL APPROACHES AND THEIR APPLICATION

ANDREAS RAUH ⚫<sup>a</sup>, MARIT LAHME ⚫<sup>a</sup>, SIMON ROHOU ⚫<sup>b</sup>, LUC JAULIN ⚫<sup>b</sup>, THACH NGOC DINH ⚫<sup>c</sup>, TAREK RAÏSSI ⚫<sup>c</sup>, AND MOHAMED FNADI ⚫<sup>d</sup>

**ABSTRACT.** Control and state estimation procedures need to be robust against imprecisely known parameters, uncertainty in initial conditions, and external disturbances. Interval methods and other set-based techniques form the basis for the implementation of powerful approaches that can be used to identify parameters of dynamic system models in the presence of the aforementioned types of uncertainty. Moreover, they are applicable to a verified feasibility and stability analysis of controllers and state estimators. In addition to these approaches which are typically used offline for analysis of system models designed with classical floating point procedures, interval and set-based methods have also been developed in recent years, which allow to directly solve the associated design tasks and to implement reliable techniques that are applicable online, i.e., during system operation. The latter approaches include set-based model predictive control, online parameter adaptation techniques for nonlinear variable-structure and backstepping controllers, interval observers, and fault diagnosis techniques. This paper provides an overview of the methodological background and reviews numerous practical applications for which interval and other set-valued approaches have been employed successfully.

*Key words and phrases:* Interval analysis, uncertain systems, parameter identification, state estimation, guaranteed stabilizing control, model predictive control, interval observers.

<figure><img src="image_2.png" /></figure>

<footnote><sup>a</sup> Distributed Control in Interconnected Systems, Department of Computing Science, Carl von Ossietzky Universität Oldenburg, D-26111 Oldenburg, Germany
*e-mail address:* andreas.rauh@uni-oldenburg.de, marit.lahme@uni-oldenburg.de</footnote>

<footnote><sup>b</sup> ENSTA, Lab-STICC (Robex), F-29806 Brest, France
*e-mail address:* simon.rohou@ensta.fr, luc.jaulin@ensta.fr</footnote>

<footnote><sup>c</sup> Cedric-Laetitia, Conservatoire National des Arts et Métiers, F-75141 Paris, France
*e-mail address:* ngoc-thach.dinh@lecnam.net, tarek.raissi@cnam.fr</footnote>

<footnote><sup>d</sup> Laboratoire d'Informatique, Signal et Image de la Côte d'Opale, Université du Littoral Côte d'Opale, F-62228, France
*e-mail address:* mohamed.fnadi@univ-littoral.fr</footnote>
# 1. INTRODUCTION

Interval analysis has become an active field of research over the past decades. Especially, the foundations for its modern use can be traced back to the works of R.E. Moore and his famous book published in 1966 [Moo66]. However, even earlier evidence from the 19th century can be found if not only works aiming at computer implementations are accounted for but if also foundations of truncated, finite-precision arithmetic are considered [Kea23, Sch79].

## 1.1. Use of Interval Methods to Deal with Bounded Uncertainty.
In its widest sense, interval analysis provides an approach to perform computations on finitely large sets in such a way that the determined solutions contain all the solutions of the mathematically stated problem with certainty [May17, JKDW01]. In such a way, it is not only possible to account for numerical errors due to the limited precision of the numbers representable in floating point arithmetic on a CPU, but also for discretization and method errors. The latter arise inevitably when solving, for example, ordinary differential equations with the help of temporal series expansions truncated after a finite order [Ned07, Ned11, KMWZ20, Loh87, Bü20]. In addition to these effects, which are small (at least) for a single evaluation step, but can cause severe damage if neglected, the same techniques also allow for computing with wider sets that result from imperfectly known system parameters and initial conditions in (technical) applications. This lack of knowledge is commonly referred to as epistemic uncertainty [KD09]. As such, interval analysis provides a mathematically sound way to handle epistemic uncertainty that can furthermore also result from other kinds of a lack of knowledge of the actual system behavior, for example, friction phenomena in dynamics of rigid bodies or fluidic applications.

In addition to the application domains mentioned above, further high-level functionalities have been developed on the basis of interval analysis. These contain techniques for solving sets of nonlinear algebraic equations with the help of Newton-type methods [Kra69, Neu90] or the identification of feasible solution sets that satisfy (nonlinear) constraints in a guaranteed manner [ATDC14]. In the latter class of problems, the SIVIA algorithm [JKDW01, WJ94] (set inversion via interval analysis) as well as contractor approaches [ATDC14] serve as the fundamental solution methodologies, commonly being combined with Newton-type methods. Moreover, extensions exist that allow for handling outliers in a suitable manner by relaxed set intersections [Jau09].

## 1.2. Notations.
Throughout this article, real-valued scalar interval variables are denoted in the form

$$x \in [x] := [\underline{x} ; \bar{x}] \qquad (1.1)$$

for which $\underline{x} \le \bar{x}$ holds. The interval bounds $\underline{x}, \bar{x} \in \mathbb{R}$ are commonly also denoted as the infimum $\inf([x])$ and supremum $\sup([x])$. To distinguish vectors and matrices from scalar quantities, boldface lower- and upper-case letters are used, respectively. Interval vectors and matrices, denoted, for example, by $[\mathbf{a}]$ and $[\mathbf{A}]$, are defined in an elementwise manner as the Cartesian product of scalar intervals.

For a set-based generalization of the evaluation of a (multi-variate, vector-valued) function $\mathbf{y} = \mathbf{f}(\mathbf{x}) : \mathbb{R}^n \to \mathbb{R}^m$, the so-called interval extension $[\mathbf{y}] = [\mathbf{f}]([\mathbf{x}])$ is defined [JKDW01, Moo66, May17]. It satisfies the property

$$[\mathbf{f}]([\mathbf{x}]) \supset \{\mathbf{f}(\mathbf{x}) | \mathbf{x} \in ([\mathbf{x}])\} . \qquad (1.2)$$
### 1.3. Challenges of the Use of Set-Based Methods.
Besides the advantageous feature that interval techniques always provide guaranteed enclosures of a mathematically formulated problem in the sense (1.2), they also come with the disadvantage of a certain degree of pessimism. This pessimism is caused by the following two fundamental issues [JKDW01, May17, CL84]:

* The dependency effect: The difference of two intervals with numerically identical bounds generally does not simplify to the value zero as shown in

$$[0; 1] - [0; 1] = [-1; 1]. \qquad (1.3)$$

This is due to the fact that the difference in the equation above does not carry any information on whether both arguments are identical quantities or not. Therefore, variables are generally treated as independent in any arithmetic operation according to

$$\{x \circ y | x \in [x] := [\underline{x}; \bar{x}], y \in [y] := [\underline{y}; \bar{y}]\}, \quad \text{where } \circ \in \{+, -, \cdot, /\}, \qquad (1.4)$$

unless simplifications can be carried out on the basis of symbolic formula manipulation, which allows to exploit knowledge that the two arguments $x$ and $y$ represent exactly the same (physical) quantity. Therefore, in classical interval analysis, there does not exist any inversion for the operations of set-based addition and multiplication.

* The wrapping effect: Transformations of sets, for example, rotations and scalings even in linear matrix-vector products as in the recursive evaluation of discrete-time dynamic systems

$$[\mathbf{x}](t_{k+1}) = \mathbf{A} \cdot [\mathbf{x}](t_k) \qquad (1.5)$$

with the exemplarily chosen initial conditions and dynamics matrix

$$\mathbf{x}(t_0) \in [\mathbf{x}](t_0) = \begin{bmatrix} [-1; 1] \\ [-1; 1] \end{bmatrix}, \quad \mathbf{A} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ -1 & 1 \end{bmatrix} \qquad (1.6)$$

may lead to an exponential increase in the volume of the solution enclosures

$$\mathbf{x}(t_k) \in 2^{\frac{k}{2}} \cdot \begin{bmatrix} [-1; 1] \\ [-1; 1] \end{bmatrix} \qquad (1.7)$$

if axis-aligned interval box enclosures were determined in each step (1.5). This growth appears despite the fact that the true volume would stay constant in the example above because the matrix $\mathbf{A}$ in (1.6) represents a volume-preserving rotation of the initial box by 45° in each step. This preservation of volume is especially visible in even sampling steps, where for example the relations

$$\mathbf{A}^2 = \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \quad \text{and} \quad \mathbf{A}^4 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \qquad (1.8)$$

hold.

### 1.4. Selected Countermeasures Against Overestimation.
In the example (1.5)-(1.7), changing the evaluation order according to

$$[\mathbf{x}](t_{k+1}) = \mathbf{A}^{k+1} \cdot [\mathbf{x}](t_0), \qquad (1.9)$$

where the integer-order powers of the point-valued matrix $\mathbf{A}$ are firstly determined, allows to solve this problem [KRK$^{+}$10]. However, there is no general solution available, especially in cases where the transformations are nonlinear or the matrices $\mathbf{A}$ may themselves include intervals in each element.
Then, a large variety of other set-representations may turn out to be more efficient, for example, the use of ellipsoids [KV97, BA21] (being mapped exactly to an ellipsoid in point-valued linear transformations), (constrained) zonotopes [IPI+22, RLRR22], zonotope bundles [AK11], or polytopes [RRJ22, RJ21b]. In nonlinear settings, the computationally inexpensive thick ellipsoid approach recently published by two of the authors provides the possibility to detect pessimism in the guaranteed outer enclosure by performing a comparison with an inner ellipsoidal enclosure belonging to the solution set with certainty [RBJ21, RJ21a]. Moreover, computationally more demanding techniques such as Taylor model arithmetic [Bü20, Bü18, Bü21, Hoe01] provide further means to reduce the wrapping effect that is omnipresent in any dynamic system simulation or during finding solutions to algebraic systems of equations. Taylor model arithmetic has also been employed in the frame of guaranteed state estimation, especially in predictor-corrector approaches, as well as in the (offline) verification of the robustness of control strategies for nonlinear uncertain systems [KRHA06, AK12, KAH07]. Although Taylor models allow for effectively reducing the wrapping effect by representing nonlinear dependencies between the individual state variables, they have the following drawbacks in an online use: First, intersections between predicted domains and information obtained by sensors is computationally expensive during state estimation. Second, strategies for limiting the Taylor model order are inevitable to limit and keep the computational complexity constant during online use. For that reasons, zonotopes as well as ellipsoids have shown to represent a good compromise with respect to a slight decrease in the state enclosure quality but a notable reduction in computational effort.

As an additional option, coordinate transformations, as derived in [KRA18, KRA19, MB10] for systems with either real or complex eigenvalues after their linearization, as well as the exploitation of specific monotonicity properties, allow for reducing both the pessimism discussed above and the computational effort required for its reduction. An example for the exploitation of monotonicity properties is the design of interval observers [REZ12, ERCZ13] that allows for computing elementwise lower and upper bounds for the state trajectories according to the following scheme [Smi95]. Here and in the remainder of the article, the word observer denotes an estimation scheme that allows for reconstructing information about not directly measurable internal state variables by means of a model of the considered dynamic system as well as information about its control signals and measured outputs.

An autonomous set of ordinary differential equations

$$ \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t)), \mathbf{x} \in \mathbb{R}^n \qquad (1.10) $$

is cooperative according to [AS03, Smi95, Hir19] if the property

$$ x_i^{(1)}(t) \geq x_i^{(2)}(t) \qquad (1.11) $$

holds at each $t \geq t_0$ for all elements $i \in \{1, \dots, n\}$ of the two vectors with $\mathbf{x}^{(1)}(t_0)$ and $\mathbf{x}^{(2)}(t_0)$ being the initial conditions which satisfy again the inequalities

$$ x_i^{(1)}(t_0) \geq x_i^{(2)}(t_0) \quad \text{for all } i \in \{1, \dots, n\}. \qquad (1.12) $$

For continuous-time systems, cooperativity can be checked by the sufficient sign conditions

$$ J_{i,j}(\mathbf{x}) \geq 0, \quad i, j \in \{1, \dots, n\}, \quad i \neq j \qquad (1.13) $$

for all off-diagonal elements of the Jacobian

$$ \mathbf{J}(\mathbf{x}) = \frac{\partial \mathbf{f}}{\partial \mathbf{x}}(\mathbf{x}) \qquad (1.14) $$
of the right-hand side of the state equations evaluated for all reachable states $\mathbf{x} = \mathbf{x}(t)$. Matrices satisfying this non-negativity property for the off-diagonal elements are also denoted as *Metzler* matrices in the literature [Smi95, RKA20a, REZ12, RE18].

In the case of linear uncertain systems, this cooperativity property simplifies to an elementwise sign condition for the off-diagonal elements of the lower interval bounds of the uncertain dynamics matrix $\mathbf{A} \in [\underline{\mathbf{A}}; \overline{\mathbf{A}}]$ in the dynamic system model

$$ \dot{\mathbf{x}}(t) = \mathbf{A} \cdot \mathbf{x}(t) . \qquad (1.15) $$

Guaranteed bounds for the state trajectories $\mathbf{x}(t)$ can be obtained in the form $[\mathbf{x}] (t) = [\mathbf{v}(t); \mathbf{w}(t)]$ so that all reachable states are included by an evaluation of the following coupled set of state equations (which are a direct consequence of [Mül27])

$$ \dot{v}_i(t) = \inf ([A_{i,i}] \cdot v_i(t)) + \inf \left( \sum_{\substack{j=1 \\ j \neq i}}^{n} [A_{i,j}] \cdot [x_j](t) \right) \qquad (1.16) $$

$$ \dot{w}_i(t) = \sup ([A_{i,i}] \cdot w_i(t)) + \sup \left( \sum_{\substack{j=1 \\ j \neq i}}^{n} [A_{i,j}] \cdot [x_j](t) \right) $$

with the resulting state bounds

$$ [x_i](t) = [v_i(t); w_i(t)] . \qquad (1.17) $$

Note that couplings between the vectors $\mathbf{v}(t)$ and $\mathbf{w}(t)$ disappear in both summations for $j \neq i$ included in (1.16) if the system is positive, i.e., $v_i(t) \ge 0$ holds for all $t \ge t_0$, see for example [REZ12, ERCZ13]. There, this property is exploited for the design of a Luenberger-type interval observer, which allows for guaranteed state estimation in the case of bounded measurement uncertainty.

Moreover, recent developments of the so-called mixed monotonicity approach and remainder form decompositions [KY23, KSY22] can be exploited as a relaxation of the aforementioned evaluation technique based on cooperativity. In this case, a cooperative subsystem is extracted from the dynamic system model, while the non-cooperative part is enclosed with the help of an additive bounded error term. In such a way, the possibility for extracting bounding systems is preserved and the need for computing with general sets is still avoided.

For scenarios, in which overestimation cannot be avoided fully, the authors of [RAA10] suggest the use of systematic criteria to quantify the worst-case level of arising pessimism in the computed state enclosures. Such criteria can be embedded into the development of automatic recommender systems that support users to select the most appropriate set-based solution procedure from a variety of different alternatives [AR12].

### 1.5. Software Libraries for Interval and Set-Based Calculus.
For the practical implementation of set-valued calculus, numerous software libraries exist. The following short list is given without any claim for completeness: INTLAB [Rum99] for MATLAB and OCTAVE, the OCTAVE interval package [Hei18], VERSOFT [Roh19], JULIAINTERVALS [SB14], C-XSC [Krä12], DYNIBEX [AdSC16], ARIADNE [CBGV12], the interval arithmetic package in MATHEMATICA [Kei93], the CORA library [Alt15], and CODAC [RD24]. Many of these libraries provide a large number of basic operations, higher level functionalities such as solving
linear or nonlinear algebraic sets of equations as well as (global) optimization problems,
solving initial value problems for ordinary differential equations, or procedures for the
verification of hybrid dynamic models arising in the field of cyber-physical systems. Note,
specific approaches that are equally applicable for offline and online use for control design
and implementation according to the following subsection are mostly left to the interested
users who have to implement these procedures by themselves.

**1.6. Distinction Between the Offline and Online Use of Set-Based Methods.** As far as control engineering applications are concerned, set-valued calculus is generally applicable during two different phases. On the one hand, an *offline* use is possible. In this frame, *offline* means the use during the design phase to verify certain system properties that should be respected by, for example, a control law or state estimation procedure to be developed. On the other hand, the *online* use aims at exploiting set-valued calculus at run-time of a dynamic system, typically coming along with constraints for real-time capability. Then, possible tasks are real-time capable gain scheduling, state estimation, or fault detection.

**1.7. Outline of this Article.** Sec. 2 gives an overview of the fundamentals of identifying parameters for dynamic systems with the help of interval and other set-based approaches. These identification results are the basis for the variable-structure control approaches reviewed in Sec. 3 as well as for the guaranteed model predictive control approaches in Sec. 4. Thereafter, we present approaches for set-based state estimation in Sec. 5, applicable to the identification of the open-circuit voltage characteristic of Lithium-ion batteries, and an ellipsoidal technique in Sec. 6 that allows for treating set-valued and stochastic uncertainty and disturbances in a unified manner when solving state estimation tasks. Finally, conclusions and a brief outlook on future work are given in Sec. 7.

## 2. SYSTEM MODELING AND VERIFIED PARAMETER IDENTIFICATION

The estimation of the parameters of both finite-dimensional sets of discrete-time difference equations as well as continuous-time ordinary differential equations is crucial for the design and implementation of all model-based control and state estimation procedures. Identification experiments are typically performed (for asymptotically stable systems) with predefined, time-dependent actuation signals so that measurements of selected outputs can be gathered at certain discrete points of time $t_k, k \in \mathbb{N}_0$. Here, $t_0 = 0$ is assumed without loss of generality. Due to the collection of measurements at discrete points of time, it is necessary that parameter identification routines allow for a propagation of state information between two subsequent measurement instants $t_{k-1}$ and $t_k, k \ge 1$. The following description of the set-based parameter identification approach is based on the more detailed discussions that can be found in [RKA20b].
## 2.1. Predictor-Corrector Framework for the Verified Identification of Parameters of Dynamic Systems.
As summarized in Fig. 1, a first class of verified parameter identification routines is characterized by an observer-based structure. This structure is closely linked to the evaluation steps that are performed in classical state estimators such as Luenberger-type observers or (Extended) Kalman Filters [Kal60, JUDW00]. Hence, prediction and correction steps are separated from each other and evaluated in an alternating manner.

In this framework, the detection of feasible parameter domains is performed by means of the following three-stage procedure [RDA11, DARA13]:

(1) A verified evaluation of the set of state equations is performed between two subsequent measurement points. In the case of discrete-time state equations, the corresponding expressions are evaluated with the help of interval arithmetic (or one of the alternative set-valued approaches mentioned in the introduction of this article), while verified solvers for initial value problems for ordinary differential equations [AC16, LS07, KMWZ20, Ned11, ARHL08, Ned99, AdSC16, Ned07, Küh99] are required if continuous-time processes are taken into consideration. In any case, problem-specific approaches for the reduction of overestimation [RKA11] — such as the exploitation of cooperativity, interval subdivision, preconditioning of state equations, symbolic simplifications, or enclosures allowing to trace the shape of the reachable state domains in a tight manner — shall be employed to minimize the pessimism in the computed state and output boundaries that is caused by both the dependency and wrapping effects [JKDW01, Loh87].

(2) Information resulting from the prediction of state intervals up to the point at which new measurements are available are combined with information obtained from sensors. The applicable techniques are:

* An intersection of the enclosures of predicted state variables with interval bounds for directly measured state intervals (in terms of an additive superposition of a point-valued measurement with the respective tolerance interval);

* The use of contractors and set inversion approaches to determine an inverse mapping from the measured outputs to the internal system states and parameters [IPI$^{+22}$. Examples are forward-backward contractors [JKDW01] initialized with the predicted state intervals that are intersected with the associated measurement outputs or interval Newton techniques in the case of nonlinear output equations or for output equations involving more than one state variable. Additionally, also observer-based formulations for solving the same task were recently published in [IPI$^{+22}$], where zonotopes are employed as the underlying uncertainty representation.

(3) Parameter subdomains are eliminated either by using contractors or after a subdivision of parameter intervals with associated multiple evaluations of the state equations. In the latter case, a parameter subinterval is guaranteed to be inconsistent with the information provided by the specified system model and the available measurements if the result of the state prediction substituted into the measurement model and the interval-valued measured data do not overlap. For a visualization of this exclusion text, see Fig. 2, where a branch-and-bound procedure on the basis of multiple simulations of the state equations over the complete horizon of gathered measured data is depicted.

## 2.2. Simulation-Based Parameter Identification Procedures.
Simulation-based parameter identification procedures, making use of multiple evaluations of the set of state
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Predictor-corrector scheme for combined state estimation and parameter identification [RKA20b].</figcaption></figure>

equations for subintervals of the complete initial domain of possible system parameters, are one of the available options for parameter identification.

<figure><img src="image_4.png" /><figcaption>(A) Measurement information of the system output $y_m(t)$, sampled at the equidistant instants $t_0, t_1, t_2, \dots$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(B) Output trajectory including interval uncertainty for a possibly feasible (undecided) system parameterization (Case 1).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(C) Output trajectory including interval uncertainty for a guaranteed feasible system parameterization (Case 2).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(D) Output trajectory including interval uncertainty for a guaranteed infeasible system parameterization (Case 3).</figcaption></figure>

<caption>FIGURE 2. Simulation-based distinction between possibly feasible, guaranteed feasible, and infeasible system parameterizations [DARA13, RKA20b].</caption>

The underlying branch-and-bound procedure<fnref n="1" /> relies on the following assumptions:

* measured data are available at discrete points of time,

* worst-case bounds for measurement tolerances are known a-priori (modeling of the sensor uncertainty),

* the model structure to be parameterized is assumed to be structurally correct,

* outer enclosures for the domains of possibly uncertain initial states are given, and

<footnote n="1">¹Branching of the parameter domain by an interval subdivision approach with a subsequent simulation-based bounding of the output trajectories.</footnote>
* interval bounds on the uncertain parameters are known in terms of conservative overapproximations.

To exclude parts of the domains of uncertain initial states and parts of the a-priori given parameter intervals, a subdivision procedure of the respective domains is executed [CR97]. To detect intervals for the parameterization of the system model which are guaranteed to be inconsistent, directly measured and simulated state intervals are intersected at identical points of time. If state variables are not directly measured, the contractor procedures mentioned in the previous subsection are equally applicable here. Parameter interval subdivisions are performed according to the following distinction of cases, where the directions in which the parameter intervals are subdivided are determined with the help of a sensitivity analysis [DARA13]. For accelerations by both additional (physics-inspired) contractors and GPU-based parallelizations, the readers are referred to [RKA20b, ARK20]:

**Case 1:** Parameter subintervals are yet undecided if their corresponding predicted output intervals overlap at least partially for all available sensors at each point of time at which measurements are available and are never outside the range of measured data.

**Case 2:** A parameter subinterval is guaranteed to be consistent if the corresponding output intervals are true subintervals of the interval-valued measured data for each of the available sensors at each point of time. This subinterval is no longer investigated in the exclusion strategy but stored in a list of guaranteed admissible boxes.

**Case 3:** A parameter subinterval is guaranteed to be inconsistent if the corresponding output interval lies outside the range of the interval-valued measured data. This has to hold true for at least one of the available sensors and for at least one point in time. The resulting parameter subinterval is excluded from further evaluations and the underlying simulation can be aborted at that point in time as which the inconsistency is firstly detected.

### 3. INTERVAL-BASED VARIABLE-STRUCTURE CONTROL

As mentioned in the previous section, verified parameter identification is the prerequisite for the implementation of guaranteed stabilizing control procedures. These include linear matrix inequality methods for systems with polytopic and norm-bounded parameter uncertainty [BEFB94, RKA20a], interval-based gain scheduling techniques [KRA20], and variable structure as well as backstepping-type controllers [SRA16, RS16, RSKA16, RSA15] without and with hard state constraints.

In this section, we review a robust control procedure that extends the sliding mode design methodology by means of the online use of interval analysis to prevent the violation of state constraints with certainty [SRA16, RS16, RSA16]. For a similar approach, applicable to the backstepping design methodology, the reader is referred to [RKA17].

#### 3.1. Goal of the Control Synthesis.
Variable-structure control procedures aim at solving two-fold control tasks. First, they ensure that a system state is guided towards a time-dependent desired trajectory within finite time. Second, they guarantee that this desired trajectory is tracked with high accuracy as soon as it has been reached. A possible application scenario is heating up a solid oxide fuel cell stack from ambient temperature to its usual operating point with temperatures above 800 K along a pre-defined trajectory and keeping the temperature in close vicinity of that operating point despite variations of the electric load
taken from the fuel cell. Accurate temperature control, including a guaranteed prevention of a rise of the temperature over a maximum admissible value is crucial to simultaneously ensure a high conversion efficiency of the chemical energy of hydrogen used as the fuel gas and simultaneously avoiding accelerated system degradation and hazardous states caused by overtemperatures [RS16, RSA15].

**3.2. Fundamental Sliding Mode Control Laws.** Both the treatment of strict inequality constraints representing the limitations of the admissible operating conditions such as the temperature limit above and bounded interval uncertainty of internal system parameters can be combined with classical first- and second-order sliding mode techniques [Eke10, FL02]. These classical approaches have the drawback that they suppose perfectly known parameters and can therefore typically not ensure the guaranteed satisfaction of state constraints. Moreover, they may be inefficient in the use of the available ranges of actuator signals [RA12].

**3.2.1. First-Order Sliding Mode Control.** For a summary of the fundamental sliding mode control design of single-input single-output systems, consider the n-th order set of ordinary differential equations to be given in the nonlinear controller canonical form

$$ 
\begin{aligned}
\dot{\mathbf{x}}(t) &= [\dot{x}_1(t) \quad \dots \quad \dot{x}_{n-1}(t) \quad \dot{x}_n(t)]^T \\
&= [x_2(t) \quad \dots \quad x_n(t) \quad a(\mathbf{x}(t), \mathbf{p}) + b(\mathbf{x}(t), \mathbf{p}) \cdot v(t)]^T
\end{aligned}
\qquad (3.1) $$

with the state vector $\mathbf{x}(t) \in \mathbb{R}^n$ and the control input $v(t) \in \mathbb{R}$. To ensure controllability, we assume $b(\mathbf{x}(t), \mathbf{p}) \neq 0$. Under this assumption, a feedback linearizing control strategy

$$ v(t) = \frac{-a(\mathbf{x}(t), \mathbf{p}) + u(t)}{b(\mathbf{x}(t), \mathbf{p})} \qquad (3.2) $$

can be designed. It aims at compensating the nonlinear expressions $a(\mathbf{x}(t), \mathbf{p})$ and $b(\mathbf{x}(t), \mathbf{p})$ which both depend on the state $\mathbf{x}(t)$ and internal parameters $\mathbf{p}$. However, in reality, these functions are imperfectly known so that the auxiliary input $u(t)$ needs to ensure both robustness and stability despite uncertainty.

For this system model, the output variable to be controlled is given by

$$ y(t) = x_1(t) . \qquad (3.3) $$

Due to the assumed structure of the system (3.1), the output variable (3.3) has the relative degree $n$ [Mar03], corresponding to the property that the n-th time derivative $x_1^{(n)}(t)$ of the system output is the lowest-order derivative explicitly depending on either of the control inputs $v(t)$ and $u(t)$. Therefore, the output $y(t)$ corresponds to a (trivial) flat system<fnref n="2" /> output [FLMR95]. This property of differential flatness ensures that an exact tracking of a sufficiently smooth, i.e., n times continuously differentiable, desired trajectory $x_{1,d}(t)$ can be achieved in the case of exactly known parameters.

On the basis of this desired trajectory, the corresponding tracking error

$$ \xi_1(t) = x_1(t) - x_{1,d}(t) \qquad (3.4) $$

<footnote n="2">²A system is differentially flat if both the control signal and the state trajectories can be expressed purely in terms of a flat output (not necessarily coinciding with the output to be controlled) and a finite number of time derivatives thereof.</footnote>
and its r-th time derivative<fnref n="3" /> are introduced. To represent both in a compact form, we exploit the notation

$$ \xi_1^{(r)}(t) = x_1^{(r)}(t) - x_{1,d}^{(r)}(t), \quad r \in \{0, 1, \dots, n\}, \qquad (3.5) $$

to define the sliding surface<fnref n="4" />

$$ s(t) = \sum_{r=0}^{n-1} \alpha_r \xi_1^{(r)}(t) \qquad (3.6) $$

in which the highest-order coefficient is normalized according to $\alpha_{n-1} = 1$. Note that in the case of perfect trajectory tracking, the equality $s(t) = 0$ holds.

To guarantee asymptotic stability for states exactly located on this sliding surface, the parameters $\alpha_r$ have to fulfill the necessary and sufficient stability conditions for a Hurwitz polynomial of the order $n-1$.

A classical first-order sliding mode control can be derived with the help of the quadratic radially unbounded candidate for a positive definite Lyapunov function<fnref n="5" />

$$ V^{(1)}(t) = \frac{1}{2} s^2(t), \qquad (3.7) $$

where $s(t)$ depends on the tracking error and its derivatives and therefore is a function of the system state, see also the definition (3.1).

For $s \neq 0$, $V^{(1)} > 0$ holds. State trajectories leading to the desired system behavior $V = 0$ can be obtained if (global) asymptotic stability of the dynamic system can be proven. This stability requirement corresponds to the (global) negative definiteness of the corresponding time derivative

$$ \dot{V}^{(1)}(t) = s(t) \cdot \dot{s}(t) = \left( \sum_{r=0}^{n-1} \alpha_r \xi_1^{(r)}(t) \right) \cdot \left( \sum_{r=0}^{n-1} \alpha_r \xi_1^{(r+1)}(t) \right) < 0, \qquad (3.8) $$

where $\dot{s}$ explicitly depends on the control input.

The stability requirement above ensures that the system state converges towards the desired reference if $s \neq 0$, either due to initial conditions not perfectly coinciding with the reference or due to the occurrence of disturbances.

Following the derivation of the variable-structure control law according to [RSKA16, RSA15], the right-hand side of the inequality (3.8) is replaced by the more conservative formulation

$$ \left( \sum_{r=0}^{n-1} \alpha_r \xi_1^{(r)} \right) \cdot \left( \sum_{r=0}^{n-1} \alpha_r \xi_1^{(r+1)} \right) < -\eta \cdot |s| = -\eta \cdot \left( \sum_{r=0}^{n-1} \alpha_r \xi_1^{(r)} \right) \cdot \operatorname{sign}(s) \qquad (3.9) $$

which guarantees global asymptotic stability for $\eta > 0$. Here, the actual choice of the constant $\eta$ significantly influences the dynamics and the maximum absolute values of the control signal in the reaching phase characterized by $s \neq 0$. As soon as $s = 0$ has been reached after a finite time duration, the control amplitudes depend (in the unperturbed case) on the actual choice of the reference trajectory $x_{1,d}$.

<footnote n="3">³Time arguments are omitted subsequently, whenever the meaning is non-ambiguous.</footnote>

<footnote n="4">⁴In this context, sliding means that the system dynamics obey a differential equation of order smaller than $n$ as soon as the state is bound to the surface $s = 0$, cf. [Utk92].</footnote>

<footnote n="5">⁵The superscript (I) for the Lyapunov function and the associated control law serve as a distinction from the second-order procedure discussed in the following subsection.</footnote>
The time-independent coefficients $\alpha_r$ have a major influence on the control amplitudes if non-modeled errors and disturbances act on the system dynamics and if the error signals $\xi_1^{(r)}$ are corrupted by non-negligible measurement noise or state reconstruction errors, preventing the perfect tracking of the desired trajectory $x_{1,d}(t)$ and making it impossible that $s \equiv 0$ holds perfectly after the end of the reaching phase.

The derivation of the control law is completed by enforcing that the second factor in (3.9) becomes proportional to the sign of the actual value of $s$ according to

$$ \sum_{r=0}^{n-2} \alpha_r \xi_1^{(r+1)} + u - x_{1,d}^{(n)} + \eta \cdot \mathrm{sign}(s) = -\beta \cdot \mathrm{sign}(s) \quad \text{with an arbitrary } \beta > 0. \quad (3.10) $$

Using $\tilde{\eta} := \eta + \beta > 0$, the control law results in

$$ u = u^{[I]} = x_{1,d}^{(n)} - \sum_{r=0}^{n-2} \alpha_r \xi_1^{(r+1)} - \tilde{\eta} \cdot \mathrm{sign}(s) . \quad (3.11) $$

### 3.2.2. Second-Order Sliding Mode Control.
For a second-order sliding mode, both $s = 0$ and $\dot{s} = 0$ have to be ensured by the feedback controller [Eke10, FL02]. This can be achieved by adding a first-order lag dynamics on the left-hand side of

$$ \gamma_1 \dot{s} + \gamma_0 s = \sum_{r=-1}^{n-1} \alpha_r \xi_1^{(r)} \quad \text{with} \quad \xi_1^{(-1)}(t) := \int_0^t \xi_1(\tau) d\tau. \quad (3.12) $$

For the sake of asymptotic stability, the coefficients $\gamma_0$ and $\gamma_1$ need to be strictly positive. For $\alpha_{-1} = 0$, this sliding surface has a proportional and differentiating characteristic, while $\alpha_{-1} \neq 0$ means that the integral $\xi_1^{(-1)}$ of the output error is additionally fed back in the definition of the sliding variable $s$.

To ensure both $s = 0$ and $\dot{s} = 0$ for the closed-loop control system, the Lyapunov function candidate, which was chosen in the previous subsection as a function solely depending on the value of $s$, is redefined in the form

$$ V^{[II]} = \frac{1}{2} \cdot (s^2 + \lambda \dot{s}^2) \quad \text{with the scaling factor } \lambda > 0. \quad (3.13) $$

This extended specification typically enhances the closed-loop system's robustness by the introduction of the temporal integral of the error signal according to (3.12) in the definition of the sliding surface. In practical applications, $\lambda$ is chosen as shown subsequently so that the stability requirements are simplified.

The computation of the time derivative of (3.13) results in

$$ \dot{V}^{[II]} = s \cdot \dot{s} + \dot{s} \cdot \left( -\frac{\lambda \gamma_0}{\gamma_1} \dot{s} + \frac{\lambda}{\gamma_1} \sum_{r=0}^{n} \alpha_{r-1} \xi_1^{(r)} \right) < 0, \quad (3.14) $$

where the specific parameter choice $\lambda = \gamma_1 > 0$ is used in the following without loss of generality [RS16]. Also according to [RS16], the stabilization of the closed-loop system towards $s = 0$ can be achieved by requiring

$$ \dot{V}^{[II]} < -\eta_1 \cdot |\dot{s}| - \eta_2 \cdot |s| \cdot |\dot{s}| = -\dot{s} \cdot \mathrm{sign}(\dot{s}) \cdot (\eta_1 + \eta_2 \cdot |s|). \quad (3.15) $$
This finally leads to the nonlinear feedback controller [Eke10, generalized form of Eqs. (22), (23)]

$$u = u^{(\text{II})} = x_{1,d}^{(n)} + \frac{1}{\alpha_{n-1}} \cdot \left( \gamma_0 \dot{s} - s - \sum_{r=0}^{n-1} \alpha_{r-1} \xi_1^{(r)} - \text{sign}(\dot{s}) \cdot (\tilde{\eta}_1 + \tilde{\eta}_2 \cdot |s|) \right) \quad (3.16)$$

with $\tilde{\eta}_i \ge \eta_i > 0$ for both $i \in \{1, 2\}$. In analogy to the previous subsection, constant design parameters are the gains $\tilde{\eta}_i$ of both variable-structure control components as well as the coefficients $\alpha_{r-1}$ together with $\gamma_0$.

**3.3. Extension by One- and Two-Sided Barrier Functions.** Both control laws $u^{(\text{I})}$ (first-order sliding mode) and $u^{(\text{II})}$ (second-order sliding mode) can be extended by one- and two-sided barrier functions, which are subsequently distinguished by the superscripts $\langle A \rangle$ and $\langle B \rangle$ [TGT09].

**3.3.1. One-Sided State Constraints.** For the case of a one-sided barrier that represents an upper bound for the output variable, it is necessary that the constraint<fnref n="6" />

$$x_1(t) < \bar{x}_{1,\text{max}}(t) := x_{1,d}(t) + \Delta x_{1,\text{max}} \quad \text{with} \quad \Delta x_{1,\text{max}} > 0 \quad (3.17)$$

is not violated for any point of time $t > 0$. As a prerequisite, we assume that the initial conditions for the state vector $\mathbf{x}(t)$ at $t = 0$ are compatible with this constraint as well. Moreover, it is necessary that the sliding surface $s = 0$ (equivalent to $x_1 = x_{1,d}$) lies within the admissible operating range defined in (3.17). Here, $\bar{x}_{1,\text{max}}$ denotes the maximum admissible upper bound for the state $x_1$, while $\Delta x_{1,\text{max}}$ is a safety margin expressing the minimum difference between the reference signal and the threshold value.

Then, the extended Lyapunov function ansatz

$$V^{(j,A)} = V^{(j)} + V^{(\text{A})} > 0 \quad \text{for } s \neq 0 \quad \text{with} \quad (3.18)$$

$$V^{(\text{A})} = \rho_V \cdot \ln \left( \frac{\sigma_V \cdot \bar{x}_{1,\text{max}}}{\bar{x}_{1,\text{max}} - x_1} \right) \quad \text{and} \quad x_1 < \bar{x}_{1,\text{max}} \quad (3.19)$$

is introduced for both alternatives $j \in \{\text{I}, \text{II}\}$. In (3.19), $\rho_V > 0$ is selected so that the singularity $\bar{x}_{1,\text{max}} - x_1 = 0$ represents a repelling potential and that control constraints are not violated for usual operating conditions. Then, the term $V^{(j)}$ has dominating influence in the neighborhood of $s = 0$, while $\sigma_V > 0$ can be used to adapt the steepness of the barrier function near its singularity.

In the case of the first-order sliding mode, the time derivative of (3.18) can be computed as

$$\dot{V}^{(j,A)} = \dot{V}^{(j)} + \dot{V}^{(\text{A})} < 0 \quad \text{with} \quad \dot{V}^{(\text{A})} = \frac{\rho_V}{\bar{x}_{1,\text{max}}} \cdot \left( \frac{-x_1 \cdot \dot{x}_{1,\text{max}} + \dot{x}_1 \cdot \bar{x}_{1,\text{max}}}{\bar{x}_{1,\text{max}} - x_1} \right), \quad (3.20)$$

where $\dot{V}^{(\text{A})}$ does not explicitly depend on the system input $u$.

In analogy to the fundamental first-order sliding mode control law $u^{(\text{I})}(t)$ derived from (3.10), the inequality

$$s \cdot \left( \sum_{r=0}^{n-2} \alpha_r \xi_1^{(r+1)} + u - x_{1,d}^{(n)} + \eta \cdot \text{sign}(s) + \frac{1}{s} \cdot \dot{V}^{(\text{A})} \right) < 0 \quad (3.21)$$

<footnote n="6">⁶For cases, where constraints on other interval state variables are specified, the reader is referred to [TGT09].</footnote>
needs to be satisfied. Following the same steps as in the derivation of $u^{(I)}(t)$ yields the control law

$$u = u^{(I,A)} = u^{(I)} - \frac{s}{s^2 + \bar{\epsilon}} \cdot \dot{V}^{(A)}, \qquad (3.22)$$

in which $\frac{1}{s}$ in (3.21) has been approximated by $\frac{s}{s^2 + \bar{\epsilon}}$ with the small constant $\bar{\epsilon} > 0$ to ensure regularity of the control law $u^{(I,A)}(t)$ on the sliding surface $s = 0$. Moreover, this modification guarantees that the barrier function becomes inactive as soon as the control goal has been reached.

In analogy to (3.22), also the second-order sliding mode control (3.16) can be extended by the barrier function (3.19). Following the same steps as in Eqs. (3.20)–(3.22) leads to

$$u = u^{(II,A)} = u^{(II)} - \frac{1}{\alpha_{n-1}} \cdot \frac{\dot{s}}{\dot{s}^2 + \bar{\epsilon}} \cdot \dot{V}^{(A)}. \qquad (3.23)$$

### 3.3.2. Two-Sided State Constraints.
For the case of two-sided state constraints, the Lyapunov functions $V^{\langle j \rangle}$, $j \in \{I, II\}$, are extended according to

$$V^{\langle j, B \rangle} = V^{\langle j \rangle} + V^{\langle B \rangle} > 0 \quad \text{for} \quad s \neq 0, \qquad (3.24)$$

where the additive term $V^{\langle B \rangle}$ is chosen in this paper so that state deviations $|x_1(t) - x_{1,d}(t)| \ge \bar{\chi}$ are guaranteed to be avoided [RS16]. Symmetric barriers can be enforced by specifying a barrier function in the form

$$V^{\langle B \rangle} = \rho_V \cdot \ln \left( \frac{\bar{\chi}^{2l}}{\bar{\chi}^{2l} - (x_1 - x_{1,d})^{2l}} \right) \quad \text{with} \quad l \in \mathbb{N}. \qquad (3.25)$$

Large values for $l$ typically lead to the fact that resulting state trajectories come closer to the edges of the admissible operating range. As before, the time derivative of the additive term (3.25) is computed, which yields

$$\dot{V}^{\langle B \rangle} = \rho_V \cdot \frac{2l \cdot (x_1 - x_{1,d})^{2l-1} \cdot (\dot{x}_1 - \dot{x}_{1,d})}{\bar{\chi}^{2l} - (x_1 - x_{1,d})^{2l}}. \qquad (3.26)$$

In full analogy to Eqs. (3.20)–(3.23), the requirement $\dot{V}^{\langle j, B \rangle} < 0$ for $s \neq 0$ (and $\dot{s} \neq 0$, resp.) leads to the extension

$$u^{(I,B)} = u^{(I)} - s^{-1} \cdot \dot{V}^{(B)} \qquad (3.27)$$

of the first-order sliding mode controller and to

$$u^{(II,B)} = u^{(II)} - (\alpha_{n-1} \cdot \dot{s})^{-1} \cdot \dot{V}^{(B)} \qquad (3.28)$$

for the second-order case, where the same regularization strategies for the rational terms $\frac{1}{s}$ and $\frac{1}{\dot{s}}$ become necessary as before. For the implementation of the control law, the term $\dot{s}$ (if required) is typically estimated at run-time by a low-pass filtered differentiation or by means of an observer.
### 3.4. Interval Extensions to Handle Bounded Parameter Uncertainty and State Estimation Errors.
To guarantee asymptotic stability despite bounded uncertainty in both the measured (resp., estimated) states $x(t)$ and system parameters $p$ included in the control law, interval techniques can be applied at run-time during the execution of the previously derived control laws if they are extended according to [RSA15,RS16]. We assume in the following that a suitable interval observer (e.g., [MB10,RE18]) provides the interval bounds $x \in [x] = [\underline{x}; \bar{x}]$ for the state vector or that all components of this vector are directly measured with bounded uncertainty. In addition, it is assumed that the system model is given as an $n$-th order set of ODEs (3.1) in nonlinear controller canonical form. If this is not the case, extended control approaches according to [SRA16] are applicable.

For the sake of controllability (and, hence, also for the existence of the following interval-based variable-structure controllers), it has to be guaranteed that

$$0 \notin [b]([\underline{x}], [\bar{p}]), \quad \forall \underline{x} \in [\underline{x}] \quad \text{and} \quad \forall \bar{p} \in [\bar{p}] \quad (3.29)$$

holds as already discussed for Eq. (3.2). The constraint (3.29) is a necessary controllability requirement that needs to be satisfied along all possible state trajectories $x(t)$ as well as for all possible system parameters. To handle the set-valued state and parameter uncertainty, the output tracking error and its $r$-th derivative are enclosed by the intervals

$$\xi_1^{(r)}(t) \in [\xi_1^{(r)}](t) = [x_1^{(r)}](t) - x_{1,d}^{(r)}(t) \quad (3.30)$$

for each $r \in \{0, 1, \dots, n\}$. These tracking error intervals can be used to generalize the first-order variable-structure controller (without and with state constraints) according to

$$[v^{(I)}] = \frac{-[a]([\underline{x}], [\bar{p}]) + x_{1,d}^{(n)} - \sum_{r=0}^{n-2} \alpha_r \cdot [\xi_1^{(r+1)}] - \tilde{\eta} \cdot \text{sign}([s])}{[b]([\underline{x}], [\bar{p}])}, \quad (3.31)$$

$$[v^{(I,A)}] = [v^{(I)}] - \frac{1}{[b]([\underline{x}], [\bar{p}])} \cdot \frac{[s]}{[s]^2 + \tilde{\epsilon}} \cdot [\dot{V}^{(A)}], \quad (3.32)$$

and

$$[v^{(I,B)}] = [v^{(I)}] - \frac{1}{[b]([\underline{x}], [\bar{p}])} \cdot \frac{[s]}{[s]^2 + \tilde{\epsilon}} \cdot [\dot{V}^{(B)}]. \quad (3.33)$$

Similarly, interval-based generalizations can be defined for all before-mentioned second-order formulations. For a detailed discussion of restrictions in the case that accounts for integrator extensions in the definition of the sliding surface $s = 0$, the reader is referred to [RS16].

In previous work such as [RKA17,RSA16,RSA15,RS16], where this approach was applied successfully to the control of the thermal behavior of high-temperature solid oxide fuel cells as well as for the control of inverted pendulum systems, the C++ toolbox C-XSC [Krä12] was used. Required derivatives, necessary to transform general nonlinear state equations into the system representation assumed in this section, can easily be obtained with the help of algorithmic differentiation, cf. [SRA16]. The template-based library FADBAD++ has shown its efficiency for this purpose, because it cannot only be applied to expressions involving classical point-valued data types but also interval variables.

It has to be pointed out that the actual implementation of the control laws (3.31)–(3.33) is performed in such a way that a point-valued control signal $v$ is chosen from the computed intervals (3.31)–(3.33) so that it guarantees asymptotic stability regardless of the sign of $b(\underline{x}, \bar{p})$. According to [RS16], this is done by testing the negative definiteness of the Lyapunov function candidate for the infima and suprema $\underline{v} := \inf\{[v]\}$ and $\bar{v} := \sup\{[v]\}$. To account for roundoff and representation errors, these values are inflated by a small constant $\epsilon > 0$ to obtain the final set of candidates $\mathcal{V} := \{\underline{v} - \epsilon, \underline{v} + \epsilon, \bar{v} - \epsilon, \bar{v} + \epsilon\}$, from which the signal
with minimum absolute value is chosen that ensures $\dot{V} < 0$ (or more generally $\dot{V}^{(j,A)} < 0$, $\dot{V}^{(j,B)} < 0$, $j \in \{I, II\}$) despite the considered interval uncertainty as a criterion for stability of the uncertain dynamic system model.

## 4. GUARANTEED MODEL PREDICTIVE CONTROL

4.1. **Goal of the Control Synthesis.** Predictive control techniques [CA13] are powerful approaches for controlling uncertain systems, where the corresponding dynamics are formulated in terms of initial value problems for finite-dimensional sets of nonlinear ordinary differential equations. Besides stabilizing the dynamics towards a desired reference trajectory, it is possible to formulate the involved optimization problem in such a way that constraints for the admissible state trajectories and system inputs are handled. Classical approaches for nonlinear model predictive control commonly do not account for interval bounds in initial states and system parameters.

The main drawback of neglecting this information about uncertainty is the lack of a robustness measure of the resulting control. This neglection may lead to scenarios where the determined control signal does only stabilize the system model for parts of the possible parameter ranges, while the resulting dynamics are unstable for others. This is especially critical for open-loop unstable system models such as the inverted pendulum system in Figure 3, for which a model predictive controller was designed in [FR24]. This controller simultaneously aims at a swing-up of the pendulum arm $S_2$ by rotation of the arm $S_1$ and holding $S_2$ in its unstable upright equilibrium by means of pure inertia coupling.

<figure><img src="image_8.png" /><figcaption>FIGURE 3. Application scenario: Swing-up of an inverted pendulum according to [FR24].</figcaption></figure>
## 4.2. Formulation of a Set-Valued Initial Value Problem for Model Predictive Control.
The issue of a missing robustness quantification of classical model predictive controllers is resolved in this section, by presenting a control procedure for interval-valued initial value problems in the form [FR24]

$$ egin{cases} 
\dot{\mathbf{x}}(t) = \mathbf{f}(t, \mathbf{x}(t), \mathbf{u}(t), \mathbf{p}) \\ 
\mathbf{x}_0 \in [\mathbf{x}_0] \subseteq \mathbb{R}^n \\ 
\mathbf{u}(t) \in [\mathbf{u}] \subseteq \mathbb{R}^m \\ 
\mathbf{p} \in [\mathbf{p}] \subseteq \mathbb{R}^p, 
\end{cases} 
\quad (4.1) $$

where the state vector is denoted by $\mathbf{x}(t)$, the vector of dynamic parameters by $\mathbf{p}$, and the control vector by $\mathbf{u}$. The sets $[\mathbf{x}_0] = [[x_{10}] \dots [x_{n0}]]^T$, $[\mathbf{u}] = [[u_1] \dots [u_m]]^T$, and $[\mathbf{p}] = [[p_1] \dots [p_p]]^T$, expressed as interval boxes, are respectively the initial condition of the state vector, the interval-bounded input (assumed without loss of generality to have constant bounds over the time horizon of interest), and the set of feasible dynamic parameters. We assume that this interval-based initial value problem has a unique solution $\mathbf{x}(t, \mathbf{x}_0, \mathbf{u}, \mathbf{p})$ at $t > 0$. A sufficient condition for this is that $\mathbf{f}: \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R}^p \to \mathbb{R}^n$ is continuous in $t$ and Lipschitz in $\mathbf{x}(t)$. For the use of the proposed solution algorithm, we further require that $\mathbf{f}$ is sufficiently smooth, i.e., of class $C^k$. For an approach applicable to linear time-invariant and linear parameter-varying systems, which is similar to the one described subsequently, the reader is referred to [dERP22, dRdSER22a, dRdSER22b].

## 4.3. Interval-Based Nonlinear Model Predictive Control.
The interval-based nonlinear model predictive control approach presented in [FR24, FAdS21] is a generalization of the optimization approach in [RH09] and is based on the computation of interval bounds for the control sequence over a receding horizon which takes into account bounded uncertainties in the parameters of the dynamic system model and in the measured data.

The control intervals are calculated such that the state trajectories converge (monotonically) to the set-point interval over time (i.e., $\mathbf{x}(t_j) \to [\mathbf{x}_r]$, as $j \to \infty$, and all the state and input constraints are satisfied (i.e., $\mathbf{x}(t_j) \in [\mathbf{x}(t_j)]$ and $\mathbf{u}(t_j) \in [\mathbf{u}(t_j)]$, $\forall j$). Here, the set-point interval either represents a bounded domain around a desired operating point (in the pendulum example, the upright position of the arm $S_2$) or a bounded tube around a desired time-dependent trajectory.

The interval algorithm for determining the model predictive control strategy encompasses two stages [FAdS21], namely,

* **Filtering and branching:** This first step provides a sequence of guaranteed input interval boxes at each time-step $t_k$ over the prediction horizon $N_p$, denoted as $[\mathbf{U}(t_k)] = [\mathbf{u}(t_k)] \times [\mathbf{u}(t_{k+1})] \times \dots \times [\mathbf{u}(t_{k+N_p-1})]$. Branching and filtering procedures allow the computation of safe input intervals<fnref n="7" /> along the receding time horizon. In this context, safety means that the state constraints that are specified in the form of admissible intervals in each time step (i.e., $\forall j, [\mathbf{x}(t_j)] \subseteq [\mathbf{x}_{\min}(t_j), \mathbf{x}_{\max}(t_j)]$ are guaranteed to be obeyed despite the parameter uncertainty, where $\mathbf{x}_{\min}(t_j)$ and $\mathbf{x}_{\max}(t_j)$ are the bounds for the admissible domain for each state variable). In addition, the desired convergence to the reference interval (i.e., $[\mathbf{x}_j] \to [\mathbf{x}_r]$) needs to be ensured.

<footnote n="7">⁷For the implementation of this branching algorithm we assume piecewise constant input intervals. Therefore, the more general notation $[\mathbf{u}(t_k)]$ is subsequently replaced by $[\mathbf{u}_k]$ to indicate control signals only changing their value at the points in time $t_k$.</footnote>
* **Interval optimization:** Since safe inputs are computed over a finite time horizon, the optimization algorithm is launched to compute the optimal inputs $[\mathbf{U}(t_k)]^*$ by minimizing as much as possible a newly formulated interval objective function to reduce both the norm of the input intervals (i.e., to reduce the control effort) and the error between the predicted and the reference set-point (or, resp., reference trajectory).

4.4. **Formulation and Minimization of Interval Cost Functions.** To determine feasible intervals of input values that satisfy the constraints on the state trajectories and system inputs, an optimization procedure is required. It aims at finding the optimal control box and — on this basis — a point-valued input $\mathbf{u}(t_k)$ that can be applied to the actuator of the considered system [FR24, FAdS21]. The continuous cost function to be minimized can be expressed over the prediction horizon $T_p = T_c \times N_p$, with the control update step size $T_c$, as

$$J(\mathbf{x}(t), \mathbf{u}(t)) = \int_{t}^{t+T_p} F(\mathbf{x}(\tau), \mathbf{u}(\tau))d\tau . \quad (4.2)$$

Here, $F$ is commonly chosen as a quadratic function that depends on the norm of the inputs and the error between the predicted states (resp, outputs) $\mathbf{x}(t)$ and the reference $\mathbf{x}_r$ according to

$$F(\mathbf{x}(t), \mathbf{u}(\tau)) = (\mathbf{x}(t) - \mathbf{x}_r)^T \mathbf{Q} (\mathbf{x}(t) - \mathbf{x}_r) + \mathbf{u}(t)^T \mathbf{R} \mathbf{u}(t), \quad (4.3)$$

where $\mathbf{Q}$ and $\mathbf{R}$ are both positive (semi-)definite weighting matrices.

Under the assumption of piecewise constant control inputs, the continuous objective function (4.2) takes the form

$$\begin{aligned} J(\mathbf{x}(t), \mathbf{u}(t)) &= \int_{t}^{t+T_p} \left( (\mathbf{x}(\tau) - \mathbf{x}_r)^T \mathbf{Q} (\mathbf{x}(\tau) - \mathbf{x}_r) + \mathbf{u}(\tau)^T \mathbf{R} \mathbf{u}(\tau) \right) d\tau \\ &= \sum_{k=1}^{N_p} \int_{t+(k-1)T_c}^{t+kT_c} \left( (\mathbf{x}(\tau) - \mathbf{x}_r)^T \mathbf{Q} (\mathbf{x}(\tau) - \mathbf{x}_r) \right) d\tau + \sum_{k=1}^{N_p} \int_{t+(k-1)T_c}^{t+kT_c} \mathbf{u}(\tau)^T \mathbf{R} \mathbf{u}(\tau) d\tau . \end{aligned} \quad (4.4)$$

As shown in [FR24, FAdS21], the value of cost function (4.4) can be enclosed with the help of techniques from interval analysis. This is achieved by applying a validated integration method for the initial value problem under consideration. This method provides a list of tight enclosures $[\mathbf{x}_0], [\mathbf{x}_1], \dots, [\mathbf{x}_K]$ for the solution of the model (4.1) on each interval between two subsequent time steps. These enclosures are obtained starting from the initial conditions $[\mathbf{x}_0]$ with bounded, piecewise constant control sequences $[\mathbf{u}_j]$ so that the obtained boxes satisfy the condition $\mathbf{x}(t) \in [\mathbf{x}_t]$ for all $t \in [t + (k-1)T_c; t + kT_c]$. By using these piecewise defined interval enclosures, bounds for the cost function are obtained according to

$$\begin{aligned} J(\mathbf{x}_t, \mathbf{u}_t) &\in T_c \sum_{j=k}^{k+N_p-1} \left( \left( [\mathbf{x}_{[j:j+1]}] - \mathbf{x}_r \right)^T \mathbf{Q} \left( [\mathbf{x}_{[j:j+1]}] - \mathbf{x}_r \right) + [\mathbf{u}_j]^T \mathbf{R} [\mathbf{u}_j] \right) \\ &\le \sup \left\{ T_c \sum_{j=k}^{k+N_p-1} \left( \left( [\mathbf{x}_{[j:j+1]}] - \mathbf{x}_r \right)^T \mathbf{Q} \left( [\mathbf{x}_{[j:j+1]}] - \mathbf{x}_r \right) + [\mathbf{u}_j]^T \mathbf{R} [\mathbf{u}_j] \right) \right\}, \end{aligned} \quad (4.5)$$

where $[\mathbf{x}_{[j:j+1]}]$ is a guaranteed enclosure of all reachable states over the temporal slice $t \in [t_j; t_j + T_c]$.

Figure 4 gives an overview of a possible extension of the interval-based nonlinear model predictive control procedure. This extension is represented by a robust pre-stabilization of
open-loop unstable plants with the help of linear matrix inequality (LMI) techniques. In such cases, the interval-based optimization procedure does not directly provide the overall control sequence but rather a feedforward control action [$u_{ff}$] that allows for optimizing the tracking behavior within the operating domain in which the underlying feedback controller possesses provable asymptotic stability properties [FR24]. For an application of this control procedure to the swing-up control of a rotary inverse pendulum sketched in Sec. 4.1, the reader is referred to [FR24].

<figure><img src="image_9.png" /><figcaption>FIGURE 4. Overall structure of the interval-based nonlinear model predictive control approach, extended by an underlying pre-stabilization of the plant dynamics.</figcaption></figure>

## 5. INTERVAL OBSERVER-BASED IDENTIFICATION OF STATE-DEPENDENT STATIC SUBSYSTEM MODELS

5.1. **Goal of the Estimator Synthesis.** Set-based estimation strategies can be used beyond the case of reconstructing internal system states in a model-based manner by means of a knowledge about control variables and measured outputs. If combined efficiently with contractor-based procedures, it is also possible to exploit the estimates provided by an interval observer to determine bounds for nonlinear static subsystem models that are embedded in the overall system representation. An example for this is the reconstruction of the nonlinear open-circuit voltage characteristic of a Lithium-ion battery on the basis of an electric equivalent circuit model, knowledge about the terminal current as an actuation signal, and the measurement of the terminal voltage with bounded uncertainty [LR24, LR23].
For this estimation task, consider the equivalent circuit model of the battery as described in [EVU09]. It is represented by the quasi-linear state equations

$$ \dot{\mathbf{x}}(t) = \mathbf{A}(\sigma(t)) \cdot \mathbf{x}(t) + \mathbf{b}(\sigma(t)) \cdot u(t) \quad (5.1) $$

$$ = \begin{bmatrix} 0 & 0 & 0 \\ 0 & \frac{-1}{C_{\text{TS}}(\sigma(t)) \cdot R_{\text{TS}}(\sigma(t))} & 0 \\ 0 & 0 & \frac{-1}{C_{\text{TL}}(\sigma(t)) \cdot R_{\text{TL}}(\sigma(t))} \end{bmatrix} \cdot \mathbf{x}(t) + \begin{bmatrix} \frac{-1}{C_{\text{Bat}}} \\ \frac{1}{C_{\text{TS}}(\sigma(t))} \\ \frac{1}{C_{\text{TL}}(\sigma(t))} \end{bmatrix} \cdot u(t) $$

with the state vector

$$ \mathbf{x}(t) = [\sigma(t) \quad v_{\text{TS}}(t) \quad v_{\text{TL}}(t)]^T . \quad (5.2) $$

These state variables denote the state of charge $\sigma$ of the battery as well as two voltages $v_{\text{TS}}$ and $v_{\text{TL}}$ across resistor-capacitor sub-networks included in the equivalent circuit to represent delay phenomena between step-like changes of the terminal current and measurable variations of the terminal voltage. According to [LR24], interval enclosures of the state variables are estimated after reformulation of both the state equations (5.1) and the output equation

$$ \tilde{v}_{\text{OC}}(\sigma(t)) = \mathbf{C}(\sigma(t)) \cdot \mathbf{x}(t) \quad (5.3) $$

$$ = \eta_{\text{OC}}(\sigma(t)) \cdot \sigma(t) = v_{\text{OC}}(\sigma(t)) - v_0 - v_2 \quad (5.4) $$

$$ = \left( v_0 \cdot \frac{e^{v_1 \cdot \sigma(t)} - 1}{\sigma(t)} + v_3 + v_4 \cdot \sigma(t) + v_5 \cdot (\sigma(t))^2 \right) \cdot \sigma(t) \quad (5.5) $$

in a quasi-linear form in which state of charge dependent matrices and vectors $\mathbf{A}(\sigma(t))$, $\mathbf{b}(\sigma(t))$, and $\mathbf{C}(\sigma(t))$, respectively, are separated from the state vector $\mathbf{x}(t)$ and the input variable $u(t)$, which is represented by the terminal current $i_{\text{T}}$. Nonlinear approximations for the state of charge dependent variations of the resistances and capacitances $R_{\text{TS}}(\sigma(t))$, $R_{\text{TL}}(\sigma(t))$, $C_{\text{TS}}(\sigma(t))$, and $C_{\text{TL}}(\sigma(t))$ for processes with short and long time constants (indices TS and TL, respectively) and nominal values for the coefficients $v_0, \dots, v_5$ of the open-circuit voltage $v_{\text{OC}}(\sigma(t))$ are identified beforehand by means of experimental data according to [RMAR16]. A conservative set-based enclosure of the resulting matrices then forms the basis for the subsequent observer interval observer design.

## 5.2. Interval Observer Synthesis.
For the design of an interval observer, cooperativity of the system model (as discussed in the introduction of this article) as well as stability of the estimation error dynamics need to be ensured. This can be achieved by choosing the observer gain matrix $\mathbf{H}$ according to [HKRA20], in the form

$$ \mathbf{H} = [ \begin{bmatrix} h_1 & 0 & 0 \end{bmatrix} ]^T \quad (5.6) $$

with the constant value $h_1 > 0$.

In such a way, interval bounds $\hat{\mathbf{x}} \in [\hat{\mathbf{x}}] := [\hat{\mathbf{x}}; \hat{\mathbf{x}}]$ can be determined for the true, non-measurable state vector $\mathbf{x}(t)$ defined in (5.2) by means of lower and upper bounding systems [HKRA20, RE18], following Müller's theorem [Mül27]. These bounding systems are given by

$$ \underline{\mathbf{A}}_{\text{O}}\hat{\mathbf{x}} + \underline{\mathbf{b}}\mathbf{u} + \underline{\mathbf{H}}\underline{\mathbf{y}}_{\text{m}} \le \dot{\hat{\mathbf{x}}} \le \overline{\underline{\mathbf{A}}_{\text{O}}\hat{\mathbf{x}}} + \overline{\underline{\mathbf{b}}\mathbf{u}} + \underline{\mathbf{H}}\overline{\underline{\mathbf{y}}_{\text{m}}}, \quad \hat{\mathbf{x}} \in [\hat{\mathbf{x}}], \quad (5.7) $$

where

$$ \underline{\mathbf{A}}_{\text{O}} = \mathbf{A} - \mathbf{H}\mathbf{C} \in [\underline{\underline{\mathbf{A}}}_{\text{O}}; \overline{\underline{\underline{\mathbf{A}}}_{\text{O}}}] \quad (5.8) $$
<figure><img src="image_10.png" /><figcaption>FIGURE 5. Estimation of the open-circuit voltage and the state of charge.</figcaption></figure>

is a Metzler matrix [GT08, KRA19, Kac02] with non-negative off-diagonal elements and

$$[\mathbf{y}_m] := [\underline{\mathbf{y}}_m ; \bar{\mathbf{y}}_m] = \mathbf{y}_m + [-\Delta\mathbf{y}_m ; \Delta\mathbf{y}_m] \quad (5.9)$$

is the measured system output with bounded uncertainty. To remove the structural limitation for the gain matrix (5.6) in which two of the elements need to be strictly equal to zero, an extended as well as cascaded observer design was recently published in [LRD24]. This extension provides additional tuning parameters for the observer design so that the estimation uncertainty is typically reduced by a reduction of the sensitivity of the state bounds against model inaccuracies and measurement errors.

The outcome of the observer-based state reconstruction is used in the following subsection for an interval-based set intersection approach. This procedure allows for identifying the characteristic of the open circuit voltage $v_{OC}(\sigma(t))$ which was structurally unknown to the interval observer and only considered by an outer interval enclosure for the matrix $\mathbf{C}(\sigma(t))$.

### 5.3. Implementation for the Estimation of the Open-Circuit Voltage Characteristic of a Lithium-Ion Battery.
According to the equivalent circuit representation of the dynamics of a Lithium-ion battery, an interval-based estimate can be determined for the open-circuit voltage according to

$$\tilde{v}_{OC}(t) = y_{m}^{*}(t) + \hat{v}_{TS}(t) + \hat{v}_{TL}(t), \quad (5.10)$$

with $\hat{v}_{TS}(t) \in [\hat{v}_{TS}(t)]$ and $\hat{v}_{TL}(t) \in [\hat{v}_{TL}(t)]$. Together with the interval enclosure $[\hat{\sigma}]$ for the state of charge, an interval box $[\Gamma] = [\sigma(t)] \times [\tilde{v}_{OC}(t)]$, $\sigma(t) \in [\underline{\sigma}(t) ; \bar{\sigma}(t)]$ and $\tilde{v}_{OC}(t) \in [\tilde{v}_{OC}(t) ; \tilde{\tilde{v}}_{OC}(t)]$, can be constructed. All required computational steps are summarized in the block diagram depicted in Fig. 5.
Successive intersections of these interval boxes $[\Gamma]$ — at different points in time — are utilized to improve the approximation of the true $v_{OC}(\sigma(t))$ characteristic according to Figure 6.

<figure><img src="image_11.png" /><figcaption>FIGURE 6. Identification of nonlinear dependencies using interval methods.</figcaption></figure>

For strategies, allowing for a reduction of the computational complexity due to a large number of interval boxes to be intersected, the interval merging routine according to [KRK$^+$10] is employed. Moreover, refinements of the estimated open-circuit voltage characteristic are possible by implementing a forward-backward contractor, as discussed from a methodological point of view in [JKDW01]. These refinements have been demonstrated recently in [LR23].

## 6. COMBINATION OF SET-BASED AND STOCHASTIC UNCERTAINTY REPRESENTATIONS: AN APPLICATION TO ITERATIVE-LEARNING OBSERVER DESIGN

6.1. **Goal of the Estimator Synthesis.** Many practical applications are characterized by the fact that uncertainties and disturbances are neither characterized purely by set-based nor by stochastic representations. Then, it is desired to have a unified representation of uncertainty that can be employed to implement state estimators. In addition, many real-life applications are also subject to nonlinearities, making an exact solution of the estimation task impossible. Widely used approximation techniques such as the Extended or Unscented Kalman Filter [JUDW00, Ste94, AM05] have the disadvantage that they cannot take into consideration the worst-case influence of linearization errors of both nonlinear state equations and measurement models.

Therefore, this section introduces a combination of stochastic and set-based uncertainty representations for the estimator design described in Sec. 6.4. This approach allows, on the one hand, for a rigorous quantification of predefined confidence levels in stochastic state estimation procedures. On the other hand, it allows for handling nonlinearities in such a way that the previously mentioned tolerance bounds are definitely not determined in an overly optimistic manner [RCD$^+$22b, RCD$^+$22a, RWH$^+$22]. Inspired by normally distributed uncertainty models, the suggested approach is based on the use of an ellipsoidal calculus.

Practical application scenarios are, for example, the state estimation of marine vessels on the basis of GPS position measurements [RWH$^+$22], state estimation for battery systems [RJ21a], or state estimators that iteratively enhance the estimation accuracy for system models that are characterized by periodically recurring trajectories [RCD$^+$22b]. An example
for the latter is given in [RCD$^+$22a] which demonstrates the successful application of the observer for the state and disturbance estimation of a Lithium-ion battery with periodically recurring input current profiles.

## 6.2. Preliminaries for Ellipsoidal State Enclosures.
In the following, we focus on an estimator synthesis in a mixed stochastic and set-based uncertainty representation for the case of periodically recurring trajectories and disturbance profiles (e.g., due to errors caused by a deterministic model mismatch due to purposefully performed model simplifications). For that purpose, assume the discrete time system model

$$ \mathbf{z}_{k+1} = \mathbf{\Phi}(\mathbf{z}_k, \mathbf{p}) \cdot \mathbf{z}_k \qquad (6.1) $$

with the bounded parameter vector $\mathbf{p} \in [\underline{\mathbf{p}}; \bar{\mathbf{p}}]$. As discussed in the following, the vector $\mathbf{z}_k \in \mathbb{R}^{\tilde{n}}$ consists of both system states and noise variables. Furthermore, assume that the ellipsoidal domain

$$ \mathcal{E}_k (\boldsymbol{\mu}_k, \boldsymbol{\Gamma}'_k, r) := \{ \mathbf{z}_k \in \mathbb{R}^n \mid (\mathbf{z}_k - \boldsymbol{\mu}_k)^T \boldsymbol{\Gamma}'_k^{-T} \boldsymbol{\Gamma}'_k^{-1} (\mathbf{z}_k - \boldsymbol{\mu}_k) \leq r^2 \} \quad (6.2) $$

specifies the confidence bound of a given percentage if the vector $\mathbf{z}_k$ is normally distributed with the ellipsoid midpoint as the expected value $\boldsymbol{\mu}_k \in \mathbb{R}^{\tilde{n}}$ and the positive definite covariance matrix $\mathbf{Q}'_k = \boldsymbol{\Gamma}'_k \boldsymbol{\Gamma}'_k^T > 0$. In (6.2), the parameter $r$ describes a magnification factor according to [WSM15] associated with the before-mentioned confidence level.

## 6.3. Ellipsoidal Calculus for the Guaranteed Production of Covariance Ellipsoids.
For the implementation of an ellipsoidal covariance prediction step (based on the references [RBJ21, RBJK21, RJ21a]), we use the definitions

$$ \boldsymbol{\Gamma}_k := r \cdot \boldsymbol{\Gamma}'_k \quad \text{and} \quad \mathbf{Q}_k := r^2 \cdot \mathbf{Q}'_k \qquad (6.3) $$

transferring the confidence ellipsoid (6.2) into the normalized standard form

$$ \mathcal{E}_k^{\boldsymbol{\Gamma}} (\boldsymbol{\mu}_k, \boldsymbol{\Gamma}_k) := \{ \mathbf{z}_k \in \mathbb{R}^n \mid (\mathbf{z}_k - \boldsymbol{\mu}_k)^T \boldsymbol{\Gamma}_k^{-T} \boldsymbol{\Gamma}_k^{-1} (\mathbf{z}_k - \boldsymbol{\mu}_k) \leq 1 \} . \qquad (6.4) $$

For the compactness of notation, reformulate the system model (6.1) into the form

$$ \mathbf{z}_{k+1} = \mathbf{\Phi}(\mathbf{z}_k, \mathbf{p}) \cdot \tilde{\mathbf{z}}_k + \tilde{\mathbf{\Phi}} \cdot \boldsymbol{\mu}_k + (\mathbf{\Phi}(\mathbf{z}_k, \mathbf{p}) - \tilde{\mathbf{\Phi}}) \cdot \boldsymbol{\mu}_k \qquad (6.5) $$

with $\mathbf{z}_k = \tilde{\mathbf{z}}_k + \boldsymbol{\mu}_k$, where

$$ \mathbf{z}_k \in \mathcal{E}_k = \mathcal{E}_k (\boldsymbol{\mu}_k, \boldsymbol{\Gamma}'_k, r) . \qquad (6.6) $$

Here, $\mathcal{E}_k$ denotes the uncertainty on the non-origin centered states $\mathbf{z}_k$,

$$ \tilde{\mathbf{z}}_k \in \tilde{\mathcal{E}}_k = \tilde{\mathcal{E}}_k (\mathbf{0}, \boldsymbol{\Gamma}'_k, r) \qquad (6.7) $$

the uncertainty on $\tilde{\mathbf{z}}_k$ after shifting the ellipsoid to the origin, and

$$ \tilde{\mathbf{\Phi}} = \mathbf{\Phi}(\boldsymbol{\mu}_k, \text{mid}([\mathbf{p}])) \qquad (6.8) $$

is the midpoint approximation of the quasi-linear system matrix with

$$ \text{mid}([\mathbf{p}]) = \frac{1}{2} \cdot (\underline{\mathbf{p}} + \bar{\mathbf{p}}) . \qquad (6.9) $$

Let also $\mathcal{E}_k$ denote an axis-aligned enclosure of the ellipsoid $\mathcal{E}_k$ defined in (6.2) in the form of an $\tilde{n}$-dimensional interval box. An alternative to the definition (6.8) is given by

$$ \tilde{\mathbf{\Phi}} = \text{mid}(\mathbf{\Phi}(\mathcal{E}_k, [\mathbf{p}])) , \qquad (6.10) $$
which is preferable if (6.8) and (6.10) strongly differ from each other in the case of large uncertainty.

Then, a confidence bound of the predicted states $z_{k+1}$ corresponding to a given confidence level with the associated magnification $r$ leads to the ellipsoid $\mathcal{E}_{k+1}(\mu_{k+1}, \Gamma'_{k+1}, r)$ with the covariance $Q'_{k+1} = \Gamma'_{k+1}(\Gamma'_{k+1})^T$ for which the parameters are computed in the following steps. For further details of the following algorithm, Consisting of the prediction steps **P1–P3**, the reader is referred to the reference [RCD$^+$22b].

**P1**: Apply

$$ \check{z}_{k+1} = \mathbf{\Phi}(z_k, \mathbf{p}) \cdot \check{z}_k \qquad (6.11) $$

to the ellipsoid $\mathcal{E}_k$ in (6.7). The outer ellipsoid enclosure of the mapping (6.11) 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, \qquad (6.12) $$

where $\alpha_{k+1} \ge 0$ is the smallest value for which the LMI

$$ \mathcal{M}_{k+1} := \mathbf{\Lambda} \begin{bmatrix} -\mathbf{Q}_k^{-1} & \mathbf{\Phi}^T(z_k, \mathbf{p}) \cdot \tilde{\mathbf{\Phi}}^{-T} \\ \tilde{\mathbf{\Phi}}^{-1} \cdot \mathbf{\Phi}(z_k, \mathbf{p}) & -\alpha_{k+1}^2 \mathcal{R}_k \end{bmatrix} \mathbf{\Lambda} \preceq 0 \qquad (6.13) $$

is satisfied for all $z_k \in \mathcal{C}_k$ and $\mathbf{p} \in [\mathbf{p}]$ with the choice

$$ \mathcal{R}_k := \mathbf{\Gamma}_k \cdot \mathbf{\Gamma}_k^T. \qquad (6.14) $$

Alternatives for the choice of $\mathcal{R}_k$ as well as advantages and disadvantages are discussed in [RBJ21]. In (6.13), the symbol $\preceq$ denotes the negative semi-definiteness of the corresponding matrix expression and $\mathbf{\Lambda}$ is a preconditioning matrix chosen according to [RBJK21]. A typical choice for the preconditioning is a block diagonal matrix $\mathbf{\Lambda} = \text{blkdiag}(\beta\mathbf{I}, \beta^{-1}\mathbf{I})$ with $\mathbf{I}$ as the identity matrix and the square root $\beta = \sqrt{\min\{\lambda_i(\mathbf{Q}_k)\}}$ of the smallest eigenvalue of $\mathbf{Q}_k$.

**P2**: Compute interval bounds for the term

$$ \mathbf{b}_k = (\mathbf{\Phi}(z_k, \mathbf{p}) - \tilde{\mathbf{\Phi}}) \cdot \mu_k \in [\mathbf{b}_k] \qquad (6.15) $$

which accounts for a non-zero ellipsoid midpoint with $z_k$, $\tilde{\mathbf{\Phi}}$, and $\mathbf{p}$ defined according to (6.6), (6.8), and (6.9). Inflate the ellipsoid bound described by the shape matrix (6.12) according to [RJ21a]

$$ \mathbf{Q}_{k+1} = (1 + \rho_{\text{O},k+1})^2 \cdot \check{\mathbf{Q}}_{k+1}, \qquad (6.16) $$

$$ \rho_{\text{O},k+1} = \sup \left\{ \|\alpha_{k+1}^{-1} \cdot \mathbf{\Gamma}_k^{-1} \cdot [\mathbf{b}_k]\| \right\}. \qquad (6.17) $$

**P3**: Compute the updated ellipsoid midpoint

$$ \mu_{k+1} = \tilde{\mathbf{\Phi}} \cdot \mu_k \qquad (6.18) $$

and its factorized shape matrix

$$ \Gamma'_{k+1} = \alpha_{k+1} \cdot (1 + \rho_{\text{O},k+1}) \cdot \tilde{\mathbf{\Phi}} \cdot \Gamma'_{k}. \qquad (6.19) $$
## 6.4. Stochastic Iterative-Learning Observer Design for Quasi-Linear State Equations.
To apply the previous set-based prediction step, initialized at $k = 1$ with a suitable ellipsoidal state enclosure for forecasting a certain confidence set, we consider the quasi-linear discrete-time state-space representation

$$ 
\begin{aligned}
\mathbf{x}_{k+1} &= \mathbf{A}(\mathbf{x}_k, \mathbf{p}) \cdot \mathbf{x}_k + \mathbf{E}(\mathbf{x}_k, \mathbf{p}) \cdot \mathbf{w}_k \\
\mathbf{y}_k &= \mathbf{C}(\mathbf{x}_k, \mathbf{p}) \cdot \mathbf{x}_k + \mathbf{v}_k
\end{aligned}
\qquad (6.20) $$

with the state vector $\mathbf{x}_k \in \mathbb{R}^n$, the measured output vector $\mathbf{y}_k \in \mathbb{R}^m$ ($m \le n$), as well as the uncorrelated process and measurement noise vectors $\mathbf{w}_k \in \mathbb{R}^{n_w}$ and $\mathbf{v}_k \in \mathbb{R}^{n_v}$, respectively. Let both noise vectors be normally distributed with the covariances $\mathbf{C}_{\mathbf{w},k}$ and $\mathbf{C}_{\mathbf{v},k}$ and vanishing mean. The iterative-learning observer is designed in this section to estimate the state vector $\mathbf{x}_k$ as well as its uncertainty (expressed by its covariance) by a Kalman Filter-like procedure. This estimator does not only operate along the time domain $k$ but also enhances the estimates successively from the trial $i$ to the trial $i + 1$ during which the underlying system operates along the same trajectories. Here, the word trial means that the system repeats the same finite-duration control task multiple times and that the estimator is meant to enhance its accuracy each time the task is re-executed. To further allow for the identification of a systematic model mismatch, a lumped correction term $\delta_k$ is added to the state equations (6.20) in the form

$$ \mathbf{x}_{k+1} = \mathbf{A}(\mathbf{x}_k, \mathbf{p}) \cdot \mathbf{x}_k + \mathbf{E}(\mathbf{x}_k, \mathbf{p}) \cdot \mathbf{w}_k + \delta_k . \qquad (6.21) $$

Strategies for learning $\delta_k$ during the execution of the iterative-learning estimator are described in [RCD+22b].

For the derivation of the learning-type framework, consider two subsequent trials $\xi = i$ and $\xi = i + 1$ with the actually measured data $\mathbf{y}_{m,k}^{\xi} = \mathbf{C}_k^{\xi} \cdot \mathbf{x}_k^{\xi} + \mathbf{v}_k^{\xi}$ corresponding to the realizations of the general outputs $\mathbf{y}_k$ in (6.20) for these two trials.

For the recursive formulation of the iterative-learning observer, we assume that both process and measurement noise are uncorrelated and normally distributed with zero mean. The superscript $p$ denotes the result of the prediction step, while the superscript $e$ refers to the estimation result as the outcome of the measurement-based innovation step.

### 6.4.1. Prediction Step.
To perform a state prediction in the case of state- and parameter-dependent disturbance inputs $\mathbf{E}(\mathbf{x}_k, \mathbf{p})$ according to [RCD+22a], define an ellipsoid for the augmented state vector

$$ \mathbf{z}_k = \left[ (\mathbf{x}_k^i)^T \quad (\mathbf{x}_k^{i+1})^T \quad (\mathbf{w}_k^i)^T \quad (\mathbf{w}_k^{i+1})^T \right]^T \qquad (6.22) $$

corresponding to the result of the preceding innovation step, augmented by the influence of the uncorrelated process noise in both trials $i$ and $i + 1$ in the form

$$ \mathcal{E}_k^{e,i|i+1} \left( \begin{bmatrix} \mu_k^{e,i} \\ \mu_k^{e,i+1} \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix}, \Gamma_k^e, r \right) \qquad (6.23) $$
with

$$ \Gamma_k^e = \begin{bmatrix} \Gamma_k^{e,i|i+1} & \mathbf{0}_{2n \times 2n_w} \\ \mathbf{0}_{2n_w \times 2n} & \begin{bmatrix} \mathbf{C}_{w,k} & \mathbf{0} \\ \mathbf{0} & \mathbf{C}_{w,k} \end{bmatrix}^{\frac{1}{2}} \end{bmatrix} \quad (6.24) $$

containing the matrix square root $\Gamma_k^{e,i|i+1}$ of the combined state covariance matrix

$$ \mathbf{C}_k^{e,i|i+1} = \Gamma_k^{e,i|i+1} \cdot (\Gamma_k^{e,i|i+1})^T, \quad (6.25) $$

in which joint distributions over the two subsequent trials are denoted by the superscript $i|i+1$, as well as of the trial-independent noise covariance matrix $\mathbf{C}_{w,k}$, and the magnification factor $r \ge 1$ as a user-defined degree of freedom. Here, in analogy to the Kalman Filter, the term innovation step refers to the operation in which the predicted state information obtained from an open-loop evaluation of the system model into the future is corrected by taking into account measured data.

Then, the application of the system model with the augmented dynamics matrix

$$ \Phi(\mathbf{z}_k, \mathbf{p}) = \begin{bmatrix} \mathbf{A}(\mathbf{x}_k^i, \mathbf{p}) & \mathbf{0}_{n \times n} & \mathbf{E}(\mathbf{x}_k^i, \mathbf{p}) & \mathbf{0}_{n \times n_w} \\ \mathbf{0}_{n \times n} & \mathbf{A}(\mathbf{x}_k^{i+1}, \mathbf{p}) & \mathbf{0}_{n \times n_w} & \mathbf{E}(\mathbf{x}_k^{i+1}, \mathbf{p}) \\ \mathbf{0}_{n_w \times n} & \mathbf{0}_{n_w \times n} & \mathbf{I}_{n_w \times n_w} & \mathbf{0}_{n_w \times n_w} \\ \mathbf{0}_{n_w \times n} & \mathbf{0}_{n_w \times n} & \mathbf{0}_{n_w \times n_w} & \mathbf{I}_{n_w \times n_w} \end{bmatrix} \quad (6.26) $$

to the ellipsoid (6.23) according to Sec. 6.3 yields the ellipsoid

$$ \mathcal{E}_{k+1}^{p,i|i+1} \left( \begin{bmatrix} \boldsymbol{\mu}_{k+1}^{p,i} \\ \boldsymbol{\mu}_{k+1}^{p,i+1} \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix} \right), \Gamma_{k+1}^{p}, r \quad (6.27) $$

In this ellipsoid, the first $2n$ components $\left[ (\boldsymbol{\mu}_{k+1}^{p,i})^T \quad (\boldsymbol{\mu}_{k+1}^{p,i+1})^T \right]^T$ of the midpoint vector represent the predicted expected value vector and the covariance $\mathbf{C}_{k+1}^{p,i|i+1}$ is obtained by extracting the upper left $(2n \times 2n)$ block of the matrix product $\Gamma_{k+1}^p \cdot (\Gamma_{k+1}^p)^T$.

### 6.4.2. Innovation Step.
The measurement-based innovation step makes use of the deviations

$$ \Delta \mathbf{y}_k^i = \mathbf{y}_{m,k}^i - \mathbf{C}_k^i \cdot \boldsymbol{\mu}_k^{p,i} \quad (6.28) $$

and

$$ \Delta \mathbf{y}_k^{i+1} = \mathbf{y}_{m,k}^{i+1} - \mathbf{C}_k^{i+1} \cdot \boldsymbol{\mu}_k^{p,i+1} \quad (6.29) $$

between the measured data in the trials $i$ and $i+1$, respectively, and the corresponding output forecasts based on the prediction step of the previous subsection. Using these output deviations, and assuming the standard detectability requirements of the Kalman Filter design to be satisfied, the expected values are updated according to [RCD$^+$22b] by

$$ \begin{bmatrix} \boldsymbol{\mu}_k^{e,i} \\ \boldsymbol{\mu}_k^{e,i+1} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\mu}_k^{p,i} \\ \boldsymbol{\mu}_k^{p,i+1} \end{bmatrix} + \tilde{\mathbf{H}}_k \cdot \begin{bmatrix} \mathbf{y}_{m,k}^i \\ \mathbf{y}_{m,k}^{i+1} \end{bmatrix} - \tilde{\mathbf{H}}_k \tilde{\mathbf{C}}_k \cdot \begin{bmatrix} \boldsymbol{\mu}_k^{p,i} \\ \boldsymbol{\mu}_k^{p,i+1} \end{bmatrix}, \quad (6.30) $$
where the combined output matrix

$$ \tilde{C}_k := \begin{bmatrix} C_k^i & \mathbf{0} \\ \mathbf{0} & C_k^{i+1} \end{bmatrix} \qquad (6.31) $$

results from pointwise evaluations of the quasi-linear system's output matrix according to

$$ C_k^i := C_k (\boldsymbol{\mu}_k^{\text{p},i}, \text{mid}([\mathbf{p}])) \qquad (6.32) $$

and

$$ C_k^{i+1} := C_k (\boldsymbol{\mu}_k^{\text{p},i+1}, \text{mid}([\mathbf{p}])) . \qquad (6.33) $$

To optimize the estimator performance, the filter gain matrix $\tilde{\mathbf{H}}_k$ needs to be tuned. This is possible by a minimization of the uncertainty on the estimated state variables. Inspired by the derivation of the Kalman Filter [Kal60], a minimization of the estimation error covariance is performed. For that purpose, this covariance is computed according to

$$ \begin{aligned} C_k^{e,i|i+1} &= \mathrm{E} \left\{ \begin{bmatrix} \mathbf{x}_k^i - \boldsymbol{\mu}_k^{e,i} \\ \mathbf{x}_k^{i+1} - \boldsymbol{\mu}_k^{e,i+1} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{x}_k^i - \boldsymbol{\mu}_k^{e,i} \\ \mathbf{x}_k^{i+1} - \boldsymbol{\mu}_k^{e,i+1} \end{bmatrix}^T \right\} \\ &= \mathrm{Cov} \left\{ \begin{bmatrix} \mathbf{x}_k^i - \boldsymbol{\mu}_k^{e,i} \\ \mathbf{x}_k^{i+1} - \boldsymbol{\mu}_k^{e,i+1} \end{bmatrix} \right\} \\ &= \mathbf{M}_k C_k^{\text{p},i|i+1} \mathbf{M}_k^T + \tilde{\mathbf{H}}_k \tilde{\mathbf{C}}_{v,k} \tilde{\mathbf{H}}_k^T, \end{aligned} \qquad (6.34) $$

where

$$ \mathbf{M}_k = \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} - \tilde{\mathbf{H}}_k \tilde{\mathbf{C}}_k . \qquad (6.35) $$

To ensure causality, i.e., the update in iteration $i$ does not depend on the future iteration $i+1$, the gain matrix is structured according to

$$ \tilde{\mathbf{H}}_k := \begin{bmatrix} \mathbf{H}_{1,k}^{i+1} & \mathbf{0} \\ \mathbf{H}_{2,k}^{i+1} & \mathbf{H}_{1,k}^{i+1} - \mathbf{H}_{2,k}^{i+1} \end{bmatrix} . \qquad (6.36) $$

As shown in [RCD+22b], the minimization task described above can be solved by differentiating the trace of the covariance $C_k^{e,i|i+1}$ with respect to the gain matrix and setting the resulting expression, being linear in both gain matrices $\mathbf{H}_{1,k}^{i+1}$ and $\mathbf{H}_{2,k}^{i+1}$, to zero. A relevant intermediate step is to determine the predicted covariance that is partitioned in a blockwise manner according to

$$ C_k^{\text{p},i|i+1} = \begin{bmatrix} \mathbf{C}_{A,k}^{\text{p}} & \mathbf{C}_{B,k}^{\text{p}} \\ (\mathbf{C}_{B,k}^{\text{p}})^T & \mathbf{C}_{C,k}^{\text{p}} \end{bmatrix} . \qquad (6.37) $$

Moreover, the residual covariance is defined as

$$ \mathbf{P} \cdot \tilde{\mathbf{\Gamma}}_k^e \cdot (\tilde{\mathbf{\Gamma}}_k^e)^T \cdot \mathbf{P}^T + \tilde{\mathbf{C}}_{v,k} = \begin{bmatrix} \mathbf{C}_{A,k} & \mathbf{C}_{B,k} \\ \mathbf{C}_{B,k}^T & \mathbf{C}_{C,k} \end{bmatrix} \qquad (6.38) $$

with the projection matrix $\mathbf{P} = [\mathbf{I}_{m \times m} \quad \mathbf{0}_{m \times (n-m)}]$, $\tilde{\mathbf{\Gamma}}_k^e$ defined in the following, and the trial-independent measurement noise covariance matrix

$$ \tilde{\mathbf{C}}_{v,k} := \begin{bmatrix} \mathbf{C}_{v,k} & \mathbf{0} \\ \mathbf{0} & \mathbf{C}_{v,k} \end{bmatrix} . \qquad (6.39) $$
As a specific feature introduced in [RCD<sup>+</sup>22a], aiming at a robustification of the innovation stage against nonlinearities, the matrix $\tilde{\Gamma}_k^e$ is obtained by the consideration of the quasi-linearity of the output equation with the help of the ellipsoid

$$ \mathcal{E}_k^{\text{p},i|i+1} \left( \left[ \begin{array}{c} \boldsymbol{\mu}_k^{\text{p},i} \\ \boldsymbol{\mu}_k^{\text{p},i+1} \end{array} \right], \Gamma_k^{\text{p},i|i+1}, r \right) \qquad (6.40) $$

that is propagated through a quasi-linear system model in the form (6.1) with the associated system matrix ($\mathbf{y}_{m,k}^i, \mathbf{y}_{m,k}^{i+1} \in \mathbb{R}^m$)

$$ \mathbf{\Phi}(\mathbf{z}_k, \mathbf{p}) = \begin{bmatrix} \mathbf{C}(\mathbf{x}_k^i, \mathbf{p}) & \mathbf{0}_{m \times n} \\ \mathbf{0}_{m \times n} & \mathbf{C}(\mathbf{x}_k^{i+1}, \mathbf{p}) \\ \mathbf{0}_{(n-m) \times m} & \mathbf{I}_{(n-m) \times (n-m)} \\ \mathbf{0}_{(n-m) \times n} & \mathbf{0}_{(n-m) \times m} \end{bmatrix} \mathbf{I}_{(n-m) \times (n-m)}, \quad (6.41) $$

where $\mathbf{\Phi}(\mathbf{z}_k, \mathbf{p}) \in \mathbb{R}^{2n \times 2n}$. The evaluation of this quasi-linear model yields an ellipsoid

$$ \mathcal{E}_k^e (\tilde{\boldsymbol{\mu}}_k^e, \tilde{\boldsymbol{\Gamma}}_k^e, r) \qquad (6.42) $$

from which the associated shape matrix $\tilde{\boldsymbol{\Gamma}}_k^e \cdot (\tilde{\boldsymbol{\Gamma}}_k^e)^T$ is computed in Eq. (6.38), followed by an extraction of the upper left $(2m \times 2m)$ block that corresponds to the actually measured output quantities. This subblock extraction is performed by the multiplication with the projection matrix $\mathbf{P}$ in the first summand of Eq. (6.38).

As shown in further detail in [RCD<sup>+</sup>22b], the optimal iterative-learning observer gains are given by

$$ \left[ \mathbf{H}_{1,k}^{i+1} \mid \mathbf{H}_{2,k}^{i+1} \right] = \left[ \left( \mathbf{C}_k^i \mathbf{C}_{A,k}^p + \mathbf{C}_k^{i+1} \mathbf{C}_{C,k}^p \right)^T \right] \left( \left( \mathbf{C}_k^i \mathbf{C}_{B,k}^p - \mathbf{C}_k^{i+1} \mathbf{C}_{C,k}^p \right)^T \right)^{-1} \cdot \left[ \begin{array}{cc} \mathbf{C}_{A,k} + \mathbf{C}_{C,k} & \star \\ \mathbf{C}_{B,k} - \mathbf{C}_{C,k} & \mathbf{C}_{A,k} - (\mathbf{C}_{B,k} + \mathbf{C}_{B,k}^T) + \mathbf{C}_{C,k} \end{array} \right]^{-1} \qquad (6.43) $$

It should be noted at this point that setting all coupling terms between the iterations $i$ and $i+1$ to zero leads exactly to the same gain matrix that would be obtained when designing a Kalman filter if all matrices in (6.21) are state-independent and only depend on precisely known parameters.

## 7. CONCLUSIONS AND FUTURE WORK

In this paper, a comprehensive overview of recent advances concerning the offline and online use of interval and other set-based approaches for control and state estimation has been given. We have equally addressed methodological approaches and selected fields of practical applications in our selection, where we have put a focus on the tasks of parameter estimation, robust and model predictive control, as well as state estimation procedures.

Future work will focus especially on widening the field of applications by an investigation of control and state estimation tasks for distributed and large-scale, interconnected systems with parameter uncertainty and disturbances that are, among others, omnipresent in control scenarios in the Industry 4.0 as well as in the domain of energy systems. To cope with such systems, we will focus especially on the development of decentralized techniques for the parameter identification of large-scale system models, which will again form the basis
for control and state estimator design. Moreover, we will further develop approaches that
allow for removing restrictive monotonicity assumptions as already shown in works such
as [WLS18, KY23, KSY22].

REFERENCES

[AC16] J. Alexandre dit Sandretto and A. Chapoutot. Validated explicit and implicit Runge-Kutta methods. *Reliable Computing*, 22(1):79-103, Jul 2016.

[AdSC16] J. Alexandre dit Sandretto and A. Chapoutot. DynIBEX: A differential constraint library for studying dynamical systems. In *Conference on Hybrid Systems: Computation and Control (HSCC 2016)*, 2016.

[AK11] M. Althoff and B.H. Krogh. Zonotope bundles for the efficient computation of reachable sets. In *2011 50th IEEE Conference on Decision and Control and European Control Conference*, pages 6814–6821, 2011. doi:10.1109/CDC.2011.6160872.

[AK12] F. Antritter and M. Kletting. Verified online state and parameter estimation for a nonlinear tank system. *IFAC Proceedings Volumes*, 45(16):452–457, 2012. 16th IFAC Symposium on System Identification. doi:10.3182/20120711-3-BE-2027.00360.

[Alt15] M. Althoff. An introduction to CORA 2015. In G. Frehse and M. Althoff, editors, *ARCH14-15. 1st and 2nd International Workshop on Applied veRification for Continuous and Hybrid Systems*, volume 34 of *EPiC Series in Computing*, pages 120–151. EasyChair, 2015. https://easychair.org/publications/paper/xMm, accessed: Sept. 1, 2023. doi:10.29007/zbkv.

[AM05] B.D.O. Anderson and J.B. Moore. *Optimal Filtering*. Dover Publications, Inc., Mineola, USA, 2005.

[AR12] E. Auer and A. Rauh. VERICOMP: A system to compare and assess verified IVP solvers. *Computing*, 94(2-4):163-172, 2012. doi:10.1007/s00607-011-0178-4.

[ARHL08] E. Auer, A. Rauh, E.P. Hofer, and W. Luther. Validated modeling of mechanical systems with SMARTMOBILE: Improvement of performance by VALENCIA-IVP. In *Proc. of Dagstuhl Seminar 06021: Reliable Implementation of Real Number Algorithms: Theory and Practice*, Lecture Notes in Computer Science, pages 1-27, 2008.

[ARK20] E. Auer, A. Rauh, and J. Kersten. Experiments-based parameter identification on the GPU for cooperative systems. *Journal of Computational and Applied Mathematics*, 371:112657, 2020. doi:10.1016/j.cam.2019.112657.

[AS03] D. Angeli and E.D. Sontag. Monotone control systems. *IEEE Transactions on Automatic Control*, 48(10):1684-1698, 2003. doi:10.1109/TAC.2003.817920.

[ATDC14] J. Alexandre dit Sandretto, G. Trombettoni, D. Daney, and G. Chabert. Certified calibration of a cable-driven robot using interval contractor programming. In F. Thomas, editor, *Computational Kinematics*, volume 15 of *Mechanisms and Machine Science Series*, pages 209-217. Springer, The Netherlands, Dordrecht, 2014. doi:10.1007/978-94-007-7214-4_24.

[BA21] Y. Becis-Aubry. Ellipsoidal constrained state estimation in presence of bounded disturbances. In *Proc. of the 2021 European Control Conference (ECC)*, pages 555-560, Delft, The Netherlands, 2021. doi:10.23919/ECC54610.2021.9655062.

[BEFB94] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan. *Linear Matrix Inequalities in System and Control Theory*. SIAM, Philadelphia, 1994. doi:10.1137/1.9781611970777.

[Bü18] F. Bünger. Shrink wrapping for Taylor models revisited. *Numerical Algorithms*, 78:1001-1017, 08 2018. doi:10.1007/s11075-017-0410-1.

[Bü20] F. Bünger. A Taylor model toolbox for solving ODEs implemented in MATLAB/INTLAB. *Journal of Computational and Applied Mathematics*, 368:112511, 2020. doi:10.1016/j.cam.2019.112511.

[Bü21] F. Bünger. Preconditioning of Taylor models, implementation and test cases. *Nonlinear Theory and Its Applications, IEICE*, 12(1):2-40, 2021. doi:10.1587/nolta.12.2.

[CA13] E.F. Camacho and C.B. Alba. *Model Predictive Control*. Advanced Textbooks in Control and Signal Processing. Springer London, 2013. doi:10.1007/978-0-85729-398-5.

[CBGV12] P. Collins, D. Bresolin, L. Geretti, and T. Villa. Computing the evolution of hybrid systems using rigorous function calculus. In *Proc. of the 4th IFAC Conference on Analysis and*
Design of Hybrid Systems (ADHS12), pages 284–290, Eindhoven, The Netherlands, 2012.
doi:10.3182/20120606-3-NL-3011.00063.

[CL84] H. Cornelius and R. Lohner. Computing the range of values of real functions with accuracy higher than second order. *Computing*, 33:331–347, 1984. doi:10.1007/BF02242276.

[CR97] T. Csendes and D. Ratz. Subdivision direction selection in interval methods for global optimization. *SIAM Journal on Numerical Analysis*, 34(3):922–938, 1997. doi:10.1137/S0036142995281528.

[DARA13] T. Dötschel, E. Auer, A. Rauh, and H. Aschemann. Thermal behavior of high-temperature fuel cells: Reliable parameter identification and interval-based sliding mode control. *Soft Computing*, 17(8):1329–1343, 2013. doi:10.1007/s00500-013-1003-0.

[dERP22] A. dos Reis de Souza, D. Efimov, T. Raïssi, and X. Ping. Robust output feedback model predictive control for constrained linear systems via interval observers. *Automatica*, 135:109951, 2022. doi:10.1016/j.automatica.2021.109951.

[dRdSER22a] A. dos Reis de Souza, D. Efimov, and T. Raïssi. Robust output feedback model predictive control of time-delayed systems using interval observers. *International Journal of Robust and Nonlinear Control*, 32(3):1180–1193, 2022. doi:10.1002/rnc.5874.

[dRdSER22b] A. dos Reis de Souza, D. Efimov, and T. Raïssi. Robust output feedback MPC for LPV systems using interval observers. *IEEE Transactions on Automatic Control*, 67(6):3188–3195, 2022. doi:10.1109/TAC.2021.3099449.

[Eke10] İ. Eker. Second-order sliding mode control with experimental application. *ISA Transactions*, 49(3):394–405, 2010. doi:10.1016/j.isatra.2010.03.010.

[ERCZ13] D. Efimov, T. Raïssi, S. Chebotarev, and A. Zolghadri. Interval state observer for nonlinear time varying systems. *Automatica*, 49(1):200–205, 2013. doi:10.1016/j.automatica.2012.07.004.

[EVU09] O. Erdinç, B. Vural, and M. Uzunoğlu. A dynamic Lithium-ion battery model considering the effects of temperature and capacity fading. In *Proc. of International Conference on Clean Electrical Power (ICCEP)*, pages 383–386. IEEE, 2009. doi:10.1109/ICCEP.2009.5212025.

[FAdS21] M. Fnadi and J. Alexandre dit Sandretto. Experimental validation of a guaranteed nonlinear model predictive control. *Algorithms*, 14(8):248, 2021. doi:10.3390/a14080248.

[FL02] L. Fridman and A. Levant. Higher order sliding modes. In J.P. Barbot and W. Perruquetti, editors, *Sliding Mode Control in Engineering*, pages 53–101. Marcel Dekker, New York, 2002.

[FLMR95] M. Fliess, J. Lévine, P. Martin, and P. Rouchon. Flatness and defect of nonlinear systems: Introductory theory and examples. *International Journal of Control*, 61:1327–1361, 1995. doi:10.1080/00207179508921959.

[FR24] M. Fnadi and A. Rauh. Exponential state enclosure techniques for the implementation of validated model predictive control. *Acta Cybernetica*, 26(4):839–854, 2024. doi:10.14232/actacyb.302013.

[GT08] M. Gennat and B. Tibken. Computing guaranteed bounds for uncertain cooperative and monotone nonlinear systems. In *Proc. of the 17th IFAC World Congress*, pages 4846–4851. IFAC Proceedings Volumes, 2008. doi:10.3182/20080706-5-KR-1001.00814.

[Hei18] O. Heimlich. Octave forge interval, 2018. https://wiki.octave.org/wiki/index.php?title=Interval_package&oldid=12214, version: 3.2.0, accessed: Sept. 01, 2023.

[Hir19] M. Hirsch. On the nonchaotic nature of monotone dynamical systems. *European Journal of Pure and Applied Mathematics*, 12:680–688, 07 2019. doi:10.29020/nybg.ejpam.v12i3.3492.

[HKRA20] E. Hildebrandt, J. Kersten, A. Rauh, and H. Aschemann. Robust interval observer design for fractional-order models with applications to state estimation of batteries. In *Proc. of the 21st IFAC World Congress*, pages 3683–3688. IFAC-PapersOnLine, 2020. doi:10.1016/j.ifacol.2020.12.2052.

[Hoe01] J. Hoefkens. *Rigorous Numerical Analysis with High-Order Taylor Models*. PhD thesis, Michigan State University, 2001. https://www.bmtdynamics.org/pub/papers/hoefkensphd/hoefkensphd.pdf (accessed: Oct. 16, 2025).

[IPI+22] S. Ifqir, V. Puig, D. Ichalal, N. Ait-Oufroukh, and S. Mammar. Zonotopic set-membership state estimation for switched systems. *Journal of the Franklin Institute*, 359(16):9241–9270, 2022. doi:10.1016/j.jfranklin.2022.08.044.

[Jau09] L. Jaulin. Robust set-membership state estimation; application to underwater robotics. *Automatica*, 45(1):202–206, 2009. doi:10.1016/j.automatica.2008.06.013.
[JKDW01] L. Jaulin, M. Kieffer, O. Didrit, and É. Walter. *Applied Interval Analysis*. Springer-Verlag, London, 2001. doi:10.1007/978-1-4471-0249-6.

[JUDW00] S.J. Julier, J.K. Uhlmann, and H.F. Durrant-Whyte. A new approach for the nonlinear transformation of means and covariances in filters and estimators. *IEEE Transactions on Automatic Control*, 45(3):477-482, 2000. doi:10.1109/9.847726.

[Kac02] T. Kaczorek. *Positive 1D and 2D Systems*. Springer-Verlag, London, 2002. doi:10.1007/978-1-4471-0221-2.

[KAH07] M. Kletting, F. Antritter, and E.P. Hofer. Robust flatness based controller design using interval methods. *IFAC Proceedings Volumes*, 40(12):876-881, 2007. 7th IFAC Symposium on Nonlinear Control Systems. doi:10.3182/20070822-3-ZA-2920.00145.

[Kal60] R.E. Kalman. A new approach to linear filtering and prediction problems. *Transaction of the ASME - Journal of Basic Engineering*, 82(Series D):35-45, 1960. doi:10.1115/1.3662552.

[KD09] A. Der Kiureghian and O. Ditlevsen. Aleatory or epistemic? Does it matter? *Structural Safety*, 31(2):105-112, 2009. Risk Acceptance and Risk Communication. doi:10.1016/j.strusafe.2008.06.020.

[Kea23] R.B. Kearfott. Interval computing. In Tsau-Young Lin, Churn-Jung Liau, and Janusz Kacprzyk, editors, *Granular, Fuzzy, and Soft Computing*, pages 701-709. Springer US, New York, NY, 2023. doi:10.1007/978-1-0716-2628-3_722.

[Kei93] J. B. Keiper. Interval arithmetic in mathematica. *Interval Computations*, 3:76-87, 1993.

[KMWZ20] T. Kapela, M. Mrozek, D. Wilczak, and P. Zgliczynski. CAPD::DynSys: A flexible C++ toolbox for rigorous numerical analysis of dynamical systems. *Communications in Nonlinear Science and Numerical Simulation*, page 105578, 11 2020. doi:10.1016/j.cnsns.2020.105578.

[Krä12] W. Krämer. XSC Languages (C-XSC, PASCAL-XSC) — Scientific computing with validation, arithmetic requirements, hardware solution and language support, 2012. www.math.uni-wuppertal.de/~xsc/, C-XSC 2.5.3, acessed: Sept. 01, 2023.

[Kra69] R. Krawczyk. Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehlerschranken. *Computing*, 4:189-201, 1969. In German. doi:10.1007/BF02234767.

[KRA18] J. Kersten, A. Rauh, and H. Aschemann. State-space transformations of uncertain systems with purely real and conjugate-complex eigenvalues into a cooperative form. In *Proc. of 23rd Intl. Conf. on Methods and Models in Automation and Robotics 2018*, Miedzyzdroje, Poland, 2018. doi:10.1109/MMAR.2018.8486085.

[KRA19] J. Kersten, A. Rauh, and H. Aschemann. Transformation of uncertain linear fractional order differential equations into a cooperative form. In *Proc. of the 11th IFAC Symposium on Nonlinear Control Systems (NOLCOS)*, pages 646-651. IFAC-PapersOnLine, 2019. doi: 10.1016/j.ifacol.2019.12.035.

[KRA20] J. Kersten, A. Rauh, and H. Aschemann. Interval methods for robust gain scheduling controllers: An LMI-based approach. *Granular Computing*, pages 203-216, 2020. doi: 10.1007/s41066-018-00147-1.

[KRHA06] M. Kletting, A. Rauh, E.P. Hofer, and H. Aschemann. Interval observer design based on Taylor models for nonlinear uncertain continuous-time systems. In *12th GAMM-IMACS International Symposium on Scientific Computing, Computer Arithmetic and Validated Numerics (SCAN 2006)*, pages 25-25, 2006. doi:10.1109/SCAN.2006.26.

[KRK+10] I. Krasnochtanova, A. Rauh, M. Kletting, H. Aschemann, E.P. Hofer, and K.-M. Schoop. Interval methods as a simulation tool for the dynamics of biological wastewater treatment processes with parameter uncertainties. *Applied Mathematical Modeling*, 34(3):744-762, 2010. doi:10.1016/j.apm.2009.06.019.

[KSY22] M. Khajenejad, F. Shoaib, and S.Z. Yong. Interval observer synthesis for locally Lipschitz nonlinear dynamical systems via mixed-monotone decompositions. In *2022 American Control Conference (ACC)*, pages 2970-2975, 2022. doi:10.23919/ACC53348.2022.9867741.

[Küh99] W. Kühn. Rigorous error bounds for the initial value problem based on defect estimation. Technical report, 1999. http://www.decatur.de/personal/papers/defect.zip (accessed: July 14, 2025).

[KV97] A.B. Kurzhanskii and I. Vályi. *Ellipsoidal Calculus for Estimation and Control*. Birkhäuser, Boston, MA, 1997.
[KY23] M. Khajenejad and S.Z. Yong. Tight remainder-form decomposition functions with applications to constrained reachability and guaranteed state estimation. *IEEE Transactions on Automatic Control*, pages 1–16, 2023. doi:10.1109/TAC.2023.3250515.

[Loh87] R. Lohner. Enclosing the solutions of ordinary initial and boundary value problems. In *E.W. Kaucher, U.W. Kulisch, and C. Ullrich, editors, Computer Arithmetic: Scientific Computation and Programming Languages*, pages 255–286, Stuttgart, 1987. Wiley-Teubner Series in Computer Science.

[LR23] M. Lahme and A. Rauh. Online identification of the open-circuit voltage characteristic of Lithium-ion batteries with a contractor-based procedure. In *IEEE International Conference on Methods and Models in Automation and Robotics MMAR*, Miedzyzdroje, Poland, 2023. doi:10.1109/MMAR58394.2023.10242524.

[LR24] M. Lahme and A. Rauh. Set-valued approach for the online identification of the open-circuit voltage of Lithium-ion batteries. *Acta Cybernetica*, 26(4):855–869, 2024. doi:10.14232/actacyb.301185.

[LRD24] M. Lahme, A. Rauh, and G. Defresne. Interval observer design for an uncertain time-varying quasi-linear system model of Lithium-ion batteries. In *Proceedings of the 2024 European Control Conference (ECC)*, Stockholm, Sweden, 2024. doi:10.23919/ECC64448.2024.10591102.

[LS07] Y. Lin and M.A. Stadtherr. Validated solutions of initial value problems for parametric ODEs. *Applied Numerical Mathematics*, 57(10):1145–1162, 2007. doi:10.1016/j.apnum.2006.10.006.

[Mar03] H.J. Marquez. *Nonlinear Control Systems*. John Wiley & Sons, Inc., New Jersey, 2003.

[May17] G. Mayer. *Interval Analysis and Automatic Result Verification*. De Gruyter Studies in Mathematics. De Gruyter, Berlin/Boston, 2017. doi:10.1515/9783110499469.

[MB10] F. Mazenc and O. Bernard. Asymptotically stable interval observers for planar systems with complex poles. *IEEE Transactions on Automatic Control*, 55(2):523–527, 2010. doi:10.1109/TAC.2009.2037472.

[Moo66] R.E. Moore. *Interval Arithmetic*. Prentice-Hall, Englewood Cliffs, N.J., 1966.

[Mül27] M. Müller. Über die Eindeutigkeit der Integrale eines Systems gewöhnlicher Differenzialgleichungen und die Konvergenz einer Gattung von Verfahren zur Approximation dieser Integrale. In *Sitzungsbericht Heidelberger Akademie der Wissenschaften*, 1927. In German.

[Ned99] N.S. Nedialkov. *Computing Rigorous Bounds on the Solution of an Initial Value Problem for an Ordinary Differential Equation*. PhD thesis, Graduate Department of Computer Science, University of Toronto, 1999.

[Ned07] N.S. Nedialkov. Interval tools for ODEs and DAEs. In *CD-Proc. of 12th GAMM-IMACS Intl. Symposium on Scientific Computing, Computer Arithmetic, and Validated Numerics SCAN 2006*, Duisburg, Germany, 2007. IEEE Computer Society. doi:10.1109/SCAN.2006.28.

[Ned11] N.S. Nedialkov. Implementing a rigorous ODE solver through literate programming. In *A. Rauh and E. Auer, editors, Modeling, Design, and Simulation of Systems with Uncertainties*, Mathematical Engineering, pages 3–19. Springer, Berlin, Heidenberg, 2011. doi:10.1007/978-3-642-15956-5_1.

[Neu90] A. Neumaier. *Interval Methods for Systems of Equations*. Cambridge University Press, Encyclopedia of Mathematics, Cambridge, 1990. doi:10.1017/CBO9780511526473.

[RA12] A. Rauh and H. Aschemann. Interval-based sliding mode control and state estimation for uncertain systems. In *2012 17th International Conference on Methods and Models in Automation and Robotics (MMAR)*, pages 595–600, 2012. doi:10.1109/MMAR.2012.6347818.

[RAA10] A. Rauh, E. Auer, and H. Aschemann. Development of a quality measure for the characterization of guaranteed solution sets to ODEs in engineering. In *8th IFAC Symposium on Nonlinear Control Systems*, volume 43, pages 48–53, Bologna, Italy, 2010. doi:10.3182/20100901-3-IT-2016.00162.

[RBJ21] A. Rauh, A. Bourgois, and L. Jaulin. Union and intersection operators for thick ellipsoid state enclosures: Application to bounded-error discrete-time state observer design. *Algorithms*, 14(3):88, 2021. doi:10.3390/a14030088.

[RBJK21] A. Rauh, A. Bourgois, L. Jaulin, and J. Kersten. Ellipsoidal enclosure techniques for a verified simulation of initial value problems for ordinary differential equations. In *Proc. of 5th International Conference on Control, Automation and Diagnosis (ICCAD'21)*, Grenoble, France, 2021. doi:10.1109/ICCAD52417.2021.9638755.
[RCD<sup>+</sup>22a] A. Rauh, T. Chevet, T.N. Dinh, J. Marzat, and T. Raïssi. Robust iterative learning observers based on a combination of stochastic estimation schemes and ellipsoidal calculus. In *25th International Conference on Information Fusion (FUSION)*, pages 1–8, Linköping, Sweden, 2022. doi:10.23919/FUSION49751.2022.9841329.

[RCD<sup>+</sup>22b] A. Rauh, T. Chevet, T.N. Dinh, J. Marzat, and T. Raïssi. A stochastic design approach for iterative learning observers for the estimation of periodically recurring trajectories and disturbances. In *Proc. of the 20th European Control Conference ECC*, London, UK, 2022. doi:10.23919/ECC55457.2022.9838288.

[RD24] S. Rohou and B. Desrochers. The Codac library: A catalog of domains and contractors. *Acta Cybernetica*, 26(4):871–887, 2024. doi:10.14232/actacyb.302772.

[RDA11] A. Rauh, T. Dötschel, and H. Aschemann. Experimental parameter identification for a control-oriented model of the thermal behavior of high-temperature fuel cells. In *CD-Proc. of IEEE International Conference on Methods and Models in Automation and Robotics MMAR*, Miedzyzdroje, Poland, 2011. doi:10.1109/MMAR.2011.6031306.

[RE18] T. Raïssi and D. Efimov. Some recent results on the design and implementation of interval observers for uncertain systems. *at-Automatisierungstechnik*, 66(3):213–224, 2018. doi:10.1515/auto-2017-0081.

[REZ12] T. Raïssi, D. Efimov, and A. Zolghadri. Interval state estimation for a class of nonlinear systems. *IEEE Trans. Automat. Contr.*, 57:260–265, 2012. doi:10.1109/TAC.2011.2164820.

[RH09] A. Rauh and E.P. Hofer. Interval methods for optimal control. In A. Frediani and G. Buttazzo, editors, *Proceedings of the 47th Workshop on Variational Analysis and Aerospace Engineering 2007*, pages 397–418, School of Mathematics, Erice, Italy, 2009. Springer–Verlag. doi:10.1007/978-0-387-95857-6_22.

[RJ21a] A. Rauh and L. Jaulin. A computationally inexpensive algorithm for determining outer and inner enclosures of nonlinear mappings of ellipsoidal domains. *International Journal of Applied Mathematics and Computer Science AMCS*, 31(3):399–415, 2021. doi:10.34768/amcs-2021-0027.

[RJ21b] S. Rohou and L. Jaulin. Exact bounded-error continuous-time linear state estimator. *Systems & Control Letters*, 153:104951, 2021. doi:10.1016/j.sysconle.2021.104951.

[RKA11] A. Rauh, I. Krasnochtanova, and H. Aschemann. Quantification of overestimation in interval simulations of uncertain systems. In *Proceedings of IEEE International Conference on Methods and Models in Automation and Robotics MMAR 2011*, Miedzyzdroje, Poland, 2011. doi:10.1109/MMAR.2011.6031328.

[RKA17] A. Rauh, J. Kersten, and H. Aschemann. Interval-based implementation of robust variable-structure and backstepping controllers of single-input single-output systems. *IFAC-PapersOnLine*, 50(1):6283–6288, 2017. 20th IFAC World Congress. doi:10.1016/j.ifacol.2017.08.865.

[RKA20a] A. Rauh, J. Kersten, and H. Aschemann. Interval and linear matrix inequality techniques for reliable control of linear continuous-time cooperative systems with applications to heat transfer. *International Journal of Control*, 93(11):2771–2788, 2020. doi:10.1080/00207179.2019.1708966.

[RKA20b] A. Rauh, J. Kersten, and H. Aschemann. Interval methods and contractor-based branch-and-bound procedures for verified parameter identification of quasi-linear cooperative system models. *Journal of Computational and Applied Mathematics*, 367:112484, 2020. doi:10.1016/j.cam.2019.112484.

[RLRR22] B.S. Rego, D. Locatelli, D.M. Raimondo, and G.V. Raffo. Joint state and parameter estimation based on constrained zonotopes. *Automatica*, 142:110425, 2022. doi:10.1016/j.automatica.2022.110425.

[RMAR16] J. Reuter, E. Mank, H. Aschemann, and A. Rauh. Battery state observation and condition monitoring using online minimization. In *Proc. of 21st Intl. Conference on Methods and Models in Automation and Robotics*, Miedzyzdroje, Poland, 2016. doi:10.1109/MMAR.2016.7575313.

[Roh19] J. Rohn. VERSOFT, 2019. http://www.cs.cas.cz/~rohn/matlab/versoft_10.zip, version: 10, accessed: Sept. 01, 2023.
[RRJ22] A. Rauh, S. Rohou, and L. Jaulin. An ellipsoidal predictor-corrector state estimation scheme for linear continuous-time systems with bounded parameters and bounded measurement errors. *Frontiers in Control Engineering*, 3:785795, 2022. doi:10.3389/fcteg.2022.785795.

[RS16] A. Rauh and L. Senkel. Interval methods for robust sliding mode control synthesis of high-temperature fuel cells with state and input constraints. In A. Rauh and L. Senkel, editors, *Variable-Structure Approaches for Analysis, Simulation, Robust Control and Estimation of Uncertain Dynamic Systems*, Mathematical Engineering, pages 53–85. Springer–Verlag, 2016. doi:10.1007/978-3-319-31539-3_3.

[RSA15] A. Rauh, L. Senkel, and H. Aschemann. Interval-based sliding mode control design for solid oxide fuel cells with state and actuator constraints. *IEEE Transactions on Industrial Electronics*, 62(8):5208–5217, 2015. doi:10.1109/TIE.2015.2404811.

[RSA16] A. Rauh, L. Senkel, and H. Aschemann. Interval methods for variable-structure control of dynamic systems with state constraints. In *2016 3rd Conference on Control and Fault-Tolerant Systems (SysTol)*, pages 466–471, 2016. doi:10.1109/SYSTOL.2016.7739793.

[RSKA16] A. Rauh, L. Senkel, J. Kersten, and H. Aschemann. Reliable control of high-temperature fuel cell systems using interval-based sliding mode techniques. *IMA Journal of Mathematical Control and Information*, 33(2):457–484, 2016. doi:10.1093/imamci/dnu051.

[Rum99] S.M. Rump. INTLAB — INTerval LABoratory. In T. Csendes, editor, *Developments in Reliable Computing*, pages 77–104. Kluver Academic Publishers, 1999. doi:10.1007/978-94-017-1247-7_7.

[RWH$^{+}$22] A. Rauh, S. Wirtensohn, P. Hoher, J. Reuter, and L. Jaulin. Reliability assessment of an Unscented Kalman Filter by using ellipsoidal enclosure techniques. *Mathematics*, 10(16):3011, 2022. doi:10.3390/math10163011.

[SB14] D.P. Sanders and L. Benet. IntervalArithmetic.jl, 2014. https://github.com/JuliaIntervals/IntervalArithmetic.jl, DOI: 10.5281/zenodo.3336308, accessed: Sept. 01, 2023.

[Sch79] J. Schnellinger. Das gekürzte Rechnen. *Zeitschrift für das Realschulwesen*, pages 257–284, 1879. In German.

[Smi95] H. L. Smith. *Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems*, volume 41. Mathematical Surveys and Monographs, American Mathematical Soc., Providence, 1995.

[SRA16] L. Senkel, A. Rauh, and H. Aschemann. Experimental and numerical validation of a reliable sliding mode control strategy considering uncertainty with interval arithmetic. In A. Rauh and L. Senkel, editors, *Variable-Structure Approaches for Analysis, Simulation, Robust Control and Estimation of Uncertain Dynamic Systems*, Mathematical Engineering, pages 87–122. Springer–Verlag, 2016. doi:10.1007/978-3-319-31539-3_4.

[Ste94] R. Stengel. *Optimal Control and Estimation*. Dover Publications, Inc., 1994.

[TGT09] K.P. Tee, S. S. Ge, and E. H. Tay. Barrier Lyapunov functions for the control of output-constrained nonlinear systems. *Automatica*, 45(4):918–927, 2009. doi:10.1016/j.automatica.2008.11.017.

[Utk92] V.I. Utkin. *Sliding Modes in Control and Optimization*. Communications and Control Engineering. Springer, Berlin, Heidelberg, 1992. doi:10.1007/978-3-642-84379-2.

[WJ94] É. Walter and L. Jaulin. Guaranteed characterization of stability domains via set inversion. *IEEE Transactions on Automatic Control*, 39(4):886–889, 1994. doi:10.1109/9.286277.

[WLS18] Z. Wang, C.-C. Lim, and Y. Shen. Interval observer design for uncertain discrete-time linear systems. *Systems & Control Letters*, 116:41–46, 2018. doi:10.1016/j.sysconle.2018.04.003.

[WSM15] B. Wang, W. Shi, and Z. Miao. Confidence analysis of standard deviational ellipse and its extension into higher dimensional euclidean space. *PLOS ONE*, 10(3):1–17, 03 2015. doi:10.1371/journal.pone.0118537.