# Contact-Robust Trajectory Planning via Parametric Sensitivity Analysis for Hybrid Robotic Systems

Tommaso Belvedere, James Zhu, Marco Cognetti, Paolo Robuffo Giordano

► To cite this version:

Tommaso Belvedere, James Zhu, Marco Cognetti, Paolo Robuffo Giordano. Contact-Robust Trajectory Planning via Parametric Sensitivity Analysis for Hybrid Robotic Systems. ICRA 2026 - IEEE International Conference on Robotics and Automation, Jun 2026, Vienne, Austria. hal-05488973v2

HAL Id: hal-05488973

https://hal.science/hal-05488973v2

Submitted on 6 Feb 2026

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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 CC BY 4.0 - Attribution - International License
# Contact-Robust Trajectory Planning via Parametric Sensitivity Analysis for Hybrid Robotic Systems

Tommaso Belvedere*, James Zhu*, Marco Cognetti, Paolo Robuffo Giordano

**Abstract**— In this paper, we combine first-order approximations of hybrid systems (i.e., the so-called saltation matrix) with previous works on parametric sensitivity for continuous systems to propose a general framework for robust trajectory optimization of hybrid systems subject to parametric uncertainties. A method for computing parametric sensitivities of both continuous dynamics and hybrid events is presented. The obtained “hybrid parametric sensitivity” is then combined with sensitivity-based tubes that encapsulate all possible perturbed states and control trajectories given a known bounded range for the uncertain parameters. The proposed method is then applied to the problem of planning robust trajectories for legged robot systems, which allows obtaining trajectories that remain feasible w.r.t. the contact constraints even in presence of uncertainties in the dynamics, guard conditions, and reset maps. We also illustrate one of the fundamental limitations of first-order approximations, that is, the fact that the sensitivity reset time is fixed, and propose an extension to the sensitivity analysis that can form the basis for future developments.

## I. INTRODUCTION

Impact-aware manipulation and agile locomotion both rely on the ability to establish and break contact with the environment. Since such events occur at non-zero velocity, reaction forces must be carefully regulated to stabilize systems across different contact modes. This is particularly challenging due to the short duration of contact events and the limitations imposed by friction, with unintended contact loss or slip having severe consequences on robot performance. In any real-world scenario, all models are uncertain to some extent and, when uncertainty is not negligible, one must rely on ‘robust’ planning, control, and/or estimation algorithms. A number of works have addressed uncertainty inherent to contact dynamics. For instance, [1] deals with contact time uncertainty by adapting reference vector fields to remove control input peaks. In legged robots, stabilization of jumping and landing maneuvers has been achieved either by embedding compliance and elasticity in the controller [2], [3] or through tailored trajectory optimization approaches [4].

In this work, we address the problem of generating motion plans that are intrinsically robust to model uncertainty for hybrid systems such as those describing robots in contact with

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

J. Zhu and M. Cognetti are with the LAAS-CNRS, CNRS, Université de Toulouse, CNRS, Toulouse, France. E-mail: {james.zhu, marco.cognetti}@laas.fr

This work was supported by the project ANR-20-CE33-0003 "CAMP", by the Chaire de Professeur Junior Grant ANR-22-CPJ1-0064-01, and by the Cluster SequoIA as part of the Chair HARMONY funded by ANR and France 2030 under the reference number ANR-23-IACL-0009.

<figure><img src="image_2.png" /><figcaption>Fig. 1. For robotic systems that undergo hybrid transitions, it is critical to follow a planned contact sequence while preventing unintended transitions. In legged robots, this requires maintaining ground reaction forces positive and within the friction cone. Baseline planning methods often saturate these constraints, making the system sensitive to small deviations that can cause slips or liftoffs, potentially leading to catastrophic failure. The proposed sensitivity-aware method estimates constraint sensitivities with respect to parametric uncertainties, enabling the planning of safe and robust motions also in presence of parametric uncertainty.</figcaption></figure>

the environment. The proposed sensitivity-aware methodology produces trajectories designed to withstand variations in system parameters and prevent unintended changes in the contact mode. We specifically study the effects of parametric uncertainty, which arises from imperfect knowledge of system parameters; for legged robots this can include mass, inertia, friction, or environmental characteristics. The particular structure of parametric uncertainty can be exploited to generate accurate estimates of the evolution of uncertain systems, which then allows planning robust motions.

An effective tool to analyze uncertain dynamical systems is sensitivity analysis, which provides a local approximation of perturbed state trajectories [5], [6]. For continuous systems, sensitivity analysis has been exploited as a computationally efficient tool to robustify robot motions by generating minimally sensitive trajectories [6], or to construct tubes for safe constrained trajectory optimization [7]–[10]. Hybrid dynamical systems such as legged robots are characterized by the interplay between continuous-time dynamics and discrete events at which the state and/or the dynamics switch instantaneously, resulting in non-smooth trajectories. Thus, uncertainties not only affect the continuous-time evolution of the system within each mode, but also discrete components such as the guard conditions (defining event triggers) and the reset maps (governing post-event transitions). This interplay between continuous and discrete elements increases the challenge of analyzing and designing robust hybrid systems, as uncertainties can alter the sequence, the timing of events, and their effects on the system trajectories. Therefore, a contribution of this work is to extend the sensitivity analysis for continuous systems [6] to the case of hybrid systems with application in legged robotics.

<footnote>*These authors contributed equally to this work.</footnote>
## II. RELATED WORKS

Several works have addressed sensitivity analysis in hybrid systems from the control perspective. For example, prior results have examined the sensitivity of hybrid systems to initial conditions, highlighting how variations propagate through both continuous and discrete state jumps [11]. Most works have focused on the concept of saltation matrix, describing the differential relationship between states before and after events (see [12] for a survey). These contributions have paved the way for applications such as local stabilization through linear hybrid feedback [13] and trajectory optimization for minimally sensitive planning [14]. Extensions to traditional sensitivity tools, such as the generalized saltation matrix introduced in [15], are also able to explicitly incorporate uncertainties in guard conditions and reset maps, enabling the use of sensitivity analysis in applications like bipedal locomotion [16].

