# On the Computation of Sensitivity Tubes

Andrea Pupa, Tommaso Belvedere, Cristian Secchi, Paolo Robuffo Giordano

► To cite this version:

Andrea Pupa, Tommaso Belvedere, Cristian Secchi, Paolo Robuffo Giordano. On the Computation of Sensitivity Tubes. IEEE Robotics and Automation Letters, 2025, 10 (9), pp.8802 - 8809. 10.1109/LRA.2025.3587562. hal-05158394

HAL Id: hal-05158394

https://hal.science/hal-05158394v1

Submitted on 11 Jul 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>
# On the Computation of Sensitivity Tubes

Andrea Pupa<fnref n="1" />, Tommaso Belvedere<fnref n="2" />, Cristian Secchi<fnref n="1" /> and Paolo Robuffo Giordano<fnref n="2" />

**Abstract**—Achieving robust robot control requires explicit treatment of model uncertainties. Closed-loop sensitivity has emerged as a powerful tool to analyze how parameter errors map into state and input deviations through so-called “sensitivity tubes”, traditionally built from ellipsoidal uncertainty sets and used to robustify system constraints. These ellipsoids, however, are themselves smooth approximations of underlying hyperboxes in the parameter space, leading to an inaccurate estimation of the parameter set. This paper extends that framework by proposing two new formulations that more precisely represent the real closed-loop behavior of the system through improved computation of the sensitivity tubes. The first constructs tubes directly from hyperboxes, exactly preserving the original parameter bounds but producing a non-differentiable description. The second employs superquadrics, which smoothly approximate the hyperbox with user-tunable fidelity while preserving differentiability, as in the ellipsoidal case. Both methods are validated through an extensive simulation campaign, where the resulting input tubes ensure actuator constraints are respected. The results demonstrate that the new tubes better enclose the perturbed trajectories with respect to ellipsoidal ones, enhancing robustness for both online and offline trajectory planning.

**Index Terms**—Planning under Uncertainty, Optimization and Optimal Control, Integrated Planning and Control

## I. INTRODUCTION

ROBOTS are now more widespread than ever across various fields. In industry, humans and robots collaborate in close proximity without the need for physical barriers. In operating rooms, robotic systems enhance surgical precision, particularly in procedures such as laparoscopic surgery. In agriculture, swarms of drones are deployed to monitor vast fields of crops efficiently. As a result, it is crucial to ensure that robots can perform the required tasks safely and reliably despite the unavoidable uncertainties affecting their models.

Various paradigms have been explored for controlling robots affected by model uncertainty. One approach involves offline identification of model parameters, but this often requires collecting large amounts of data, which can itself be unsafe when uncertainties are substantial [1]. Alternatively, online

adaptive methods [2] can compensate for uncertainties or disturbances; however, they may degrade transient performance, compromise system stability, and prove unsuitable for tasks where parameters change abruptly. Robust control techniques [3], such as sliding mode control [4] or passivity-based methods [5], [6], offer another solution. However, these methods rely on specific model properties, often require high actuation bandwidth, and typically result in performance degradation even when parameters are well estimated (e.g., they are tuned for the worst case).

A complementary approach, which can be integrated with the aforementioned methods, involves planning and executing motions that result in inherently reduced sensitivity to uncertainties. In this context, closed-loop sensitivity has recently emerged as a valuable tool for handling uncertainties in system parameters [7]. This technique focuses on analyzing how variations in model parameters affect the behavior of a system when operating under feedback control, and can be used to generate intrinsically robust motion plans. The main advantage of this approach is that it does not require the design or use of a specific controller, allowing it to be successfully applied to many different cases. For instance, it has been effectively employed and validated experimentally in quadrotor trajectory tracking [8], energy tank initialization for manipulators [9], and sampling-based motion planning [10]. A significant development in this area is the introduction of sensitivity tubes, as proposed in [11]. These tubes represent an envelope that bounds the perturbed trajectories of the system states and inputs under the effect of parametric uncertainties. This has been demonstrated to be a powerful, direct approach to account for such uncertainties in safety critical offline planning [11] and robust Model Predictive Control [12]. To construct these tubes, one assumes that the model parameters lie within a bounded range centered around a nominal value—typically modeled as a hyperbox. The corresponding sensitivity matrices are then used to estimate how deviations from the nominal parameters influence the state and input trajectories. This approach effectively captures the worst-case impact of uncertainty and allows for robust constraint satisfaction in trajectory optimization. In [11], however, the parameter set is approximated by the ellipsoid inscribed in the hyperbox for obtaining a tractable and differentiable formulation of the tubes suitable for gradient-based trajectory optimization. Unfortunately, this ellipsoidal approximation can be quite restrictive since, as the number of parameters increases, the ellipsoid-to-hyperbox volume ratio tends to zero, excluding many valid parameter combinations near the corners. This under-approximates system perturbations and degrades predictions of uncertainty effects.

To address this issue, this work extends [11] by introducing

Manuscript received: March, 4, 2025; Revised May, 28, 2025; Accepted June, 27, 2025.

This paper was recommended for publication by Editor L. Pallottino upon evaluation of the Associate Editor and Reviewers' comments. This work was supported by the project ANR-20-CE33-0003 "CAMP" and by the National Recovery and Resilience Plan (NRRP), Mission 04 Component 2 Investment 1.5 - NextGenerationEU, Call for tender n. 3277 dated 30/12/2021. Award Number: 0001052 dated 23/06/2022.

Digital Object Identifier (DOI): see top of this page.

<footnote n="1">¹ A. Pupa and C. Secchi are with the Department of Science and Methods of Engineering, University of Modena and Reggio Emilia, Italy. E-mail: {andrea.pupa, cristian.secchi}@unimore.it.</footnote>

<footnote n="2">²T. Belvedere and P. Robuffo Giordano are with CNRS, Univ Rennes, Inria, IRISA, Campus de Beaulieu, 35042 Rennes Cedex, France. E-mail: {tommaso.belvedere, prg}@irisa.fr.</footnote>
three methods for computing sensitivity tubes. The first computes tubes directly assuming parameter deviations lie within a hyperbox, avoiding set approximation but resulting in non-differentiable tubes. The second method is the one described in [11], where the hyperbox is approximated by an ellipsoid, but with an extended formal proof of the mapping from the parameter space to the state space. Lastly, the third method uses superquadrics to approximate the hyperbox with arbitrary precision [13], [14] while maintaining differentiability, offering a compromise between the previous two methods.