Although successful, all these works have focused on the effect of uncertainties only in the initial conditions and/or the hybrid event, which represent a limited class of problems in sensitivity analysis. Moreover, works like [14], [16] only aim at minimizing (a norm of) the saltation matrix and do not consider the (more complex) problem of robust constraint satisfaction over a future planned trajectory. In this respect, our first contribution is the development of a unified approach for computing first-order approximations of hybrid systems w.r.t. uncertain parameters affecting the complete hybrid dynamics, including equations of motion, guard conditions, and reset maps (Sect. IV). We then propose an extension to the sensitivity definition (Sect. IV-B) that, while increasing complexity in the modeling of parameter deviations, provides a more accurate approximation in the vicinity of switching events. Finally, we apply the proposed hybrid parametric sensitivity to robust trajectory planning for two common legged locomotion problems (Sect. V), generating motions that are intrinsically robust to parametric uncertainties.

## III. PRELIMINARIES

This section recaps the concepts of state sensitivity (Sect. III-A) and tubes of perturbed trajectories (Sect. III-B) for continuous systems, providing the foundation for extending these notions to hybrid systems in Sect. IV.

### A. Sensitivities of Continuous Systems

Consider a dynamical system

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

where $\mathbf{x} \in \mathbb{R}^{n_x}$ and $\mathbf{u} \in \mathbb{R}^{n_u}$ are the state and input vectors, $\mathbf{p} \in \mathbb{R}^{n_p}$ represents uncertain model parameters with nominal value $\mathbf{p}_n$, and $\mathbf{f}: \mathbb{R}^{n_x} \times \mathbb{R}^{n_u} \times \mathbb{R}^{n_p} \times \mathbb{R} \rightarrow \mathbb{R}^{n_x}$ is a smooth dynamics map. Let also

$$ \mathbf{u} = \mu(\mathbf{x}, \mathbf{y}_d(t), \mathbf{p}_n) \quad (2) $$

be a control action<fnref n="1" /> to drive system (1) to track a reference output $\mathbf{y}_d(t)$, and let $\varphi: \mathbb{R}^{n_x} \times \mathbb{R}^{n_p} \times \mathbb{R} \rightarrow \mathbb{R}^{n_x}$ be the

flow of (1-2) so that, for an initial state $\mathbf{x}(0) = \mathbf{x}_0$, the state trajectory is $\mathbf{x}(t) = \varphi(\mathbf{x}_0, \mathbf{p}, t)$.

The sensitivity matrix $\mathbf{\Pi}(t) \in \mathbb{R}^{n_x \times n_p}$, defined as

$$ \mathbf{\Pi}(t) = \left. \frac{\partial \varphi(\mathbf{x}_0, \mathbf{p}, t)}{\partial \mathbf{p}} \right|_{\mathbf{p}=\mathbf{p}_n}, \quad (3) $$

represents a first-order approximation of the effects of a parameter deviation on the state trajectory around the nominal trajectory $\mathbf{x}_n(t)$ that is obtained when $\mathbf{p} = \mathbf{p}_n$. A closed-form expression for (3) is, in general, not available but, as shown in [6], it is possible to obtain the dynamics of $\mathbf{\Pi}(t)$ in closed-form as

$$ \dot{\mathbf{\Pi}}(t) = \left( \frac{\partial \mathbf{f}}{\partial \mathbf{x}} + \frac{\partial \mathbf{f}}{\partial \mathbf{u}} \frac{\partial \mu}{\partial \mathbf{x}} \right) \mathbf{\Pi}(t) + \frac{\partial \mathbf{f}}{\partial \mathbf{p}}, \quad \mathbf{\Pi}(0) = \mathbf{0}, \quad (4) $$

where all the derivatives are evaluated at the nominal state/input trajectory at time $t$. Note that the term $\frac{\partial \mu}{\partial \mathbf{x}}$ captures the sensitivity of the feedback controller to changes in the state (for an open-loop system this term would be zero).

In case the initial state $\mathbf{x}_0$ is also uncertain, it is also possible to define the closed-loop state sensitivity w.r.t. $\mathbf{x}_0$ as $\mathbf{\Psi}(t) = \frac{\partial}{\partial \mathbf{x}_0} \varphi(\mathbf{x}_0, \mathbf{p}, t)$ and obtain a closed-form expression for its dynamics as discussed in [7].

### B. Tubes of Perturbed Trajectories

One use of sensitivity analysis is to construct reference trajectories $\mathbf{y}_d(t)$ that, when tracked, intrinsically ensure the feasibility of constraints (thus safety) even in the presence of parametric uncertainty. In [10], the state sensitivity is used to build an envelope, called a tube, centered at the nominal state trajectory and containing all possible perturbed state evolutions for a given range of parametric variation. The trajectory can then be optimized to enforce constraint satisfaction even in presence of perturbations by modifying the constraints with suitable margins obtained from the corresponding tubes.

Assume that the uncertain parameters $\mathbf{p}$ belong to an ellipsoid centered at the nominal value $\mathbf{p}_n$:

$$ \mathcal{E}_{\mathbf{p}} = \{ \mathbf{p} \in \mathbb{R}^{n_p} : (\mathbf{p} - \mathbf{p}_n)^T \mathbf{W} (\mathbf{p} - \mathbf{p}_n) \le 1 \} \quad (5) $$

with $\mathbf{W} = \text{diag}\{\Delta \mathbf{p}_{\text{max},i}\}, i = 0, ..., n_p$, and $\Delta \mathbf{p}_{\text{max},i}$ being a maximum deviation for the $i$-th parameter. The maximum predicted deviation of any scalar function $c(\mathbf{x}, \mathbf{p}, t)$ (i.e. the tube radius in the direction of the constraint) at a given time $t$ can be expressed as follows:

$$ \mathbf{\Pi}_c(t) = \frac{\partial c}{\partial \mathbf{x}} \mathbf{\Pi}(t) + \frac{\partial c}{\partial \mathbf{p}} \quad (6) $$

$$ \rho_c(t) = \sqrt{\mathbf{\Pi}_c(t) \mathbf{W} \mathbf{\Pi}_c(t)^T}. \quad (7) $$

Note that since $c$ is a scalar function, $\mathbf{\Pi}_c$ has dimension $1 \times n_p$, yielding a scalar radius $\rho_c(t)$. The tube radius (7) can be evaluated along the whole trajectory, which allows obtaining (for a particular constraint $c$) the envelope of perturbed trajectories w.r.t. the nominal case under the parametric uncertainty model (5). This tube computation has already been successfully exploited in previous works, e.g., for planning a robust collision-free trajectory [7] or for reducing undesired input saturations due to model uncertainty [17].

<footnote n="1">¹For simplicity, the formulation presented here considers static controllers, but it is straightforward to extend the derivations to controllers with internal states and their sensitivity as shown in, e.g., [7].</footnote>
IV. SENSITIVITIES OF HYBRID SYSTEMS

In the following, we extend the notion of parametric state sensitivity to the case of hybrid systems with state-triggered jumps. To simplify the discussion and maintain clarity in notation, we focus on the case of a single jump between two system modes. However by following [12], our approach can be readily extended to more general scenarios [18].

Consider a smooth guard condition $g(x, p, t)$: at time $\tau$, an event is triggered by $g(x(\tau), p, \tau) = 0$ after which the smooth reset map $T(x, p, \tau)$ is applied to the state. We denote ante- and post- quantities before and after the event, respectively. Without loss of generality, we will assume that a feedback control law (2) is applied to the dynamics. Moreover, we remove explicit dependence on the input $u$ for sake of clarity. Therefore, the closed-loop system is governed by the two, possibly different, ante- and post- Lipschitz continuous dynamics

$$ 
\begin{aligned}
\dot{x}^a(t) &= f^a(x, p, t) \quad t \in [t_0, \tau) \\
\dot{x}^p(t) &= f^p(x, p, t) \quad t \in [\tau, t_f].
\end{aligned}
\quad (8) $$

Let $x(t)$, $t \in [t_0, t_f]$, be the closed-loop state trajectory defined as

$$ x(t) = \begin{cases} x^a(t_0) = x_0 \\ x^a(t) = \varphi^a(x_0, p, t) \quad t \in [t_0, \tau) \\ x^p(\tau) = T(x^a(\tau), p, \tau) \\ x^p(t) = \varphi^p(x^p(\tau), p, t) \quad t \in [\tau, t_f], \end{cases} \quad (9) $$

where we will refer to $x_n(t)$ (obtained when $p = p_n$) as the nominal trajectory. For reasons that will be clarified later, assume that the transversality condition

$$ \frac{\partial g}{\partial t} + \frac{\partial g}{\partial x} f^a \neq 0 \quad (10) $$

on the guard condition holds at $(x_n^a(\tau), p_n, \tau)$. If the system parameters are perturbed, that is $p = \bar{p} = p_n + \Delta p$, the state reset might be triggered at a different time $\bar{\tau} = \tau + \Delta\tau$ determined by the guard condition $g(\bar{x}^a(\bar{\tau}), \bar{p}, \bar{\tau}) = 0$. As in [15], assume that the dynamics vector fields $f^a, f^p$ are defined in sets $D^a, D^b \subseteq \mathbb{R}^{n_x}$ that extend beyond the nominal guard condition $g(x_n(\tau), p_n, \tau) = 0$, so that all expressions are well defined also in the perturbed case. The perturbed trajectory can then be described as:

$$ \bar{x}(t) = \begin{cases} \bar{x}^a(t_0) = x_0 \\ \bar{x}^a(t) = \varphi^a(x_0, \bar{p}, t) \quad t \in [t_0, \bar{\tau}) \\ \bar{x}^p(\bar{\tau}) = T(\bar{x}^a(\bar{\tau}), \bar{p}, \bar{\tau}) \\ \bar{x}^p(t) = \varphi^p(\bar{x}^p(\bar{\tau}), \bar{p}, t) \quad t \in [\bar{\tau}, t_f]. \end{cases} \quad (11) $$

Similarly to (3), we seek a linear approximation of the closed-loop state trajectory around the nominal parameters, which allows exploiting the first-order approximation $\bar{x}(t) \simeq x_n(t) + \Pi(t)\Delta p$ for small enough $\Delta p$.

**Proposition 1.** Consider the hybrid state trajectory $x(t)$ in (9). Under the above-mentioned assumptions, its hybrid state sensitivity w.r.t. uncertain parameters $p$ can be expressed as:

$$ \Pi(t) = \begin{cases} \Pi^a(t_0) = 0 \\ \Pi^a(t) = \frac{\partial \varphi^a(x_0, p, t)}{\partial p} \quad t \in [t_0, \tau) \\ \Pi^p(\tau) = \Pi_0^p \\ \Pi^p(t) = \frac{\partial \varphi^p(x^p(\tau), p, t)}{\partial p} \quad t \in [\tau, t_f] \end{cases} \quad (12) $$

with

$$ \Pi_0^p = \left( f^p - \frac{\partial T}{\partial x} f^a - \frac{\partial T}{\partial t} \right) \frac{\partial g}{\partial p} + \frac{\partial g}{\partial x} \Pi^a + \frac{\partial T}{\partial x} \Pi^a + \frac{\partial T}{\partial p} \quad (13) $$

where all derivatives are evaluated at $(x_n^a(\tau), p_n, \tau)$, $f^a := f^a(x_n^a(\tau), p_n, \tau)$, and $f^p := f^p(x_n^p(\tau), p_n, \tau)$.

*Proof.* In the two separate time intervals $[t_0, \tau)$, $(\tau, t_f]$, the system is governed by the continuous dynamics (8). Thus, the state sensitivity can be computed as illustrated in Sect. III-A. What remains to be found is the computation of the sensitivity reset at time $t = \tau$, which resets the initial condition for the sensitivity in $(\tau, t_f]$. Thanks to the transversality condition (10), it is possible to apply the Implicit Function Theorem to consider $\tau$ as an implicit function of the state and, as a byproduct, of the parameters $g(x_n^a(\tau), p, \tau(p))$. By expanding the guard condition around the nominal trajectory, one has

$$ 0 = \frac{\partial g}{\partial x} \left( \frac{\partial \varphi^a}{\partial t} \Delta \tau + \frac{\partial \varphi^a}{\partial p} \Delta p \right) + \frac{\partial g}{\partial t} \Delta \tau + \frac{\partial g}{\partial p} \Delta p. \quad (14) $$

We can then solve for $\Delta\tau$, as long as the transversality condition (10) holds:

$$ \Delta\tau = - \left( \frac{\partial g}{\partial t} + \frac{\partial g}{\partial x} f^a \right)^{-1} \left( \frac{\partial g}{\partial x} \Pi^a + \frac{\partial g}{\partial p} \right) \Delta p. \quad (15) $$