Summarizing, the main contributions of this work are:

* Two new methods to compute the sensitivity tubes: one based on hyperboxes, suitable for gradient-free optimization algorithms (such as MPPI), and one based on superquadrics, suitable for gradient-based optimization.

* The formal proofs of these two new methods as well as the one previously proposed in [11];

* A comprehensive validation of the new methods in a simulated environment, including also a comparison with [11], demonstrating better performance in encapsulating the real behavior of the system.

The paper is organized as follows: Sec. II provides the theoretical background on closed-loop sensitivity. The first part of Sec. III formally defines the sensitivity tube radius and presents the three methods for computing sensitivity tubes; while Sec. III-A–III-C detail the corresponding proofs for each method and may be skipped on an initial reading. Lastly, Sec. IV presents the experimental validation across different simulation environments, and Sec. V discusses conclusions and future work.

## II. CLOSED-LOOP SENSITIVITY

For the reader's convenience, we now recall the main notion related to the closed-loop sensitivity and associated quantities. Consider a generic robot modeled by the following ordinary differential equation:

$$ \dot{x} = f(x, u, p), \quad (1) $$

where $x \in \mathbb{R}^{n_x}$ represents the state and $u \in \mathbb{R}^{n_u}$ is the control input. Vector $p \in \mathbb{R}^{n_p}$, instead, represents a set of model or environmental parameters whose true values are not perfectly known. For example, in a robot arm, such parameters may represent the position of the center of mass of a link, the friction coefficient of the motors, or the elastic properties of the material it is interacting with.

In the context of trajectory tracking, consider an output function $y(x) \in \mathbb{R}^{n_y}$, e.g., the manipulator pose, that needs to follow a desired reference trajectory $y_d(t) \in \mathbb{R}^{n_y}$, and an associated controller with generic expression