Then, the first-order expansion of the post-event state yields

$$ \begin{aligned} \bar{x}^p(\bar{\tau}) &\approx x_n^p(\bar{\tau}) + \frac{\partial \varphi^p(\bar{\tau})}{\partial p} \Delta p + \frac{\partial \varphi^p(\bar{\tau})}{\partial t} \Delta \tau \\ &\approx x_n^p(\bar{\tau}) + \Pi_0^p \Delta p + f^p(x_n^p(\bar{\tau}), p_n, \bar{\tau}) \Delta \tau. \end{aligned} \quad (16) $$

At the same time, from (11), $x^p(\tau)$ depends on the ante-event state through the reset map. Thus, it can also be expanded as

$$ \bar{x}^p(\bar{\tau}) \approx x_n^p(\bar{\tau}) + \frac{\partial T}{\partial x} (f^a \Delta\tau + \Pi^a \Delta p) + \frac{\partial T}{\partial t} \Delta\tau + \frac{\partial T}{\partial p} \Delta p. \quad (17) $$

Finally, by equating (16) and (17), and substituting (15), one obtains (13). $\square$

**Remark 1.** The only difference between parametric state sensitivity (12) for hybrid systems and the continuous case presented in Sect. III-A lies in the reset map (13) at time $t = \tau$. For this reason, the hybrid state sensitivity can also be leveraged for robust trajectory optimization, for instance, by minimizing a norm of the sensitivity or by imposing robust constraints through sensitivity-based tubes. Note that, while the obtained expressions account for the perturbation of the switching time, the sensitivity reset still happens at the nominal time $\tau$ rather than at the perturbed time $\bar{\tau}$, which induces a structural approximation error. An extended discussion on this point is provided in Sect. IV-B.

**Remark 2.** It can be shown that the saltation matrix [12] can be reformulated in terms of the sensitivity reset (13) if the parametric sensitivity $\Pi(t)$ is replaced by the sensitivity w.r.t. the initial conditions $\Psi(t)$. Moreover, following [7], it is also possible to combine these sensitivities to obtain a general framework that accounts for uncertainties in both model parameters and initial states.
## A. Illustrative Example

We now validate the hybrid sensitivity expression and discuss its use in approximating the evolution of a hybrid system when its parameters are perturbed. As an illustrative case study, consider a ball bouncing on a rigid surface under the effect of a scalar force $F$. Let $\mathbf{x} = (x, \dot{x})$ denote the state of the system consisting of the distance from the surface and its time derivative. The second-order robot system dynamics is $\ddot{x} = -p_1 F$, the reset map is $\mathbf{T}(\mathbf{x}, \mathbf{p}) = (x, -\dot{x}p_2)^T$, and the guard condition is $g(\mathbf{x}, \mathbf{p}) = x - p_3$. The parameters $\mathbf{p} = (p_1, p_2, p_3)$ (considered uncertain) are a multiplicative uncertainty on the force ($p_1$), the coefficient of restitution after impact ($p_2$), and the position of the rigid surface ($p_3$).

To validate the parametric sensitivity computation, we start by numerically simulating the system evolution for a nominal combination of parameters $\mathbf{p}_n = (1, 1, 0)$ and then compute the state sensitivity $\mathbf{\Pi}(t)$ associated with the nominal trajectory $\mathbf{x}_n(t)$. Then, we simulate the evolution $\bar{\mathbf{x}}(t)$ for some perturbed parameters $\bar{\mathbf{p}}$ and compare it with the first-order approximation $\bar{\mathbf{x}}(t) = \mathbf{x}_n(t) + \mathbf{\Pi}(t)\Delta\mathbf{p}$. The result of this procedure for two parametric perturbations is depicted in Fig. 2. One can note how, despite the relatively large deviation in the resulting position, the sensitivity is able to provide an accurate approximation of the perturbed trajectory for times sufficiently far from the nominal event time $\tau$. On the other hand, in a small interval before or after the event (depending on whether it is anticipated or not), there exists a discrepancy between the first-order approximation and the actual perturbed trajectories. This is an intrinsic limitation of such linear approximations [11], even at their second-order formulation [19]. In fact, these analyses are not able to deal with the split distribution arising near a hybrid event [12]. As discussed in Remark 1, this can be attributed to the fact that the actual perturbed time $\bar{\tau}$ explicitly depends on the parameter deviation $\Delta\mathbf{p}$, not known in advance, while the sensitivity reset in (13) is evaluated at the nominal time $\tau$.

## B. Parametric Sensitivity Extension

A possible way to overcome the above-mentioned limitation is to apply the reasoning of *ante* and *post* trajectory extensions [11], [13] to the sensitivity evaluation. Having assumed that the dynamics is well defined in a region that extends beyond the nominal guard condition, we can extend the nominal trajectories $\mathbf{x}_n^a(t)$ and $\mathbf{x}_n^p(t)$ in (9) for a sufficiently large $\Delta t_{\text{ext}} > 0$:

$$ 
\begin{aligned}
\mathbf{x}_{\text{ext}}^{a}(t) & \quad t \in [t_0, \tau + \Delta t_{\text{ext}}] \\
\mathbf{x}_{\text{ext}}^{p}(t) & \quad t \in [\tau - \Delta t_{\text{ext}}, t_f],
\end{aligned} 
\quad (18) $$

with boundary conditions obtained from the nominal trajectory $\mathbf{x}_{\text{ext}}^a(0) = \mathbf{x}_0$, $\mathbf{x}_{\text{ext}}^p(\tau) = \mathbf{T}(\mathbf{x}_{\text{ext}}^a(\tau), \mathbf{p}_n, \tau)$. Similarly, it is possible to extend the parametric sensitivities for the *ante-* and *post-* event trajectories:

$$ 
\begin{aligned}
\Pi_{\text{ext}}^{a}(t) & \quad t \in [t_0, \tau + \Delta t_{\text{ext}}] \\
\Pi_{\text{ext}}^{p}(t) & \quad t \in [\tau - \Delta t_{\text{ext}}, t_f]
\end{aligned} 
\quad (19) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Evolution of the ball position over time in the nominal (black) and perturbed cases (cyan), along with their first-order approximation obtained using standard sensitivity (orange) and extended sensitivity (violet). The two nominal extended trajectories are depicted with dashed gray lines. The approximate switching time $\tau + \Delta\tau$ is also shown close to the actual perturbed switching time $\bar{\tau}$. Initial condition $\mathbf{x}(0) = (0.5, 0)$, $F = 10$ N perturbed parameters $\bar{\mathbf{p}}^1 = (0.8, 0.95, -0.02)$ and $\bar{\mathbf{p}}^2 = (1.2, 1.05, 0.02)$.</figcaption></figure>