$$ \left\{ \begin{array}{l} \dot{\xi} = g(\xi, x, y_d(t), p_n, k_c, t), \\ u = h(\xi, x, y_d(t), p_n, k_c, t), \end{array} \right. \quad (2) $$

where $\xi \in \mathbb{R}^{n_\xi}$ are the internal states of the controller, $p_n \in \mathbb{R}^{n_p}$ are a nominal value for the model parameters, i.e., the estimated nominal value of $p$ used by the controller that may differ from the real $p$ due to the uncertainties, and $k_c \in \mathbb{R}^{n_k}$ are the controller gains. When $p \neq p_n$, the so-called perturbed case, the controller cannot perfectly realize

the tracking task, e.g. the error may only converge to zero after some transient, or converge to a bounded set around the desired trajectory, and safety constraints might be violated even if the nominal trajectory is feasible. It is therefore crucial to develop a framework that, being aware of model uncertainties, can not only minimize their effects on the task execution but also (and perhaps most importantly) guarantee constraint satisfaction in the real (perturbed) case.

In order to address this problem, it is necessary to characterize how perturbations in the parameters $p$ affect the system state $x$, the input $u$, and the output $y$ trajectories in closed-loop during motion. To this aim, it is possible to exploit the closed-loop sensitivity notion introduced in [7] to describe how much the (closed-loop) system is locally sensitive to parametric uncertainties. Given the closed-loop system (1)–(2), one defines the so-called state sensitivity matrix:

$$ \Pi(t) = \left. \frac{\partial x(t)}{\partial p} \right|_{p=p_n} \in \mathbb{R}^{n_x \times n_p}, \quad (3) $$

and input sensitivity matrix:

$$ \Theta(t) = \left. \frac{\partial u(t)}{\partial p} \right|_{p=p_n} \in \mathbb{R}^{n_u \times n_p}, \quad (4) $$

which quantify how variations in the parameters $p$ (w.r.t. the nominal $p_n$) affect the states and inputs during motion. In the general case, a closed-form expression for $\Pi(t)$ and $\Theta(t)$ is not available. However, these quantities can be easily obtained by forward integrating the following system of differential equations:

$$ \left\{ \begin{array}{l} \dot{\Pi}(t) = \frac{\partial f}{\partial x}\Pi + \frac{\partial f}{\partial u}\Theta + \frac{\partial f}{\partial p}, \quad \Pi(t_0) = 0, \\ \dot{\Pi}_{\xi}(t) = \frac{\partial g}{\partial x}\Pi + \frac{\partial g}{\partial \xi}\Pi_{\xi}, \quad \Pi_{\xi}(t_0) = 0, \\ \Theta(t) = \frac{\partial h}{\partial x}\Pi + \frac{\partial h}{\partial \xi}\Pi_{\xi}, \end{array} \right. \quad (5) $$

where $\Pi_{\xi} \in \mathbb{R}^{n_\xi \times n_p}$ is the sensitivity of the controller states.

## III. SENSITIVITY TUBES COMPUTATION

To compute sensitivity tubes, we assume that the actual parameters lie within a compact and known set $\mathcal{P}$. The aim is to assess how variations in these parameters w.r.t. a nominal value affect the closed-loop trajectory. Such an analysis enables the construction of the so-called sensitivity tubes, a set of envelopes that bundles all possible perturbed trajectories resulting from parameter deviations, which is essential information for robust trajectory planning.

Rather than describing state/input tubes based on the geometric shape of $\mathcal{P}^1$, the focus of the proposed analysis is in obtaining a general measure: the maximum deviation $\alpha$ that the perturbed system can exhibit along a given direction $n \in \mathbb{R}^{n_x}$. This deviation $\alpha$ represents the tube radius along that direction. It is worth underlining that, although the computed value of $\alpha$ will depend on the geometric shape of $\mathcal{P}$, e.g. box versus ellipsoids, the definition of $\alpha$ as the maximum deviation along $n$ is independent from the geometric shape $\mathcal{P}$.

<footnote n="1">¹Which is, for example, the case in [15], where the authors characterize tubes using ellipsoids, an approach that cannot be directly generalized to other geometric shapes.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Representation of the mapping obtained with the sensitivity matrix $\Pi$. On the left the hyperbox $\mathcal{B}_p$, the ellipsoid $\mathcal{E}_p$, and the superquadric $\mathcal{Q}_p$ in the parameter space. On the right, the corresponding shapes $\mathcal{B}_{\mathcal{X}}$, $\mathcal{E}_{\mathcal{X}}$, and $\mathcal{Q}_{\mathcal{X}}$ in the state space.</figcaption></figure>

To carry out the following analysis, the set $\mathcal{P}$ is assumed to be a symmetric hyperbox<fnref n="2" />. That is, each $p_i$ can deviate from its nominal value $p_{n,i}$ by at most $\Delta p_{\max,i} > 0$, such that $p_i \in [p_{n,i} - \Delta p_{\max,i}, p_{n,i} + \Delta p_{\max,i}]$. This is a reasonable assumption as, in many applications, physical parameters can be associated with an uncertainty range due to measurement errors or task specifications, e.g., a minimum/maximum payload weight. Consequently, the set of the parameters can be defined as:

$$ \mathcal{P} := \mathcal{B}_p = \{\mathbf{p} \in \mathbb{R}^{n_p} : -\Delta \mathbf{p}_{\max} \le \mathbf{p} - \mathbf{p}_n \le \Delta \mathbf{p}_{\max}\}, \quad (6) $$

with inequalities intended element-wise. In the following, three different methods to compute $\alpha$ are proposed and analyzed. Such methods exploit the closed-loop state sensitivity matrix $\Pi$ to map finite perturbations from parameter space to state space. This mapping is graphically represented in Fig. 1. We note that, although the following derivations focus on the state space, the proposed method can be seamlessly applied to any quantity whose sensitivity matrix is available, e.g., the input space by using $\Theta$, or any submanifolds arising from lower-dimensional functions of the state and/or inputs.

Let $\Delta x(t) = x(t) - x_n(t)$: for small enough $\Delta p$, $\Delta x(t)$ can be linearly approximated as

$$ \Delta x(t) \approx \Pi(t) \Delta p. \quad (7) $$

Starting from this linearization, the tube radius along an arbitrary direction can be computed for all the geometric shapes considered, namely the hyperbox, the ellipsoids, and the superquadric as:

$$ \alpha_{\mathcal{B}} = \|n^{\top} \Pi W^{\frac{1}{2}}\|_{1} \quad (8) $$

$$ \alpha_{\mathcal{E}} = \sqrt{n^{\top} \Pi W \Pi^{\top} n} \quad (9) $$

$$ \alpha_{\mathcal{Q}} = \sqrt[\frac{2k}{2k-1}]{n^{\top} \Pi W^{\frac{1}{2}} (W^{\frac{1}{2}} \Pi^{\top} n)^{\circ \frac{1}{2k-1}}} \quad (10) $$

where $k \ge 1$, superscript $\circ$ denotes the element-wise power [16], e.g. $\Delta p^{\circ j \top} = [\Delta p_1^j, \dots, \Delta p_{n_p}^j]$, and $W = \text{diag}\{\Delta p_{\max,i}^2\}$. The detailed proofs of (8)–(10) are presented in the following Sec. III-A–Sec. III-C.

The first method derives directly from the definition of $\Delta p$.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Graphical representation of the tube radius $\alpha$ along a direction $n$. Regardless of the analyzed geometric shape, $\alpha$ denotes the distance from the center of the shape to the tangent plane whose normal is $n$.</figcaption></figure>

Since equation (6) represents a hyperbox centered in $p_n$, it can be demonstrated (see Section III-A), that the maximum state deviation occurs when the parameter deviations align with one of the hyperbox corners. However, mapping such a corner results in a non-differentiable expression for $\alpha$, which is undesirable for gradient-based trajectory optimization algorithms.

To ensure differentiability, the original hyperbox can be approximated by its inscribed ellipsoid as proposed in [11]. Using matrix $\Pi$, this ellipsoid in parameter space is mapped into a corresponding ellipsoid in state space. The tube radius is then obtained by evaluating the maximum variation of the ellipsoid along a generic direction $n$.

However, the original set of parameter variations defined in (6) is an hyperbox: as the number of parameters increases, the portion of volume of the box occupied by the inscribed ellipsoid decreases which in turn implies an increased probability that the real parameter deviation will lie outside the ellipsoid. This results in tubes (and tube radii) that can underestimate the effects of uncertainties and lead to incorrect predictions in trajectory planning. Of course, one could approximate the hyperbox with the circumscribed ellipsoid. However, this solution would clearly lead to conservative behaviors.

The third method addresses this issue by utilizing superquadrics [13], specifically superellipsoids, a shape that more closely resembles the original hyperbox. As defined in Sect. III-C, the superellipsoid reduces to a standard ellipsoid when of order $k = 1$, while it converges to the hyperbox as $k \to \infty$.

It is worth emphasizing that, regardless of the analyzed geometric shape, the tube radius $\alpha$ always maintains the same geometric meaning. As illustrated in Fig. 2, from a geometric perspective, $\alpha$ represents the distance from the center of the analyzed shape, i.e., $p_n$, to the tangent plane with normal vector $n$.

**Remark 1.** For all three methods, the tube radius expression has been derived for a generic direction $n$. Often, one needs to evaluate how each component of the state deviates from its nominal value. To evaluate this for the $i$-th component, one can set $n = e_i$, where $e_i$ denotes the $i$-th vector of the canonical basis ($e_1, \dots, e_{n_x}$) of $\mathbb{R}^{n_x}$. This enables, for instance, a direct check of whether the deviation of the $i$-th state component remains within its specified bounds in the worst-case scenario, thereby allowing robust constraint enforcement.

The detailed proofs of (8)–(10) for the three considered cases are now provided. Readers not interested in the technical

<footnote n="2">²While symmetry is not strictly required, it simplifies the mathematical formulation.</footnote>
details may skip these sections and proceed directly to the results in Sec. IV.

## A. Hyperbox

Assuming that the uncertain parameters belong to the set (6), the sensitivity tube radius can be computed by finding which vertex of the hyperbox generates the maximum state deviation along a given direction $n$.

**Proposition 1.** The maximum state deviation occurs when the real parameters lie at one of the vertices, depending on the closed-loop sensitivity matrix $\Pi$, and the tube radius $\alpha_B$ along a generic direction $n$ can be computed using equation (8).

*Proof.* The goal is to obtain the maximum deviation $\alpha_B$ such that:

$$ \alpha_B = n^\top \Delta x \wedge -\Delta p_{\max} \le \Delta p \le \Delta p_{\max}, \quad (11) $$

where, for ease of notation, the dependency of time has been dropped. Mathematically this can be expressed by the following optimization problem:

$$ \alpha_B^* = \max_{\Delta p} n^\top \Pi \Delta p \\ \text{s.t.} \quad -\Delta p_{\max} \le \Delta p \le \Delta p_{\max}, \quad (12) $$

where (7) has been used. Analyzing the cost function, it is possible to rewrite it as a sum of scalar elements:

$$ \alpha_B = \sum_{j=1}^{n_p} \sum_{i=1}^{n_x} n_i \Pi_{ij} \Delta p_j = \sum_{j=1}^{n_p} a_j \Delta p_j, \quad (13) $$

where $a_j$ is the j-th element of the vector $a = n^\top \Pi$. By denoting as $a^+$ and $a^-$ the positive and negative elements of $a$, it is possible to split (13) as:

$$ \alpha_B = \sum_{i \in A^+} a_i^+ \Delta p_i + \sum_{j \in A^-} a_j^- \Delta p_j, \quad (14) $$

where $A^+$ and $A^-$ are the set of the indexes of positive and negative elements of $a$, respectively.

It is evident that the first sum is a monotonically increasing function if $\Delta p_i > 0$, while the second one is a monotonically increasing function if $\Delta p_j < 0$. Thus, it is possible to find the optimal solution, i.e. the maximum, $\alpha_B^*$, as:

$$ \alpha_B^* = \sum_{i \in A^+} a_i^+ \Delta p_{\max,i} + \sum_{j \in A^-} a_j^- (-\Delta p_{\max,j}). \quad (15) $$

Indeed, the maximum contribution of each $\Delta p_i$ and $\Delta p_j$ is obtained at the boundary of the constraint, i.e. at one of the vertices of the hyperbox. Lastly, it is possible to combine the two terms in (15):

$$ \alpha_B^* = \sum_{j=1}^{n_p} |a_j| \Delta p_{\max,j} = \sum_{j=1}^{n_p} \sum_{i=1}^{n_x} n_i \Pi_{ij} |\Delta p_{\max,j}| \quad (16) \\ = \|n^\top \Pi W \Pi^\top\|_1, $$

where the definition of 1-norm has been used. □

Despite its simplicity, the presence of the absolute value makes equation (8) non-differentiable, thus not suitable for gradient-based trajectory optimization algorithms (but it can be useful for gradient-free methods such as, e.g., MPPI). It is worth underlining that, for the hyperbox case, one may

avoid directly computing the tube radius and may exploit equation (13) to build a set of $2^{n_p}$ constraints. However, this exponentially increases the number of constraints, with a consequent increase in computational complexity.

## B. Ellipsoid

One way to obtain a differentiable expression for $\alpha$ is to exploit ellipsoids. This method, firstly proposed in [11], is here revisited and extended with a formal proof, which was lacking in [11].

The parameter space ellipsoid centered at $p_n$ and inscribed<fnref n="3" /> in (6) is defined as:

$$ \mathcal{E}_p = \{ \Delta p \in \mathbb{R}^{n_p} : \Delta p^\top W^{-1} \Delta p \le 1 \} \simeq \mathcal{P}. \quad (17) $$

Exploiting the sensitivity matrix $\Pi$, it is possible to map (17) into the corresponding ellipsoid in the state space. For convenience, a proper linear scaling is first applied:

$$ \sigma_p = W^{-\frac{1}{2}} \Delta p \quad (18) $$

$$ P = \Pi W \frac{1}{2}, \quad (19) $$

such that the parameter ellipsoid defined in (17) becomes a sphere in the scaled parameter space:

$$ \mathcal{S}_p = \{ \sigma_p \in \mathbb{R}^{n_p} : \sigma_p^\top \sigma_p \le 1 \}. \quad (20) $$

We underline that the linearization in (7) holds also for the scaled space and scaled sensitivity. Then, combining (18)-(19) leads to:

$$ \Delta x \approx \Pi \Delta p = P \sigma_p, \quad (21) $$

which can be pseudoinverted to compute:

$$ \sigma_p = P^\top \Delta x + N_P w, \quad (22) $$

where $N_P \in \mathbb{R}^{n_p \times n_p}$ is the nullspace projector of the scaled sensitivity matrix $P$, and $w \in \mathbb{R}^{n_p}$ is a general solution projected in the nullspace.

By plugging (22) in (20), one can obtain the corresponding ellipsoid in the state space. Two possible cases arise depending on the dimensions of matrix $P$. When $P$ is tall or square, i.e. $n_p \le n_x$, $N_P = 0$, leading to a unique solution. However, if $P$ is fat, i.e. $n_p > n_x$, $N_P \ne 0$. Therefore, the parameter space ellipsoid can be mapped into an infinite number of state space ellipsoids, depending on $w$.

**Proposition 2.** The maximum state space ellipsoid inscribed in the hyperbox is defined as:

$$ \mathcal{E}_x = \{ \Delta x \in \mathbb{R}^{n_x} : \Delta x^\top K_{\Pi}^\dagger \Delta x \le 1 \}, \quad (23) $$

where $K_{\Pi} = \Pi W \Pi^\top$ represents the kernel matrix.

*Proof.* By plugging (22) in (20), one can obtain all the corresponding ellipsoids in the state space

$$ (P^\top \Delta x + N_P w)^\top (P^\top \Delta x + N_P w) \le 1. \quad (24) $$

Remembering that $P^\top N_P = 0$, since $P^\top$ maps any vector to the column space of $P$ while $N_P$ annihilates all components in the column space, equation (24) becomes:

$$ \Delta x^\top (P P^\top)^\top \Delta x + w^\top N_P^\top N_P w \le 1. \quad (25) $$

<footnote n="3">³Although this leads to a more conservative result, the equations can also be directly applied to the circumscribed ellipsoid.</footnote>
Since $w^\top N_P^\top N_P w \ge 0 \quad \forall w \in \mathbb{R}^{n_p}$, it is clear that by choosing any $w \notin \ker(P)$, the volume of the ellipsoid would be reduced. Therefore, the largest ellipsoid is given by

$$ \Delta x^\top K_\Pi^\dagger \Delta x \le 1, \quad (26) $$

with $K_\Pi = PP^\top = \Pi W \Pi^\top$. $\square$

Regarding the tube radius, as detailed in [17], the maximum deviation of the ellipsoid with a kernel matrix $K_\Pi$ along a generic direction $n$ can be computed as:

$$ \alpha_{\mathcal{Q}} = \sqrt{n^{\top} K_{\Pi} n}, \quad (27) $$

which coincides with equation (9).

### C. Superquadrics

As discussed previously, as the number of parameters increases, the ellipsoid covers a progressively smaller fraction of the hyperbox. To avoid this, in this paper we propose to exploit the superquadrics, a generalization of quadrics. In particular, among all the superquadrics a suitable approximation is given by the superellipsoids, defined by the following equation<fnref n="4" />:

$$ \mathcal{Q}_p = \left\{ \Delta p \in \mathbb{R}^{n_p} : \sum_{i=1}^{n_p} \left( \frac{\Delta p_i}{\Delta p_{\max,i}} \right)^{2k} \le 1 \right\}, \quad (28) $$

where $2k$ is the degree of the superellipsoid, with $k \ge 1$. In particular, as $k$ increases, the superellipsoid more closely approximates the shape of the hyperbox, especially around its corners; in contrast, when $k=1$, equation (28) reduces to (17), yielding a regular ellipsoid. For simplicity, it is possible to express (28) in a vector form as:

$$ \mathcal{Q}_p = \left\{ \Delta p \in \mathbb{R}^{n_p} : \Delta p^{\circ k^\top} W^{-k} \Delta p^{\circ k} \le 1 \right\} \approx \mathcal{P}, \quad (29) $$

where, as previously detailed, the superscript $\circ$ denotes the element-wise power.

As in the ellipsoid case, we apply a linear scaling such that the superquadric becomes a $2k$-norm ball [18]

$$ S_p^k = \left\{ \sigma_p \in \mathbb{R}^{n_p} : \sigma_p^{\circ k^\top} \sigma_p^{\circ k} \le 1 \right\}, \quad (30) $$

which is then mapped into state space by exploiting the sensitivity matrix. According to (22), the general solution to (21) is $\sigma_p = (P^\dagger \Delta x + N_P w)^\circ$, from which one can obtain the state space superquadric:

$$ \mathcal{Q}_x = \left\{ \Delta x \in \mathbb{R}^{n_x} : (P^\dagger \Delta x + N_P w)^{\circ k^\top} (P^\dagger \Delta x + N_P w)^{\circ k} \le 1 \right\}. \quad (31) $$

Due to the element-wise power, it is easy to demonstrate that $w \in \ker(P)$ is not the solution that maximizes the volume of the superquadric. This makes it not possible to explicitly define the state space superquadric a priori, depending instead on the unknown $w$, and use it directly to construct tubes. Nonetheless, having defined the tube radius as the maximum directional deviation (cfr. Sect.III), an explicit expression for $\mathcal{Q}_x$ is not needed, as the tube radius $\alpha$ can be computed as follows.

**Proposition 3.** The maximum deviation induced by the superquadric along a given direction $n$ is

$$ \alpha_{\mathcal{Q}} = \sqrt[\frac{2k}{2k-1}]{n^{\top} P (P^{\top} n)^{\circ \frac{1}{2k-1}}}. \quad (32) $$

Furthermore, as expected, $\alpha_{\mathcal{Q}}$ converges to (8) as $k \to \infty$.

*Proof.* The goal is to obtain the maximum $\alpha_{\mathcal{Q}}$ such that

$$ n^{\top} \Delta x = n^{\top} P \sigma_p = \alpha_{\mathcal{Q}} \wedge \sigma_p \in S_p^k. \quad (33) $$

Since the $2k$-norm ball is a convex set, it is possible to focus on the boundary $\partial S_p^k$. Thus, the optimal solution can be found by solving the following optimization problem:

$$ \begin{array}{ll} \alpha_{\mathcal{Q}}^* = \max_{\sigma_p} & n^{\top} P \sigma_p \\ \text{s.t.} & \sigma_p^{\circ k^\top} \sigma_p^{\circ k} = 1. \end{array} \quad (34) $$

As only equality constraints are present, problem (34) can be solved exploiting the Lagrangian Multipliers method. Firstly, the Lagrangian is defined as

$$ \mathcal{L}(\sigma_p, \lambda) = -n^{\top} P \sigma_p + \lambda (\sigma_p^{\circ k^\top} \sigma_p^{\circ k} - 1), \quad (35) $$

and its gradient is

$$ \begin{bmatrix} \nabla_{\sigma_p} \mathcal{L} \\ \nabla_{\lambda} \mathcal{L} \end{bmatrix} = \begin{bmatrix} -P^{\top} n + 2k\lambda \sigma_p^{\circ 2k-1} \\ \sigma_p^{\circ k^\top} \sigma_p^{\circ k} - 1 \end{bmatrix}. \quad (36) $$

At this point, we compute the optimal $\alpha^* = n^{\top} \Delta x^*$ by imposing $\nabla \mathcal{L} = 0$, where for ease of notation $\sigma_p^* = \sigma_p$. Therefore:

$$ \left\{ \begin{array}{l} 2k\lambda\sigma_p^{\circ 2k-1} = P^{\top} n, \\ \sigma_p^{\circ k^\top} \sigma_p^{\circ k} = 1. \end{array} \right. \quad (37) $$

$$ \sigma_p^{\circ k^\top} \sigma_p^{\circ k} = 1. \quad (38) $$

Left-multiplying (37) for $\sigma_p^\top$, and noting that $\sigma_p^\top \sigma_p^{\circ 2k-1} = \sigma_p^{\circ k^\top} \sigma_p^{\circ k}$, by substituting (38) one obtains

$$ \lambda^* = \frac{\sigma_p^\top P^\top n}{2k} = \frac{n^{\top} P \sigma_p}{2k} = \frac{\alpha_{\mathcal{Q}}^*}{2k}, \quad (39) $$

which plugged back in (37) yields

$$ \alpha_{\mathcal{Q}}^* \sigma_p^{\circ 2k-1} = P^{\top} n. \quad (40) $$

By transposing and applying the $(2k-1)$-th root to equation (40), one can obtain:

$$ \alpha_{\mathcal{Q}}^* \frac{1}{2k-1} \sigma_p^\top = (n^\top P)^{\circ \frac{1}{2k-1}}. \quad (41) $$

Finally, multiplying (40) by (41) results in

$$ \alpha_{\mathcal{Q}}^* \frac{2k}{2k-1} \sigma_p^\top \sigma_p^{\circ 2k-1} = (n^\top P)^{\circ \frac{1}{2k-1}} P^\top n, \quad (42) $$

which, recalling (38) and applying the $2k/(2k-1)$-th root, leads to the optimal solution $\alpha_{\mathcal{Q}}^*$:

$$ \begin{align} \alpha_{\mathcal{Q}}^* &= \sqrt[\frac{2k}{2k-1}]{(n^{\top} P)^{\circ \frac{1}{2k-1}} P^{\top} n} \tag{43} \\ &= \sqrt[\frac{2k}{2k-1}]{n^{\top} P (P^{\top} n)^{\circ \frac{1}{2k-1}}}. \end{align} $$

It is worth highlighting that equation (43) returns exactly the ellipsoid solution in (9) when choosing $k=1$, while it converges to the one of the hyperbox (8) when $k \to \infty$. Indeed:

$$ \lim_{k \to \infty} \sqrt[\frac{2k}{2k-1}]{n^{\top} P (P^{\top} n)^{\circ \frac{1}{2k-1}}} = n^{\top} P S = \|n^{\top} \Pi W^{\frac{1}{2}}\|_{1}, \quad (44) $$

<footnote n="4">⁴An alternative definition can be: $\sum_{i=1}^{n_p} \left| \frac{\Delta p_i}{\Delta p_{\max,i}} \right|^k \le 1$</footnote>
where $S = \{\text{sign}(P^\top n)_i\} \in \mathbb{R}^{n_p}$. $\square$

Therefore, by increasing the order $k$ it is possible to obtain an expression that is still differentiable while reducing portion of unoccupied volume of the original hyperbox. Crucially, by computing the volume ratio of the superquadric and hyperbox [18], it is possible to choose $k$ to meet a desired box covering ratio $r \in (0, 1)$ via

$$ \frac{V_Q}{V_B} = \frac{\Gamma\left(1 + \frac{1}{2k}\right)^{n_p}}{\Gamma\left(1 + \frac{n_p}{2k}\right)} \geq r. \quad (45) $$

## IV. EVALUATION

The proposed methods have been validated in two different scenarios via extensive numerical simulations that reproduce the effects of different parameter variations. The first case study involves a mass linearly moving freely in 3D space, while the second case study considers the dynamic behavior of a UR10e manipulator. In both cases, the sensitivity has been exploited to build an optimization problem computing a trajectory for point-to-point motion in configuration space. The optimization aims at minimizing the state sensitivity at the final time for improving accuracy and it leverages the input tubes to ensure that the actual (perturbed) inputs remain within the actuator limits, avoiding unwanted saturations which would result in a deviation from the planned motion. The optimization problem is formulated as:

$$ \begin{array}{ll} \min_{a} & \|\Pi(a, t_f) W \Pi(a, t_f)^\top\|_{\gamma} \\ \text{s.t.} & A_{\text{eq}}a = b_{\text{eq}}, \quad (46) \\ & u_{\text{min}} + \alpha_{\Theta}(a) \le u(a) \le u_{\text{max}} - \alpha_{\Theta}(a), \end{array} $$

where $\|\cdot\|_{\gamma}$ represents a generic matrix norm, e.g., 2-norm or Frobenius norm. $a \in \mathbb{R}^{n_a}$ are the coefficients that parametrize the trajectory, i.e. the coefficients of a polynomial parametrization, while $A_{\text{eq}} \in \mathbb{R}^{n_{\text{eq}} \times n_a}$ and $b_{\text{eq}} \in \mathbb{R}^{n_a}$ are the matrix and the vector encoding the boundary conditions, respectively. Lastly, $u_{\text{min}}$ and $u_{\text{max}}$ are the upper and lower input bounds, and $\alpha_{\Theta}$ is the stack of tube radii computed with the input sensitivity matrix $\Theta$ in the direction of the $i$-th unit vector of the canonical basis $e_i$, $\forall i \in [1, n_u]$.

All the simulations have been performed in Matlab and can be found at the following Github repository<fnref n="5" />.

### A. 3D Mass

The goal of this simulation is to present a comparison between the two methods proposed in this paper, i.e. the tube computed with the hyperbox and the superquadric, showing that having a non-differentiable constraint may lead to failure of the optimization process. The optimization has been obtained exploiting CasADi and Ipopt solver.

Defining $z \in \mathbb{R}^3$ and $\dot{z} \in \mathbb{R}^3$ as the cartesian position and velocity of the mass respectively, and taking $x = [z \ \dot{z}]^\top$ as the system state, the 3D mass can be modeled as

$$ f = \left[ \begin{matrix} \dot{z} \\ m^{-1}(\mathbf{u} - b\dot{z}) \end{matrix} \right] \quad (47) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3: Evolution of the input for the 3D Mass with the superquadric. The yellow lines represent the sensitivity tube, while the red lines illustrates the upper and lower bound of the input.</figcaption></figure>

with $m \in \mathbb{R}$ and $b \in \mathbb{R}$ being the mass and the damping of the system, respectively, which are taken as the uncertain parameters.

To control the system, a PD + Friction compensation controller is employed:

$$ u = m_n \left( \ddot{z}_{\text{des}} + K_p \ddot{z} + K_d \dot{\dot{z}} \right) + b_n \dot{z}. \quad (48) $$

Here, $m_n \in \mathbb{R}$ and $b_n \in \mathbb{R}$ are the nominal mass and damping, and $K_p \in \mathbb{R}^{3 \times 3}$ and $K_d \in \mathbb{R}^{3 \times 3}$ are the positive definite diagonal proportional and derivative gain matrices, respectively. The terms $\ddot{z} = z_{\text{des}} - z$ and $\dot{\dot{z}} = \dot{z}_{\text{des}} - \dot{z}$ represent the error and its derivative, with $z_{\text{des}}$, $\dot{z}_{\text{des}}$, $\ddot{z}_{\text{des}} \in \mathbb{R}^3$ being the desired position, velocity and acceleration of the trajectory.

The mass has to move from an initial configuration $z_0 = [0, 0, 0, 0]^\top$ to a desired final configuration $z_f = [1.0, 3.0, -1.0]^\top$ in $t_f = 5$ s, while both the initial and final velocity and acceleration are zero, minimizing the square of the Frobenius norm of the sensitivity at the end, i.e. $\gamma = F$. The nominal parameters are equal to $p_n = [m_n, b_n]^\top = [10.0, 5.0]^\top$, while the real ones may vary within a range of 25% around their nominal values. The input limits are assumed equal in all the directions and symmetric, i.e. $u_{\text{max},i} = u_{\text{min},i} = 16.0$ N. The optimization problem in (46) is solved computing $\alpha$ with the two different approaches proposed, i.e. the hyperbox defined in (8) and the superquadric defined in (10) with order $k = 50$.

Using the hyperbox constraint, the optimization problem is not capable of converging after 10,000 iterations. However, it still finds a feasible suboptimal solution with a cost function value of $\|\Pi(t_f)W\Pi(t_f)^\top\|_F = 4.9\text{e}-5$. With the superquadric constraint, instead, the optimization problem converges in only 31 iterations, obtaining a better solution of $\|\Pi(t_f)W\Pi(t_f)^\top\|_F = 3.5\text{e}-8$. It is worth underlining that the solution obtained with the superquadric produces a lower cost while still being feasible for the hyperbox constraint. This means that, for this case, theoretically the optimization problem could converge to the same solution.

To further validate the sensitivity tubes, 100 simulations were performed with randomly varying parameters lying within the hyperbox. These simulations demonstrate that the tubes are capable of encapsulating all the perturbed state/input evolutions. In particular, Fig. 3 illustrates the input evolution for the superquadric case (the hyperbox is similar), showing

<footnote n="5">⁵https://github.com/apupa/sensitivity_tubes.git</footnote>
that the input never violates the constraint.

## B. Dynamic Manipulator

The second simulation involves the dynamic model of a UR10e, a 6 DoF collaborative manipulator. The goal of this simulation is twofold. From one hand, it shows that the sensitivity approach can be applied also to more complex dynamics the full UR10e dynamic model. On the other hand, the simulation shows that the sensitivity tube computed through the superquadric is capable of encapsulating almost all the system evolution, outperforming the ellipsoid approximation exploited in the previous works. The optimization has been performed with fmincon using the interior-point algorithm. Let $x = [q \dot{q}]^T$ be the state of the system, with $q$ being the joint configuration. The robot is modeled exploiting the Euler-Lagrange formulation:

$$f = \begin{bmatrix} \dot{q} \\ \ddot{q} \end{bmatrix} = \left[ M(q,p)^{-1} (\tau - C(q, \dot{q}, p)\dot{q} - g(q,p)) \right]. \quad (49)$$

where $M(q,p) \in \mathbb{R}^{6 \times 6}$ represents the inertia matrix at configuration $q$, $C(q, \dot{q}, p) \in \mathbb{R}^{6 \times 6}$ is the Coriolis matrix, and $g(q,p) \in \mathbb{R}^6$ is the gravity vector. Lastly, $\tau \in \mathbb{R}^6$ are the torque inputs. The set of parameters that are considered uncertain are the mass of each link. Since the robot manipulator has to follow a desired trajectory in the joint space, the following control law is employed:

$$\tau = M(q,p_n) (\ddot{q}_{des} + K_p \tilde{q} + K_d \dot{\tilde{q}}) + C(q, \dot{q}, p_n) \dot{q} + g(q, p_n), \quad (50)$$

where the dependency of time has been dropped for ease of notation. $\ddot{q}_{des} \in \mathbb{R}^6$ is the feedforward acceleration, while $\tilde{q} = q_{des} - q$ and $\dot{\tilde{q}} = \dot{q}_{des} - \dot{q}$ are the position and velocity error, respectively. Lastly, $K_p \in \mathbb{R}^{6 \times 6}$ and the $K_d \in \mathbb{R}^{6 \times 6}$ are the respective gains.

The robot is tasked to move from an initial configuration $q_0 = [0, 0, -\frac{3}{5}\pi, \frac{3}{5}\pi, -\frac{\pi}{2}, -\frac{\pi}{2}, 0.0]^T$ to a desired final configuration $q_f = [\frac{\pi}{2}, -\frac{3}{5}\pi, \frac{3}{5}\pi, -\frac{\pi}{2}, -\frac{\pi}{2}, 0.0]^T$ in $t_f = 1$ s, while both the initial and final velocity and acceleration are zero. The nominal mass of each link can be found on the official UR website, while it is assumed that the real ones may vary within a range of 20% of their nominal one. The input limit is assumed equal in all the directions and symmetric, i.e. $u_{max,i} = -u_{min,i} = 23$ Nm.

First, it is essential to quantify how using a superquadric approximation of the set impacts the estimation of the tube radius. This is done by defining the ratio $r_\alpha = \alpha_Q / \alpha_B$ as a function of the superquadric order $k$. Unlike the volume ratio in (45), $r_\alpha$ varies with the chosen direction $n$ and the scaled sensitivity matrix $P$, making a closed-form expression difficult to obtain. Consequently, a statistical analysis on the manipulator's input tubes is performed: for each of the 50 randomly generated trajectories, $r_\alpha$ is evaluated across all input directions along the entire path. Figure 4 shows boxplots of $r_\alpha$ for various values of $k$ (left axis), while the dotted line traces the corresponding volume ratio from (45) (right axis). It is worth underlining that the ellipsoidal approximation employed in [11], i.e. $k = 1$, underestimates the tube radius of

<figure><img src="image_5.png" /><figcaption>Fig. 4: Boxplots of the tube radius ratio (left) and the evolution of the volume ratio (right) for different values of $k$. The ellipsoidal approximation introduced in [11] clearly underestimates the maximum deviation and achieves low volume coverage.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5: Evolution of the manipulator joints (top) and inputs (bottom). The yellow lines represent the superquadric-based sensitivity tubes, while the green lines the ellipsoid-based ones. Lastly, the red lines illustrates the input upper and lower bound of the input. Only the inputs that saturate are shown for compactness, showing that the constraints are never violated.</figcaption></figure>

approximately 50%, while the volume covered is lower than 10%.

Setting the order of the superquadric $k = 10$, the optimization problem minimizes the 2-norm in (46) in 163 iterations and the required time is around 3 minutes. Also in this case, to further validate the sensitivity tubes, 500 simulations were performed with randomly varying parameters lying within the hyperbox. These simulations demonstrate that the superquadric tubes are capable of encapsulating most of the perturbed state/input evolutions, shown in Fig. 5, while the ellipsoid tubes underestimate the variation.

To further validate the effectiveness of the superquadric approximation, the same optimization problem has been solved using ellipsoids, as in [11]. Fig. 6 illustrates that, in this case, the tubes fail to fully encapsulate the actual input evolutions. Consequently, some runs incur an unplanned input saturation with a negative impact on the joint evolution, showing significant deviations outside of the tubes.

Intuitively, this underestimation of the tube radius reduces the sensitivity tubes' ability to capture the system's closed-loop evolution. For this reason, another statistical campaign was carried out using 20 randomly generated trajectories.
<figure><img src="image_7.png" /><figcaption>Fig. 6: Evolution of the manipulator joints (top) and inputs (bottom) with only the ellipsoid tubes. The solid red lines represent instances where the inputs have been saturated in order to respect the constraints. Only the inputs that saturate are shown for compactness.</figcaption></figure>

Specifically, for each trajectory, 200 simulations were performed with the parameters sampled uniformly within the hyperbox, and it was analyzed wether the perturbed inputs remains within the tube computed using both ellipsoid and superquadric of order $k = 10$. To minimize numerical errors and remove linearization errors, the perturbed evolution was considered invalid if it never goes out of the hyperbox-based tube and if, at any time step, it exceeded the ellipsoid or superquadric tube by more than 2% of the maximum tube radius. The results show that with the ellipsoid 61% of the perturbed inputs go outside the tube, i.e. they are invalid, while with the superquadric only 4%.

## V. CONCLUSIONS AND FUTURE WORK

In this work, two new ways for computing sensitivity tubes have been proposed: the first method exploits superquadrics for approximating a bounding hyperbox in parameter space with an arbitrary precision while retaining differentiability of the resulting tube radius, thus making it seamlessly usable in any gradient-based optimization algorithm. The second method, instead, exploits directly the definition of the bounding hyperbox for evaluating the maximum state deviation, which is proven to always lie at a vertex of the box, but loses full differentiability of the tube radius. For both methods, a formal proof is given for obtaining the corresponding tube radius, and the whole framework is validated in a simulation campaign on two case studies.

The proposed methods characterize the state/input deviations by exploiting the respective sensitivity matrices and assuming a linear behavior around the nominal parameters (7). Future work will focus on quantifying the linearization error and identifying an upper bound for its value. Consequently, the tube size could be increased to achieve robust guarantees. Moreover, the assumption of uniformly distributed parameter variations in a range could be replaced with other assumptions such as a Gaussian distribution centered at a nominal value. In this case, superquadrics could be designed to capture a desired probability of variations. Lastly, in real scenarios parameter

variations are often correlated, e.g., a change in mass may also affect the position of the center of mass. The proposed method could be extended to account for such dependencies as well.

## REFERENCES

[1] C. Zimmer, M. Meister, and D. Nguyen-Tuong, "Safe active learning for time-series modeling with gaussian processes," in *Advances in Neural Information Processing Systems*, vol. 31, 2018.

[2] K. Astrom, "Adaptive feedback control," *Proceedings of the IEEE*, vol. 75, no. 2, pp. 185-217, 1987.

[3] K. Zhou and J. C. Doyle, *Essentials of robust control*. Prentice hall Upper Saddle River, NJ, 1998, vol. 104.

[4] Y. Shtessel, C. Edwards, L. Fridman, A. Levant et al., *Sliding mode control and observation*. Springer, 2014, vol. 10.

[5] R. Lozano, B. Brogliato, O. Egeland, and B. Maschke, *Passivity-Based Control*. London: Springer London, 2000, pp. 227-277.

[6] A. van der Schaft, *L2-Gain and Passivity Techniques in Nonlinear Control*, 3rd ed. Springer Publishing Company, Incorporated, 2016.

[7] P. Robuffo Giordano, Q. Delamare, and A. Franchi, "Trajectory generation for minimum closed-loop state sensitivity," in *2018 IEEE International Conference on Robotics and Automation (ICRA)*. IEEE, 2018, pp. 286-293.

[8] A. Srour, S. Marcellini, T. Belvedere, M. Cognetti, A. Franchi, and P. Robuffo Giordano, "Experimental validation of sensitivity-aware trajectory planning for a quadrotor uav under parametric uncertainty," in *2024 International Conference on Unmanned Aircraft Systems (ICUAS)*. IEEE, Jun. 2024, p. 572-578.

[9] A. Pupa, P. Robuffo Giordano, and C. Secchi, "Optimal energy tank initialization for minimum sensitivity to model uncertainties," in *2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*. IEEE, 2023, pp. 8192-8199.

[10] S. Wasiela, P. Robuffo Giordano, J. Cortés, and T. Simeon, "A sensitivity-aware motion planner (samp) to generate intrinsically-robust trajectories," in *2023 IEEE International Conference on Robotics and Automation*, 2023.

[11] A. Afifi, T. Belvedere, A. Pupa, P. Robuffo Giordano, and A. Franchi, "Safe and robust planning for uncertain robots: A closed-loop state sensitivity approach," *IEEE Robotics and Automation Letters*, vol. 9, no. 11, p. 9962-9969, Nov. 2024.

[12] T. Belvedere, M. Cognetti, G. Oriolo, and P. Robuffo Giordano, "Sensitivity-aware model predictive control for robots with parametric uncertainty," *IEEE Transactions on Robotics*, vol. 41, pp. 3039-3058, 2025.

[13] A. H. Barr, "Superquadrics and angle-preserving transformations," *IEEE Computer Graphics and Applications*, vol. 1, no. 01, pp. 11-23, 1981.

[14] G. Vezzani, U. Pattacini, and L. Natale, "A grasping approach based on superquadric models," in *2017 IEEE International Conference on Robotics and Automation (ICRA)*. IEEE, 2017, pp. 1579-1586.

[15] F. Messerer and M. Diehl, "An efficient algorithm for tube-based robust nonlinear optimal control with optimal linear feedback," in *2021 60th IEEE Conference on Decision and Control (CDC)*. IEEE, 2021, pp. 6714-6721.

[16] R. Reams, "Hadamard inverses, square roots and products of almost semidefinite matrices," *Linear Algebra and its Applications*, vol. 288, pp. 35-43, 1999.

[17] D. H. Van Hessem, C. W. Scherer, and O. H. Bosgra, "Lmi-based closed-loop economic optimization of stochastic process operation under state and input constraints," in *2001 40th IEEE Conference on Decision and Control (CDC)*. IEEE, 2001, pp. 4228-4233.

[18] A. Slavík and D. Šulc, "Maximal volumes of n-dimensional balls in the p-norm," *Archiv der Mathematik*, vol. 114, no. 3, pp. 305-312, 2020.