and boundary conditions $\Pi_{\text{ext}}^a(0) = 0$, $\Pi_{\text{ext}}^p(\tau) = \Pi_0^p$ from (13). Based on this, we propose the alternative approximation

$$ \bar{\mathbf{x}}(t) \approx \tilde{\mathbf{x}}_{\text{ext}}(t) = \begin{cases} \mathbf{x}_{\text{ext}}^{a}(t) + \mathbf{\Pi}_{\text{ext}}^{a}(t)\Delta\mathbf{p} & t \in [t_0, \tau + \Delta\tau) \\ \mathbf{x}_{\text{ext}}^{p}(t) + \mathbf{\Pi}_{\text{ext}}^{p}(t)\Delta\mathbf{p} & t \in [\tau + \Delta\tau, t_f] \end{cases} \quad (20) $$

with $\Delta\tau$ defined as in (15) and evaluated at the nominal switching time $\tau$.

Compared to (12), the extended expression (20) adapts the switching time (as a function of $\Delta\mathbf{p}$) instead of keeping it fixed at $\tau$. As shown in Fig. 2, this allows for much more accurate predictions than the previous linear approximation $\tilde{\mathbf{x}}(t)$, with only a minor error from estimating $\Delta\tau$ via (15). However, computing ellipsoidal tubes based on this extended sensitivity requires knowledge of the unknown parameter deviation $\Delta\mathbf{p}$, possibly requiring a switching mechanism to account for the worst-case sensitivity. Also, being discontinuous and requiring to extend the *ante-* and *post-* event trajectories, utilizing this extended sensitivity in trajectory optimization is less straightforward. For these reasons, in the following, we focus on the simpler, albeit less accurate, linear sensitivity approximation (12–13) and show how it can effectively be used in the context of robust trajectory generation. It is nevertheless worth mentioning that (20) could be applied in sensitivity-based data augmentation for efficient learning (inspired by, e.g., [20], [21]), in which approximate samples are generated through sensitivity to decrease the computational requirements of dataset generation and domain randomization in learning algorithms.
**Remark 3.** While similar in spirit to the approximation shown in [22, Appendix A], the proposed expression (20) with extended sensitivities is different as it makes use of (19) instead of performing a linear approximation of the dynamics around the time of the event, which would coincide with a 1-step Euler integration over the estimated time interval $\Delta\tau$.

## V. EXPERIMENTS & RESULTS

In this section, we leverage the hybrid sensitivity derivation of (12) to design robust trajectories for legged robots. As a case study, we focus on the stable landing problem, where a legged robot transitions from an aerial phase to a stance phase while avoiding spurious contact transitions.

Legged robots hold an advantage over wheeled platforms in their ability to execute highly dynamic maneuvers such as running, leaping, and flipping. At the same time, these capabilities expose them to greater risks of tripping, slipping, and falling. To mitigate such risks, it is essential to ensure that the robot strictly follows the intended contact sequence, with feet making contact when and only when planned. To address this challenge, we propose a robust trajectory optimization framework that suppresses undesired contact transitions under bounded parametric uncertainty. We validate the approach in simulation on two canonical models of legged locomotion: the spring-loaded inverted pendulum (SLIP) and a planar quadruped.

### A. Vertical SLIP Hopper

We first consider a vertical spring-loaded inverted pendulum (SLIP) hopper composed of a point-mass body and a massless spring leg. The leg is equipped with an actuator that can preload the spring during the aerial phase and apply additional force during stance. The SLIP provides a simple yet representative model of legged locomotion that captures the essential hybrid dynamics of flight and stance phases [23]. As such, it serves as a useful benchmark for testing whether the proposed trajectory optimization framework can suppress undesired transitions and guarantee consistent landings under parameter variations. The system state is given by $\mathbf{x} = (z, \dot{z}) \in \mathbb{R}^2$, where $z$ denotes the body height and $\dot{z}$ the vertical velocity. The model includes parametric uncertainty in the body mass $m$ and spring stiffness $k$, represented by $\mathbf{p} = (m, k)$, with nominal values $\mathbf{p}_n = (1.0, 100.0)$. The nominal spring length at rest is fixed at $\ell_0 = 0.25$ m.

The trajectory optimization problem is formulated within a hybrid direct collocation framework [24]. The aerial mode is discretized into $N$ steps of duration $\delta t^a$, with associated states $\mathbf{x}_{0:N}^a$ and inputs $\mathbf{u}_{0:N-1}^a$ representing the spring leg preload in the air. Similarly, the stance mode is discretized into $N$ steps of duration $\delta t^s$, with states $\mathbf{x}_{0:N}^s$ and inputs $\mathbf{u}_{0:N-1}^s$ representing the leg force applied in stance. The cost function to be minimized is

$$J = \sum_{k=0}^{N-1} (\mathbf{u}_k^a \mathbf{R}^a \mathbf{u}_k^{aT} + \mathbf{u}_k^s \mathbf{R}^s \mathbf{u}_k^{sT}) + \mathbf{x}_N^s \mathbf{Q}^s \mathbf{x}_N^{sT} \quad (21)$$

where $\mathbf{R}^a, \mathbf{R}^s \succeq 0$ weight the control effort in each mode, and $\mathbf{Q}^s \succeq 0$ penalizes deviation from the desired terminal state at the end of stance. The decision variables are

the control inputs $\mathbf{u}_{0:N-1}^a$ and $\mathbf{u}_{0:N-1}^s$ and the time steps $\delta t_{0:N-1}^a, \delta t_{0:N-1}^s$. By treating the time steps as optimization variables, the framework allows the transition time $\tau = (N-1)\delta t^a$ to be modulated. The state trajectory is constrained to begin from the initial condition $\mathbf{x}_0$ and to evolve according to the discretized version of the hybrid dynamics (8) in each mode, with $\mathbf{f}^a$ governing the dynamics in the aerial phase and $\mathbf{f}^s$ the stance dynamics. The sensitivity matrices $\mathbf{\Pi}_{0:N}^a, \mathbf{\Pi}_{0:N}^s$ are computed using (13).<fnref n="2" />

For this experiment, the system is initialized at a height of 0.5 m with zero vertical velocity, $\mathbf{x}_0 = (0.5, 0)$, and is tasked with finishing with zero velocity in stance. The weighting matrices are selected as $\mathbf{R}^a = \mathbf{R}^s = 0.1\mathbf{I}$ and $\mathbf{Q}^s = \text{diag}(0, 10)$. The phase duration is set to $N = 15$.

The optimization is subject to constraints of the form

$$\mathbf{c}^a(\mathbf{x}_k^a, \mathbf{u}_k^a, k\delta t^a) > \mathbf{0} \quad \forall k \in [0, N] \quad (22)$$

$$\mathbf{c}^s(\mathbf{x}_k^s, \mathbf{u}_k^s, (N-1)\delta t^a + k\delta t^s) > \mathbf{0} \quad \forall k \in [0, N] \quad (23)$$

where $\mathbf{c}^a, \mathbf{c}^s$ are vectors of concatenated scalar constraints. The following constraints are applied: (i) the final state of the aerial phase $\mathbf{x}_k^s$ for the foot should be at zero height (guard transition constraint); and (ii) the ground should not be penetrated, i.e., $z_k^s > 0$, and the ground reaction force should be positive $F_{z,k}^s > 0$ during the stance phase.

By solving a standard optimization problem that minimizes the cost function (21) subject to the constraints (22)–(23), a baseline trajectory can be generated that satisfies all constraints by selecting appropriate inputs and time steps. However, as illustrated in Fig. 3, the baseline solution saturates the ground penetration constraint ($z^s \to 0$) and the stance ground reaction force constraint ($F_z \to 0$). Therefore, under parametric uncertainty small perturbations can likely induce undesired slip or premature liftoff transitions.

To address this limitation, we introduce a robust reformulation of the constraints using the computed sensitivities. From (7), the constraint tubes for each component of $\mathbf{c}^a$ and $\mathbf{c}^s$ can be evaluated to form radius vectors $\rho_{c,k}^a, \rho_{c,k}^s$. These quantities allow us to define robustified constraints:

$$\mathbf{c}^a(\mathbf{x}_k^a, \mathbf{u}_k^a, k\delta t^a) > \rho_{c,k}^a \quad (24)$$

$$\mathbf{c}^s(\mathbf{x}_k^s, \mathbf{u}_k^s, (N-1)\delta t^a + k\delta t^s) > \rho_{c,k}^s \quad (25)$$

As shown in Fig. 3, the robust trajectory avoids saturating both penetration and ground reaction force constraints, thereby providing a safety margin against uncertainty. Fig. 4 illustrates that this robustness is achieved through increased actuator effort, most notably by preloading the spring during the aerial phase and applying larger forces immediately after touchdown. This results in a large ground reaction force.

To evaluate performance, we simulated 100 open-loop trajectories under uniformly distributed parameter perturbations drawn from an ellipsoid with $\Delta\mathbf{p}_{\max} = (0.1, 10)$, corresponding to up to a 10% deviation in either parameter.

<footnote n="2">²Note that in this collocation framework, all of the state, input, timestep, and sensitivity variables are concatenated into a single vector that can be optimized by generic solvers like `Ipopt`. However, the high dimensionality of the $\mathbf{\Pi}$ matrices slow down the convergence of the solver. Future work could leverage constrained shooting methods [25] that would not require the sensitivity matrices to be included as optimization variables.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of the nominal baseline and robust optimization results for the vertical SLIP hopper. The baseline solution (purple) satisfies the ground penetration and ground reaction force constraints only by saturating them, yielding a trajectory that is not robust to parameter deviations. In contrast, the proposed sensitivity-based formulation (blue) produces a trajectory that avoids saturation, providing robustness against uncertainty.</figcaption></figure>

<table><caption>TABLE I

FAILURE RESULTS FOR 100 SIMULATED TRAJECTORIES WITH PERTURBED PARAMETERS</caption><thead><tr><th>Failure Mode</th><th>Baseline</th><th>Proposed</th></tr></thead><tbody><tr><td>Body Contact With Ground</td><td>50</td><td>7</td></tr><tr><td>Foot Liftoff</td><td>41</td><td>5</td></tr><tr><td>Total</td><td>91</td><td>12</td></tr></tbody></table>

The results, summarized in Tab. I, show that the baseline (non-robust) trajectory failed in 91% of trials, with violations roughly split between ground penetration and premature liftoff. In contrast, the proposed robust formulation reduced the failure rate to 12%. The remaining failures are attributable to nonlinear effects of the hybrid dynamics that exceed the first-order sensitivity approximation.

Figure 4 illustrates the evolution of body height and ground reaction force for both formulations. The nominal trajectories are shown in bold, the sensitivity-based tubes in shaded regions, and the perturbed trials as thin lines. The distribution of perturbed trajectories aligns closely with the predicted tubes. The baseline solution clearly leads to frequent constraint violations, while the robust approach maintains feasible trajectories within the safety margins.

## B. Planar Quadruped

We next consider a more complex legged robot model to study robust trajectory optimization in the context of stable landings. Specifically, we employ a planar quadruped model, which enforces sagittal symmetry where the front and rear pairs of legs move identically and the robot's body is restricted to planar motion with no roll or yaw rotation. The body is modeled with a mass *m* and moment of inertia *I*, while the legs are assumed to be massless, see Fig. 1 for a visual representation.

The system state is defined as $\mathbf{q} = (x, z, \theta, \alpha_f, \beta_f, \alpha_r, \beta_r)$ with corresponding velocities $\dot{\mathbf{q}} = (\dot{x}, \dot{z}, \dot{\theta}, \dot{\alpha}_f, \dot{\beta}_f, \dot{\alpha}_r, \dot{\beta}_r)$, and full state vector $\mathbf{x} = (\mathbf{q}, \dot{\mathbf{q}})$. *x* and *z* denote the body position, *θ* is the pitch angle, and $\alpha_f, \beta_f$ and $\alpha_r, \beta_r$ are the hip and knee angles of the front and rear legs, respectively.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Simulation results for the nominal inputs, body height and ground reaction force under parameter perturbations. Bold curves denote the nominal trajectories, shaded regions indicate sensitivity-based tubes, and thin lines represent 100 perturbed trials. Easily seen from the zoomed-in inset figures, the baseline solution (purple) frequently violates the constraints, whereas the robust formulation (blue) maintains feasible trajectories within the predicted tubes.</figcaption></figure>

Recent studies [14], [26], [27] have explored the use of parallel-spring elasticity in legged robots, highlighting its potential to improve energy efficiency during locomotion. However, the resonant dynamics introduced by such springs necessitate careful control to avoid instability. Motivated by these works, we augment the planar quadruped model by adding torsional spring-dampers to each hip and knee joint, with a nominal stiffness of *k* and damping coefficient of *b*. The rest positions of the hip and knee joints are at angles $-3\pi/4$ and $-\pi/4$ from the horizontal, respectively. The body mass and moment of inertia, along with the joint spring and damping coefficients are set as uncertain parameters, $\mathbf{p} = (m, I, k, b)$ with nominal values $\mathbf{p}_n = (5.0, 0.1, 15.0, 0.2)$.

The optimization structure follows that of (21), with the key difference that the planar quadruped now involves three phases and two hybrid transitions. The robot begins in the aerial domain with initial configuration $\mathbf{q}_0 = (0, 0.25, 0, -3\pi/4, -\pi/4, -3\pi/4, -\pi/4)$ corresponding to a body height of 0.25 m and the joints at their nominal positions. The initial velocity is $\dot{\mathbf{q}}_0 = (0.5, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0)$ with positive forward ve-
locity and negative rotational velocity. Due to the initial pitch rate, the front foot contacts the ground first, transitioning the system into a front-stance domain. When the rear foot subsequently touches down, the system enters the doublestance domain. The cost function penalizes input magnitudes in each domain and encourages the robot to reach a terminal state with all joints at their nominal positions and zero velocity. For this example, we set the number of discretization steps per domain to $N = 10$.

Controlling a planar quadruped in this three-phase setting introduces challenges beyond those encountered in simpler models like the SLIP hopper. In particular, we focus on two classes of failure modes: a vanishing ground reaction force and saturation of the friction cone. Similar to the SLIP example, the torsional spring dynamics can drive the ground reaction force to zero as the joints rebound after touchdown. Additionally, in the planar quadruped, friction cone constraints become critical and can be expressed as $(\mu F_z)^2 - F_x^2 > 0$, where $F_x$ and $F_z$ are the horizontal and vertical contact forces, and the coefficient of friction is set to $\mu = 0.5$. Violation of this constraint leads to slipping, which can destabilize the robot and result in a severe failure.

Following the setup of the previous example, a baseline trajectory is generated with ground reaction force and friction cone constraints enforced in the front-stance and doublestance domains. A sensitivity-aware trajectory is also generated using the robustified constraints. The maximum parameter variations are set to $\Delta \mathbf{p}_{\max} = (0.2, 0.02, 2.0, 0.02)$.

The baseline trajectory satisfies the ground reaction force and friction cone constraints in the nominal case but only by saturating them, as shown in Fig. 1. In contrast, the robust sensitivity-aware trajectory avoids constraint saturation. To evaluate performance under uncertainty, 625 perturbed trajectories were simulated for each method, and Tab. II shows that the baseline approach fails in 98.6% of trials, with failures concentrated on the front foot, where both ground reaction force and friction cone constraints are frequently violated. The sensitivity-aware trajectory, on the other hand, experiences failures in only 1.9% of trials, with 12 occurrences of slipping at the rear foot. These rare failures can again be attributed to nonlinearities in the hybrid dynamics.

Fig. 5 illustrates the evolution of the ground reaction force and friction cone constraints for the front and rear feet under both methods. For the baseline trajectory, the ground reaction force drops to zero at the end of the motion due to spring rebound, and the friction cone constraints are saturated at rear touchdown and near the end of the trajectory. The proposed sensitivity-aware method successfully avoids these saturations, maintaining feasible trajectories under perturbations.

The two examples demonstrate the inability of baseline methods to handle parametric uncertainties, with failure rates exceeding 90% in both experiments. Meanwhile, the proposed sensitivity-aware method achieves significantly better performance, even under substantial parameter perturbations. Importantly, the approach effectively models the evolution of uncertainties along the trajectory, enabling performance and robustness without being overly conservative.

<table><caption>TABLE II
FAILURE RESULTS FOR 625 SIMULATED TRAJECTORIES WITH PERTURBED PARAMETERS</caption><thead><tr><th>Failure Mode</th><th>Baseline</th><th>Proposed</th></tr></thead><tbody><tr><td>Front Liftoff</td><td>181 (29.0%)</td><td>0</td></tr><tr><td>Rear Liftoff</td><td>6 (0.1%)</td><td>0</td></tr><tr><td>Front Slip</td><td>333 (53.3%)</td><td>0</td></tr><tr><td>Rear Slip</td><td>96 (15.4%)</td><td>12 (1.9%)</td></tr><tr><td>Total</td><td>616 (98.6%)</td><td>12 (1.9%)</td></tr></tbody></table>

## VI. CONCLUSION & FUTURE WORK

This work presented a framework for computing parametric sensitivities in hybrid systems and demonstrated its application to robust trajectory optimization. The proposed methodology combines sensitivity analysis for both the continuous dynamics and the discrete events, including guard conditions and reset maps, and provides a unified view of the closed-loop parametric sensitivity of the system over the whole trajectory. The hybrid sensitivity has then been successfully applied to the problem of robust trajectory optimization for two canonical legged locomotion systems, in which we leveraged the hybrid sensitivity-based tubes to achieve constraint satisfaction in the presence of non-negligible uncertainties, showing how the proposed method is capable of outperforming vanilla trajectory optimization techniques in avoiding contact related failures. Future works will entail the validation on a real quadruped robot. Having highlighted the limitations of standard sensitivity analysis due to the fixed reset time, we also proposed an alternative sensitivity definition based on extended trajectories. This new definition provides a more accurate approximation of the hybrid system, but it also introduces a discontinuity that makes its application to the trajectory optimization framework less straightforward. For this reason, future works will focus on the computation of closed-form expressions for tubes of trajectories using this extended sensitivity and investigate its use in sampling-based trajectory optimization or data-efficient learning.

## REFERENCES

[1] J. van Steen, G. v. d. Brandt, N. v. d. Wouw, J. Kober, and A. Saccon, “Quadratic programming-based reference spreading control for dual-arm robotic manipulation with planned simultaneous impacts,” *IEEE Transactions on Robotics*, vol. 40, pp. 3341–3355, 2024.

[2] D. J. Lynch, K. M. Lynch, and P. B. Umbanhowar, “The soft-landing problem: Minimizing energy loss by a legged robot impacting yielding terrain,” *IEEE Robotics and Automation Letters*, vol. 5, no. 2, pp. 3658–3665, 2020.

[3] J. Ding, V. Atanassov, E. Panichi, J. Kober, and C. D. Santina, “Robust quadrupedal jumping with impact-aware landing: Exploiting parallel elasticity,” *IEEE Transactions on Robotics*, vol. 40, pp. 3212–3231, 2024.

[4] S. H. Jeon, S. Kim, and D. Kim, “Online optimal landing control of the mit mini cheetah,” in *2022 International Conference on Robotics and Automation (ICRA)*, 2022, pp. 178–184.

[5] A. Ansari and T. Murphey, “Minimum sensitivity control for planning with parametric and hybrid uncertainty,” *The International Journal of Robotics Research*, vol. 35, no. 7, p. 823–839, 2016.

[6] 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)*, 2018, pp. 286–293.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Simulation results for the ground reaction forces and friction cone constraints under parameter perturbations. Bold curves denote the nominal trajectories, shaded regions indicate sensitivity-based tubes, and thin lines represent 100 perturbed trials. For the baseline method, the front foot is shown in purple and the rear foot in red; for the proposed method, the front foot is in blue and the rear foot in orange. In the baseline case, the front and rear ground reaction forces drop to zero at the end of the trajectory, and the front-foot friction cone is saturated at rear touchdown. Additionally, the friction cone for both feet is saturated as the ground reaction forces vanish near the end of the motion. In contrast, the proposed method maintains robustness to parameter perturbations, avoiding constraint saturation.</figcaption></figure>

[7] 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, pp. 9962-9969, 2024.

[8] 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.

[9] J. Zhu, T. Simeon, and M. Cognetti, "Robust sensitivity-aware chance-constrained mpc for efficient handling of multiple uncertainty sources," *IEEE Robotics and Automation Letters*, vol. 10, no. 10, pp. 10330-10337, 2025.

[10] A. Pupa, T. Belvedere, C. Secchi, and P. Robuffo Giordano, "On the computation of sensitivity tubes," *IEEE Robotics and Automation Letters*, vol. 10, no. 9, pp. 8802-8809, 2025.

[11] A. Saccon, N. van de Wouw, and H. Nijmeijer, "Sensitivity analysis of hybrid systems with state jumps with application to trajectory tracking," in *53rd IEEE Conference on Decision and Control*, 2014, p. 3065-3070.

[12] N. J. Kong, J. Joe Payne, J. Zhu, and A. M. Johnson, "Saltation matrices: The essential tool for linearizing hybrid dynamical systems," *Proceedings of the IEEE*, vol. 112, no. 6, pp. 585-608, 2024.

[13] M. Rijnen, A. Saccon, and H. Nijmeijer, "On optimal trajectory tracking for mechanical systems with unilateral constraints," in *2015 54th IEEE Conference on Decision and Control*, 2015, p. 2561-2566.

[14] J. Zhu, J. J. Payne, and A. M. Johnson, "Convergent ilqr for safe trajectory planning and control of legged robots," in *2024 IEEE International Conference on Robotics and Automation (ICRA)*, 2024, p. 8051-8057.

[15] J. J. Payne, N. J. Kong, and A. M. Johnson, "The uncertainty aware salted kalman filter: State estimation for hybrid systems with uncertain guards," in *2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*, 2022, p. 8821-8828.

[16] M. Tucker, N. Csomay-Shanklin, and A. D. Ames, "Robust bipedal locomotion: Leveraging saltation matrices for gait optimization," in *2023 IEEE International Conference on Robotics and Automation (ICRA)*, 2023, p. 12218-12225.

[17] 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 *Int. Conf. on Unmanned Aircraft Systems (ICUAS)*. IEEE, 2024, p. 572-578.

[18] R. Goebel, R. G. Sanfelice, and A. R. Teel, "Hybrid dynamical systems," *IEEE Control Systems Magazine*, vol. 29, no. 2, p. 28-93, 2009.

[19] S. Geng and I. A. Hiskens, "Second-order trajectory sensitivity analysis of hybrid systems," *IEEE Transactions on Circuits and Systems I: Regular Papers*, vol. 66, no. 5, p. 1922-1934, 2019.

[20] D. Krishnamoorthy, "A sensitivity-based data augmentation framework for model predictive control policy approximation," *IEEE Transactions on Automatic Control*, vol. 67, no. 11, p. 6090-6097, 2022.

[21] A. Tagliabue and J. P. How, "Efficient deep learning of robust policies from mpc using imitation and tube-guided data augmentation," *IEEE Transactions on Robotics*, vol. 40, p. 4301-4321, 2024.

[22] I. Hiskens and M. Pai, "Trajectory sensitivity analysis of hybrid systems," *IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications*, vol. 47, no. 2, p. 204-220, 2000.

[23] T. A. McMahon and G. C. Cheng, "The mechanics of running: How does stiffness couple with speed?" *Journal of Biomechanics*, vol. 23, pp. 65-78, 1990.

[24] M. Kelly, "An introduction to trajectory optimization: How to do your own direct collocation," *SIAM Review*, vol. 59, no. 4, pp. 849-904, 2017.

[25] T. A. Howell, B. E. Jackson, and Z. Manchester, "Altro: A fast solver for constrained trajectory optimization," in *2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*, 2019, pp. 7674-7679.

[26] J. Ding, V. Atanassov, E. Panichi, J. Kober, and C. D. Santina, "Robust quadrupedal jumping with impact-aware landing: Exploiting parallel elasticity," *IEEE Transactions on Robotics*, vol. 40, pp. 3212-3231, 2024.

[27] M. J. Pollayil, C. Della Santina, G. Mesesan, J. Englberger, D. Seidel, M. Garabini, C. Ott, A. Bicchi, and A. Albu-Schaffer, "Planning natural locomotion for articulated soft quadrupeds," in *2022 International Conference on Robotics and Automation (ICRA)*, 2022, pp. 6593-6599.