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

HAL
open science

# Analytical derivatives of strain-based dynamic model for hybrid soft-rigid robots

Anup Teejo Mathew, Frédéric Boyer, Vincent Lebastard, Federico Renda

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

## To cite this version:

Anup Teejo Mathew, Frédéric Boyer, Vincent Lebastard, Federico Renda. Analytical derivatives of strain-based dynamic model for hybrid soft-rigid robots. The International Journal of Robotics Research, 2025, 10.1177/02783649251346209. hal-05230423

HAL Id: hal-05230423

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

Submitted on 29 Aug 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analytical Derivatives of Strain-Based Dynamic Model for Hybrid Soft-Rigid Robots

The International Journal of Robotics Research
©The Author(s) 2024
sagepub.co.uk/journalsPermissions.nav
DOI:
www.sagepub.com/

SAGE

Anup Teejo Mathew<fnref n="1" />,<fnref n="2" />, Frederic Boyer<fnref n="3" />, Vincent Lebastard<fnref n="3" />, Federico Renda<fnref n="1" />,<fnref n="2" />

## Abstract

Algorithms that use derivatives of governing equations have accelerated rigid robot simulations and improved their accuracy, enabling the modeling of complex, real-world capabilities. However, extending these methods to soft and hybrid soft-rigid robots is significantly more challenging due to the complexities in modeling continuous deformations inherent in soft bodies. A considerable number of soft robots and the deformable links of hybrid robots can be effectively modeled as slender rods. The Geometric Variable Strain (GVS) model, which employs the screw theory and the strain parameterization of the Cosserat rod, extends the rod theory to model hybrid soft-rigid robots within the same mathematical framework. Using the Recursive Newton-Euler Algorithm, we developed the analytical derivatives of the governing equations of the GVS model. These derivatives facilitate the implicit integration of dynamics and provide the analytical Jacobian of the statics residue, ensuring fast and accurate computations. We applied these derivatives to the mechanical simulations of six common robotic systems: a soft cable-driven manipulator, a hybrid serial robot, a fin-ray finger, a hybrid parallel robot, a contact scenario, and an underwater hybrid mobile robot. Simulation results demonstrate substantial improvements in computational efficiency, with speed-ups of up to three orders of magnitude (1000x). We validate the model by comparing simulations done with and without analytical derivatives. Beyond static and dynamic simulations, the techniques discussed in this paper hold the potential to revolutionize the analysis, control, and optimization of hybrid robotic systems for real-world applications.

## Keywords

Differentiable Simulation, Soft Robotics, Cosserat Rod, Mathematical Modeling, Screw Theory

## 1 Introduction

Rigid-body algorithms have been developed for the mechanical analysis of multi-body systems, enabling the modeling of dynamic response and static equilibrium of robots Featherstone (2008); Murray et al. (1994). Over the years, these algorithms have undergone remarkable improvements, enabling faster-than-real-time simulations Newbury et al. (2024). One of the most significant advancements has been the incorporation of gradient information into these algorithms Carpentier et al. (2019); Howell et al. (2023); Giftthaler et al. (2017a); Geilinger et al. (2020). The ability to accurately and efficiently compute the derivatives of governing equations with respect to the system's state, model parameters, and control variables has been pivotal for implicit integration, design optimization, trajectory optimization, and optimal control. The impact of these advancements is exemplified by leading humanoid and quadrupedal robots, such as those developed by Boston Dynamics and Unitree Bjelonic et al. (2022). These robots leverage the gradient information for trajectory optimization and model-predictive control (MPC), enabling real-time control, enhanced stability, and adaptability in complex environments Sukhija et al. (2023); Neunert et al. (2018); Wensing et al. (2024).

Several methods exist for computing the derivatives of governing equations, each with its trade-offs Newbury et al. (2024). The most straightforward approach is the numerical scheme of finite difference. The finite difference

method can offer simplicity and ease of parallelization but often falls short in accuracy Todorov et al. (2012). Automatic differentiation (AutoDiff) computes derivatives by recognizing that even complex functions are built from fundamental operations and functions. It systematically applies the chain rule to these operations within the algorithm, allowing the program to automatically calculate derivatives alongside the original calculations Tedrake and the Drake Development Team (2019); Giftthaler et al. (2017b). However, AutoDiff involves intermediate computations that are generally hard to simplify. Analytical methods take a manual approach by directly applying the chain rule to recursive algorithms like the Recursive Newton-Euler Algorithm (RNEA) Carpentier and Mansard (2018); Singh et al. (2022). By exploiting the inherent structure and spatial algebra at the core of rigid-body dynamics algorithms, analytical derivatives (AD) can

<footnote n="1">¹Department of Mechanical Engineering, Khalifa University, Abu Dhabi, UAE</footnote>

<footnote n="2">²Khalifa University Center for Autonomous Robotic Systems (KUCARS), Khalifa University, Abu Dhabi, UAE</footnote>

<footnote n="3">³LS2N Laboratory, Institut Mines Telecom Atlantique, Nantes 44307, France</footnote>

<footnote>**Corresponding author:**
Anup Teejo Mathew, Department of Mechanical Engineering, Khalifa University, Abu Dhabi, UAE.</footnote>

<footnote>Email: anup.mathew@ku.ac.ae, anupteejo@gmail.com</footnote>
simplify computations and achieve greater computational efficiency than automatic differentiation methods. Efficient implementation of analytical derivatives can lead to faster and more resource-efficient computations with improved accuracy and provide deeper insight into the mathematical structure of the derivatives. However, deriving analytical derivatives manually can be complex and time-consuming, making them challenging to implement Singh et al. (2022).

Deriving analytical derivatives in soft robots is significantly more challenging than in rigid body systems. The primary difficulty stems from the complex nature of deformable bodies, which undergo large, continuous deformations, making it highly challenging to derive closed-form equations for their dynamics. The general class of soft robots, which cannot be modeled as a system rods, require numerical methods based on 3D continuum mechanics such as Finite Element Methods (FEM) Duriez (2013) or Material Point Method (MPM) Hu et al. (2018). Analytical derivatives of FEM in robotics have been explored in several works, including Hoshyari et al. (2019); Hafner et al. (2019); Geilinger et al. (2020); Bächer et al. (2021); Jatavallabhula et al. (2021); Du et al. (2022). On the other hand, MPM is often referred to as naturally differentiable due to its particle-grid representation and the smooth interpolation between particles and the grid, which makes gradient computation more efficient Spielberg et al. (2023); Huang et al. (2021). FEM and MPM use maximal coordinate representations for rigid bodies, increasing the degrees of freedom (DoF) and the computational cost for hybrid soft-rigid robots.

A large portion of soft robots and the compliant links of hybrid soft-rigid robots can be effectively modeled as slender, rod-like structures, making them well-suited for analysis using the Cosserat rod theory Armanini et al. (2023). The Cosserat rod is a 1D continuum mechanics object that can model deformations of slender bodies, including twisting, bending in two axes, stretching, and shearing in two axes Cao and Tucker (2008). Four distinct research communities have leveraged Cosserat rod theory to address their specific challenges, each producing specialized numerical methods tailored to their needs. Ranking them by their order of appearance, these communities are: the geometrically exact FEM community, the ocean engineering community, the computer graphics community, and the robotics community. The geometrically exact FEM community has focused primarily on developing FEM software that can predict the movements and stresses of mechanisms undergoing large deformations Simo and Vu-Quoc (1988); Meier et al. (2017); Eugster and Harsch (2023). Meanwhile, the ocean engineering community has applied Cosserat rod models to the simulation of towed submarine cables, addressing the unique challenges posed by underwater environments Burgess (1992); Tjavaras et al. (1998). On the other hand, the computer graphics community has prioritized computational speed for interactive simulations of filament-like structures, such as hair, using the Discrete Elastic Rod (DER) approaches Bergou et al. (2008); Gazzola et al. (2018). Finally, the robotics community has concentrated on the simulation and control of soft or continuum robots to safely interact with their surroundings Rucker and Webster (2011); Till et al. (2019); Boyer et al. (2022).

To cater specifically to the needs of robotics, a novel parameterization of the configuration space of Cosserat rods has been proposed, focusing on strain fields rather than the traditional pose-based approach used in FEM and DER. This strain-based approach, referred to as the Geometric Variable Strain (GVS) method Renda et al. (2020); Boyer et al. (2020), is geometrically exact and aligns well with the demands of soft robotic applications. The GVS approach provides a reduced-order formulation with the least number of generalized coordinates for modeling their dynamic and compliant behavior. Among the various models based on the Cosserat rod, the GVS model stands out for its ability to merge the screw theory-based formulation of rigid robots with the strain parameterization of the rod. Its Lagrangian mechanics formulation with minimal generalized coordinates enables efficient analysis of hybrid soft-rigid robotic systems within a unified mathematical framework. The implementation of the GVS model, based on the approximate Magnus expansion of the rod's strain field, makes the soft body computationally equivalent to multidimensional, discrete rigid joints Mathew et al. (2024). The model has been extensively validated through comparisons with analytical models, FEM, and other rod models Boyer et al. (2023); Mathew et al. (2022a), as well as experimental studies Alkayas et al. (2025); Feliu-Talegon et al. (2024). In Mathew et al. (2022a, 2024), we implemented the GVS model for hybrid soft robots using a built-in implicit integration scheme in MATLAB called ode15s. When analytical derivatives (Jacobian) are not supplied, ode15s internally compute the Jacobian of the governing equations using a finite difference scheme. However, this can lead to longer computational times or cause the solver to stall due to errors introduced by the numerical approximation, particularly in high-DoF systems with constraints. The objective of this work is to develop and implement an analytical derivative for the GVS model, aiming to improve computational efficiency and robustness in the simulation of hybrid soft-rigid robots.

**Related works:** The Piecewise Constant Strain (PCS) model, which discretizes the continuous deformation of a Cosserat rod into a finite number of segments with constant strain, is a subclass of the GVS model. Based on the Comprehensive Motion Transformation Matrix (a Lie group of coordinate transformations of displacement, velocity, and acceleration), an analytical gradient of the PCS model was derived in Ishigaki et al. (2024). A differentiable soft robot simulation environment called DisMech was introduced based on the implicit DER method Choi et al. (2024). DisMech employs a finite difference scheme to compute the necessary gradients for the implicit integration. Recently, a new algorithm for the implicit dynamics of robots with rigid bodies and Cosserat rods has been proposed Boyer et al. (2023). By applying an exact symbolic differentiation of the robot's RNEA inverse dynamics (named IDM in Boyer et al. (2023)), a new RNEAAlgorithm, called the tangent inverse dynamics model (TIDM), has been derived. This algorithm is then fed with unit inputs to numerically calculate the Jacobian of the inverse dynamics. This calculation is performed at the continuous level, i.e. directly on the partial differential equations (PDEs), and before spatial integration, for which spectral methods have been employed instead
of Magnus expansion. Although the approach adheres to the true definition of geometrically exact methods - where discretization of time and space occurs only after all analytical calculations - due to its implicit nature, it does not directly provide Jacobian matrices in analytical matrix form, which can be advantageous for control and fast simulation of robots with MATLAB. In contrast, the approach here presented lies in its analytic and explicit formulation, which was made possible by the Magnus expansion.

**Contributions of this work:** By leveraging the “pseudo-rigid joint” formulation of GVS and building upon established methods for analytical derivatives in rigid body systems Carpentier and Mansard (2018); Singh et al. (2022), we developed analytical derivatives for soft and hybrid soft-rigid robots with slender soft bodies. While Mathew et al. (2024) formulated the governing equations of the GVS model, this paper focuses on deriving their analytical gradients and developing an efficient algorithm for their computation. We implemented the derivatives in two implicit integration algorithms: ode15s of MATLAB and Newmark-β scheme for dynamics and provided the analytical Jacobian for efficient static equilibrium computation. Our method significantly improved the computational speed of implicit integration for dynamic simulations, achieving speed-ups of up to three orders of magnitude compared to traditional methods without analytical derivatives. Similarly, for static analysis, we observed substantial improvements in computational efficiency. Six common robotic systems are considered for the simulation study. Additionally, we implemented the theory developed in this paper into our GUI-based MATLAB toolbox, SoRoSim, making it a fully differentiable robotics toolbox capable of modeling generic hybrid soft-rigid robots with slender soft rods. Mathew et al. (2025). To the best of the authors’ knowledge, this is the first work to derive and implement analytical derivatives for the mechanical analysis of hybrid soft-rigid robotic systems of this kind.

**Summary of Contributions:**

* Derivation of the analytical gradient of the GVS model.

* Integration of the gradient information for static and dynamic simulations.

* Analysis of different robotic systems validating the analytical gradient and demonstrating computational efficiency.

* Integration of the theory into the SoRoSim toolbox, making it a differentiable tool for modeling hybrid soft-rigid robots.

**Organization of the paper:** A summary of the GVS model is presented in Section 2. Section 3 details the implementation of analytical derivatives in dynamic and static algorithms for fast and accurate computations. In Section 4, we focus initially on a single soft body, applying the RNEA and *ID* framework to derive the analytical derivatives of the governing equations. This framework is extended to hybrid multi-body systems in Section 5. We address two typical constraints in multi-body systems: joint actuation via joint coordinates and closed-chain (CC) systems. Section 6 discusses the derivation of analytical derivatives for systems subjected to three common external

forces: point forces, contact loads, and hydrodynamic forces. Simulations and validations across six robotic systems are presented across these sections, demonstrating the effectiveness of the approach.

Table 1 lists all the important symbols used in this paper. Readers are encouraged to refer to the supplementary videos to visualize the dynamic simulation results presented in this work.

<table><caption>**Table 1.** List of symbols and their descriptions</caption><thead><tr><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>n<sub>dof</sub></td><td>Total degrees of freedom</td></tr><tr><td>n<sub>a</sub></td><td>Number of actuators</td></tr><tr><td>N</td><td>Number of Cosserat rods (including rigid joints)</td></tr><tr><td>n<sub>p</sub></td><td>Number of computational points</td></tr><tr><td>q</td><td>Generalized coordinates</td></tr><tr><td>q̇</td><td>Generalized velocities</td></tr><tr><td>q̇̇</td><td>Generalized accelerations</td></tr><tr><td>x</td><td>Generalized robot state. x = [q<sup>T</sup> q̇<sup>T</sup>]<sup>T</sup></td></tr><tr><td>h</td><td>Time step</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>X</td><td>Curvilinear abscissa of the soft body, X ∈ [0, L]</td></tr><tr><td>FD</td><td>Forward dynamics</td></tr><tr><td>ID</td><td>Inverse dynamics</td></tr><tr><td>M</td><td>Generalized mass matrix</td></tr><tr><td>F</td><td>Generalized external and Coriolis force</td></tr><tr><td>τ</td><td>Generalized internal force</td></tr><tr><td>B</td><td>Generalized actuation matrix</td></tr><tr><td>u</td><td>Actuator strength</td></tr><tr><td>K</td><td>Generalized stiffness matrix</td></tr><tr><td>D</td><td>Generalized damping matrix</td></tr><tr><td>ξ</td><td>Screw strain</td></tr><tr><td>Φ<sub>ξ</sub></td><td>Strain basis</td></tr><tr><td>(●)<sub>α</sub></td><td>Quantity at joint or computational point α</td></tr><tr><td>(●)<sup>B</sup></td><td>Quantities computed during the forward pass</td></tr><tr><td>(●)<sup>C,S</sup></td><td>Quantities computed during the backward pass</td></tr><tr><td>S<sub>α</sub></td><td>Joint motion subspace matrix</td></tr><tr><td>F<sub>k</sub></td><td>Local wrench</td></tr><tr><td>M<sub>k</sub></td><td>Screw inertia matrix</td></tr><tr><td>η<sub>k</sub></td><td>Velocity twist</td></tr><tr><td>η̇<sub>k</sub></td><td>Acceleration twist</td></tr><tr><td>G</td><td>G = [0<sup>T</sup>a<sub>g</sub><sup>T</sup>]<sup>T</sup>, where a<sub>g</sub> is the acceleration due to gravity</td></tr><tr><td>g</td><td>Homogeneous transformation matrix in SE(3)</td></tr><tr><td>J</td><td>Geometric (kinematic) Jacobian matrix</td></tr><tr><td>Ad(.)</td><td>Adjoint maps in SE(3)</td></tr><tr><td>Ad<sup>*</sup>(.)</td><td>adjoint operators in se(3)</td></tr><tr><td>ad(.)</td><td></td></tr><tr><td>ad<sup>*</sup>(.), ad̄(.)</td><td></td></tr></tbody></table>

## 2 Summary of the GVS Model

The GVS model introduces generalized coordinates (q) using a variable strain parameterization of the Cosserat rod:
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Graphical summary of the GVS model: (a) Schematic of the implemented recursive scheme. (b) Block diagram showing the summary of the GVS FD algorithm.</figcaption></figure>

$$ \xi_i(X, q_i) = \Phi_{\xi_i}(X)q_i + \xi_i^*(X) \quad (1) $$

where $X \in [0, L_i]$ is the curvilinear abscissa of the rod, $\xi \in \mathbb{R}^6$ is the screw strain, $\Phi_{\xi_i} \in \mathbb{R}^{6 \times n_{dof}}$ ($n_{dof}$ is the degrees of freedom) is a strain basis, and $\xi^*$ is the reference strain. The subscript $i$ indicates the $i^{th}$ Cosserat rod. In the GVS formulation, a rigid joint is equivalent to a fictitious Cosserat rod spanning from $X=0$ to $X=1$ Mathew et al. (2024). For a rigid joint, the strain twist is equivalent to the joint twist in $\mathbb{R}^6$ and is independent of $X$.

For a hybrid robot with $N$ Cosserat rods (including rigid joints and soft bodies), the state of the robot $x$ is governed by the generalized coordinates $q \in \mathbb{R}^{n_{dof}}$ ($n_{dof} = \sum_{i=1}^N n_{dof_i}$) and its time derivative $\dot{q}$. Using $q, \dot{q}$ and basis functions, $\xi$ and $\dot{\xi}$ can be computed at each rigid joint and at any computational point on the soft bodies.

Using these, a recursive scheme that employs the exponential map of the Lie algebra of $SE(3)$ is implemented to compute the robot's forward and differential kinematics (Appendix A). This yields $g_i \in SE(3)$, the homogeneous transformation matrix mapping from the inertial frame to the local frame, $J_i \in \mathbb{R}^{6 \times n_{dof}}$, the geometric Jacobian in the local frame, and its time derivative $\dot{J}_i$ at discrete computational point along the rod domains.

$$ g_{i,\alpha+1} = g_{i,\alpha} \exp(\hat{\Omega}_{i,\alpha}) \quad (2a) $$

$$ J_{i,\alpha+1} = \mathrm{Ad}^{-1}_{\exp(\hat{\Omega}_{i,\alpha})} \left( J_{i,\alpha} + \left[ \mathbf{0}_{6 \times n_i^{-}} \mathbf{S}_{\alpha} \mathbf{0}_{6 \times n_i^{+}} \right] \right) \quad (2b) $$

$$ \dot{J}_{i,\alpha+1} = \mathrm{Ad}^{-1}_{\exp(\hat{\Omega}_{i,\alpha})} \left( \dot{J}_{i,\alpha} + \left[ \mathbf{0}_{6 \times n_i^{-}} \dot{\mathbf{S}}_{\alpha} + \mathrm{ad}_{\eta_{i,\alpha}} \mathbf{S}_{\alpha} \mathbf{0}_{6 \times n_i^{+}} \right] \right) \quad (2c) $$

where, $\hat{\Omega}_{i,\alpha}(q_i)$ is an approximation of the Magnus expansion of $\xi_i$, representing an equivalent joint twist from $\alpha$ to $\alpha+1$, $\mathbf{S}_{\alpha}(q_i) \in \mathbb{R}^{6 \times n_{dof}}$ is the joint motion subspace matrix, and $\dot{\mathbf{S}}_{\alpha}(q_i, \dot{q}_i)$ is its time derivative. $\eta_{i,\alpha} = J_{i,\alpha} \dot{q}$ is the screw velocity in the local frame. The matrix dimensions $n_i^{-} = \sum_{k=0}^{i-1} n_{dof_k}$ and $n_i^{+} = \sum_{k=i+1}^{N} n_{dof_k}$ ($\hat{\cdot}$) is an isomorphism from $\mathbb{R}^6$ to $\mathfrak{SE}(3)$ and $\mathrm{Ad}(\cdot)$ represents

the Adjoint map in $SE(3)$. Note that these equations also apply to rigid joints, providing a mapping from $X=0$ to $X=1$. Readers may refer to Appendices A and C for their analytical expressions.

Based on D'Alembert's principle, the free dynamics of soft links and rigid bodies are projected onto the generalized coordinate space using the geometric Jacobian (Kane method Kane and Levinson (1983)), yielding the generalized mass matrix $M$, the generalized internal force $\tau$, and the generalized external and Coriolis forces $F$. The soft body is typically discretized into computational points for numerical spatial integration. For this numerical integration, the continuous domain of $X$ is discretized into $n_p$ computational points, and the integral of a distributed quantity $\bar{Q}$ is approximated as a weighted sum of the integrands.

$$ \int_{0}^{L_i} \bar{Q}_i dX = \sum_{k=1}^{n_p} w_k \bar{Q}_{ik} \quad (3) $$

Finally, the Forward Dynamics (FD) computes $\ddot{q}$ by solving,

$$ M(q)\ddot{q} = \tau(q, \dot{q}, u) + F(q, \dot{q}) \quad (4) $$

where, $u \in \mathbb{R}^{n_a}$ is the vector of applied actuation forces.

Figure 1 presents a graphical summary of the model implementation. The dynamic simulation computes the temporal evolution of the robot's state by time-integrating $[\dot{q}^T, \ddot{q}^T]^T$, while static simulation involves the computation of $q$ for which the following equation is satisfied.

$$ \tau(q, 0, u) + F(q, 0) = 0 \quad (5) $$

For further details on the model and its computational implementation, readers are encouraged to refer to Mathew et al. (2024).

## 3 Derivatives for Dynamics and Statics

Analytical derivatives of governing equations, (4) and (5), are crucial for efficiently solving dynamic and static
problems. Here, we explain how the derivatives can be utilized to ensure accurate and fast computations.

## 3.1 Implicit Euler Integration Using Inverse Dynamics

The temporal evolution of the robot state $\mathbf{x} = [\mathbf{q}^T, \dot{\mathbf{q}}^T]^T$ can be computed numerically by integrating $\dot{\mathbf{x}}$ using an explicit or implicit scheme *Butcher* (2016).

The explicit Euler method is given by:

$$ \mathbf{x}_{n+1} = \mathbf{x}_n + h\dot{\mathbf{x}}(t_n, \mathbf{x}_n) \quad (6) $$

where $t$ is the time, $h$ is the step size, and the subscripts $n$ and $n+1$ indicate current and next time steps respectively.

For an explicit integrator, $\dot{\mathbf{x}}$ is a function of the current robot state. In contrast, the implicit Euler method is:

$$ \mathbf{x}_{n+1} = \mathbf{x}_n + h\dot{\mathbf{x}}(t_{n+1}, \mathbf{x}_{n+1}) \quad (7) $$

Since $\mathbf{x}_{n+1}$ appears on both sides of the equation, it necessitates a solution through iteration. The equation can be solved for the unknown $\mathbf{x}_{n+1}$ by making it a residual:

$$ \mathbf{R}_{dyn}(\mathbf{x}_{n+1}) = \mathbf{x}_{n+1} - \mathbf{x}_n - h\dot{\mathbf{x}}(t_{n+1}, \mathbf{x}_{n+1}) \quad (8) $$

The Jacobian of this equation is given by:

$$ \mathbf{J}_{dyn}(\mathbf{x}_{n+1}) = \mathbf{I} - h\left.\frac{\partial\dot{\mathbf{x}}}{\partial\mathbf{x}}\right|_{n+1} \quad (9) $$

With an initial guess of $\mathbf{x}_{n+1_{guess}} = \mathbf{x}_n$, an iterative solver such as the Newton-Raphson method uses the residual and the Jacobian to march towards the solution. We used the ode15s function in MATLAB, a variable-step, variable-order ODE solver that uses an advanced form of the implicit Euler method. Similar to the Euler scheme, ode15s utilizes the state derivative ($\dot{\mathbf{x}}$) and its Jacobian for time integration. For a Lagrangian system governed by (4),

$$ \dot{\mathbf{x}} = \begin{bmatrix} \dot{\mathbf{q}} \\ FD \end{bmatrix} \quad (10a) $$

$$ \frac{\partial \dot{\mathbf{x}}}{\partial \mathbf{x}} = \begin{bmatrix} 0_{n_{dof} \times n_{dof}} & \mathbf{I}_{n_{dof}} \\ \frac{\partial FD}{\partial \mathbf{q}} & \frac{\partial FD}{\partial \dot{\mathbf{q}}} \end{bmatrix} \quad (10b) $$

While the Jacobian can be computed numerically or analytically, the analytical computation can provide faster and more accurate results. To obtain the derivatives of FD for the computation of the Jacobian, we take a partial derivative of (4) with respect to $\mathbf{q}$:

$$ \frac{\partial \mathbf{M}}{\partial \mathbf{q}} \ddot{\mathbf{q}} + \mathbf{M} \frac{\partial FD}{\partial \mathbf{q}} = \frac{\partial \tau}{\partial \mathbf{q}} + \frac{\partial F}{\partial \mathbf{q}} $$

By rearranging terms, we get:

$$ \frac{\partial FD}{\partial \mathbf{q}} = \mathbf{M}^{-1} \left( \frac{\partial \tau}{\partial \mathbf{q}} - \frac{\partial ID}{\partial \mathbf{q}} \right) \quad (11) $$

where $ID$ is the inverse dynamics given by,

$$ ID(\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}}) = \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} - \mathbf{F}(\mathbf{q}, \dot{\mathbf{q}}) \quad (12) $$

with $\ddot{\mathbf{q}}$ treated as an independent variable.

Similarly, we get:

$$ \frac{\partial FD}{\partial \dot{\mathbf{q}}} = \mathbf{M}^{-1} \left( \frac{\partial \tau}{\partial \dot{\mathbf{q}}} - \frac{\partial ID}{\partial \dot{\mathbf{q}}} \right) \quad (13) $$

## 3.2 Newmark-β Method

The Newmark $\beta$ method is another advanced implicit scheme, widely used in structural mechanics where it has been improved over the years thanks to the HHT or the Generalized-$\alpha$-method among others Cardona and Géradin (1994); Chung and Hulbert (1993). Recently, it has been shown that the Newmark scheme is symplectic and variational Kane et al. (2000), which may explain the reasons for its success in the FEM community. In the Newmark scheme, $\ddot{\mathbf{q}}_{n+1}, \dot{\mathbf{q}}_{n+1}$, and $\mathbf{q}_{n+1}$ are solved using equation (4) together with:

$$ \begin{aligned} \dot{\mathbf{q}}_{n+1} &= \frac{\gamma}{\beta h} (\mathbf{q}_{n+1} - \mathbf{q}_n) + \left(1 - \frac{\gamma}{\beta}\right) \dot{\mathbf{q}}_n + h \left(1 - \frac{\gamma}{2\beta}\right) \ddot{\mathbf{q}}_n \\ \ddot{\mathbf{q}}_{n+1} &= \frac{1}{\beta h^2} (\mathbf{q}_{n+1} - \mathbf{q}_n) - \frac{1}{\beta h} \dot{\mathbf{q}}_n + \left(1 - \frac{1}{2\beta}\right) \ddot{\mathbf{q}}_n \end{aligned} \quad (14) $$

where $\beta$ and $\gamma$ are parameters that can be adjusted to control the stability and accuracy of the method. For instance, $(\beta, \gamma) = (1/4, 1/2)$ ensures second-order accuracy with no damping. Using equation (14), we can convert (4) into a residue of $\mathbf{q}_{n+1}$:

$$ \mathbf{R}_{dyn}(\mathbf{q}_{n+1}) = \tau(\mathbf{q}_{n+1}) - ID(\mathbf{q}_{n+1}) \quad (15) $$

where the Jacobian is given by:

$$ \mathbf{J}_{dyn}(\mathbf{q}_{n+1}) = \frac{\partial \tau}{\partial \mathbf{q}} + \frac{\gamma}{\beta h} \frac{\partial \tau}{\partial \dot{\mathbf{q}}} - \left( \frac{\partial ID}{\partial \mathbf{q}} + \frac{\gamma}{\beta h} \frac{\partial ID}{\partial \dot{\mathbf{q}}} + \frac{1}{\beta h^2} \frac{\partial ID}{\partial \ddot{\mathbf{q}}} \right) \quad (16) $$

The solution for $\mathbf{q}_{n+1}$ can be obtained iteratively, starting with an initial guess $\mathbf{q}_{n+1_{guess}} = \mathbf{q}_n$ and using a time step $h$. The process can be computationally efficient if the analytical Jacobian is provided. The Newmark-$\beta$ method can be efficient, particularly for index-3 Differential Algebraic Equations (DAEs) with a large number of kinematic constraints. We implemented from scratch a custom MATLAB code based on dynamic equations (15) and (16). Unlike ode15s, the Newmark scheme works with a fixed time step, which has been chosen for each example, so that the position offset between the custom code output and ode15s, is of the order of millimeters.

## 3.3 Computation of Static Equilibrium

The role of the Jacobian in the statics simulation is straightforward. The residual and Jacobian for solving the static equilibrium can be derived from (5) and (12) with $\dot{\mathbf{q}} = 0$ and $\ddot{\mathbf{q}} = 0$.

$$ \mathbf{R}_{sta}(\mathbf{q}) = \tau - ID \quad (17a) $$

$$ \mathbf{J}_{sta}(\mathbf{q}) = \frac{\partial \tau}{\partial \mathbf{q}} - \frac{\partial ID}{\partial \mathbf{q}} \quad (17b) $$

Numerical solvers such as `fsolve` in MATLAB, use (17a) with an initial guess $\mathbf{q}_{guess}$ and march towards a solution using the Jacobian (17b). If the analytical Jacobian is not provided, the solver uses numerical techniques, such as finite differences, to compute it. Therefore, with the analytical derivative of the residue, the problem can be solved quickly and efficiently.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Schematic of the RNEA for a single soft body showing the forward pass for kinematics and the backward pass for joint wrench computation. The joint motion subspace ($\mathcal{S}_{\alpha}$) projects the resultant wrench acting on the joint into the generalized coordinate space.</figcaption></figure>

# 4 Derivatives of Soft Body Dynamics

We begin by deriving the derivative for a single soft body and then extend this approach to a hybrid multi-body system. For simplicity, we omit the subscript $i$ in this section. By virtue of the Magnus expansion, the soft body is computationally equivalent to $n_p - 1$ rigid joints governed by the same $q$. The computation of the partial derivative of FD requires determining the partial derivative of ID, according to equations (11) and (13). The most efficient algorithm for the computation of ID is the RNEA Featherstone (2008). It is a two-pass algorithm with a forward pass that calculates link kinematics and a backward pass that determines the joint wrench needed to produce this motion. The schematic of RNEA for a soft body is shown in Figure 2. According to RNEA, the ID of the discretized soft body is given by:

$$ID = \sum_{\alpha=1}^{n_p-1} ID_{\alpha} = \sum_{\alpha=1}^{n_p-1} S_{\alpha}^{T} \mathcal{F}_{\alpha}^{C} \quad (18)$$

where $\mathcal{F}_{\alpha}^{C}$ is the resultant of all the dynamic (inertial and Coriolis) and external wrenches (such as gravity) acting on all points $k > \alpha$, transformed to the frame of $\alpha$ (Figure 2).

$$\mathcal{F}_{\alpha}^{C} = \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}}^{*} \mathcal{F}_{k} \quad (19)$$

where $\mathrm{Ad}_{(\cdot)}^{*} \in \mathbb{R}^{6 \times 6}$ is coadjoint map of $SE(3)$ and $\mathcal{F}_k$ is the resultant point wrench in the local frame of $k$, given by the sum of the inertial and Coriolis forces minus the external force. For a distributed wrench ($\overline{\mathcal{F}}_k$), an equivalent point wrench, given by $\mathcal{F}_k = w_k \overline{\mathcal{F}}_k$, where $w_k$ is the quadrature weight for integration, is considered. For convenience, we drop the overbar above all the distributed quantities using this rule. We have,

$$\mathcal{F}_k = \mathcal{M}_k \dot{\eta}_k + \mathrm{ad}_{\eta_k}^* \mathcal{M}_k \eta_k - \mathcal{M}_k \mathrm{Ad}_{g_k}^{-1} \mathcal{G} \quad (20)$$

where, $\mathcal{M}_k \in \mathbb{R}^{6 \times 6}$ is the screw inertia matrix, $\eta_k$ and $\dot{\eta}_k$ and the local velocity and acceleration twists, $\mathcal{G} = [\mathbf{0}^T \mathbf{a}_g^T]^T$, where $\mathbf{a}_g$ is the acceleration due to gravity in the global frame, and $\mathrm{ad}_{(\cdot)}^*$ is the coadjoint operator (Appendix A).

With this we proceed to compute the partial derivative of ID with respect to the states of the robot. A summary of the results is presented below. For detailed derivations, readers are encouraged to refer to Appendix C and D.

## 4.1 Partial Derivative of ID With Respect to Generalized Coordinates

Substituting (19) in (18), we can see that the partial derivative of $ID_{\alpha}$ with respect to $q$ is given by:

$$\frac{\partial ID_{\alpha}}{\partial q} = \frac{\partial S_{\alpha}^{T}}{\partial q} \mathcal{F}_{\alpha}^{C} + S_{\alpha}^{T} \sum_{k=\alpha+1}^{n_p} \frac{\partial \mathrm{Ad}_{g_{\alpha k}}^{*}}{\partial q} \mathcal{F}_{k} + S_{\alpha}^{T} \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}}^{*} \frac{\partial \mathcal{F}_{k}}{\partial q} \quad (21)$$

The analytical formula for the first term is provided in the Appendix C. The second and third terms involve sum over $k$ from $\alpha + 1$ to $n_p$. The derivations of their summands are given in the Appendix D. We get,

$$\frac{\partial \mathrm{Ad}_{g_{\alpha k}}^{*}}{\partial q} \mathcal{F}_{k} = \sum_{\beta=\alpha}^{k-1} \mathrm{Ad}_{g_{\alpha \beta}}^{*} \overline{\mathrm{Ad}}_{\mathrm{Ad}_{g_{\beta k}}^{*}}^{*} \mathcal{F}_{k} S_{\beta} \quad (22)$$

and

$$\frac{\partial \mathcal{F}_k}{\partial q} = \mathcal{N}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{R}_\beta + \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{Q}_\beta \quad (23)$$

where,

$$\mathcal{N}_k = \overline{\mathrm{Ad}}_{\mathcal{M}_k \eta_k}^* + \mathrm{ad}_{\eta_k}^* \mathcal{M}_k - \mathcal{M}_k \mathrm{ad}_{\eta_k} \quad (24a)$$

$$\mathbf{R}_{\beta} = \mathrm{ad}_{\eta_{\beta}^{+}} \mathbf{S}_{\beta} + \frac{\partial \mathbf{S}_{\beta}}{\partial q} \dot{\mathbf{q}} \quad (24b)$$

$$\mathbf{Q}_{\beta} = \mathrm{ad}_{\eta_{\beta}^{+}} \mathbf{S}_{\beta} + \mathrm{ad}_{\eta_{\beta}^{+}} \mathbf{R}_{\beta} + \mathrm{ad}_{\eta_{\beta}} \frac{\partial \mathbf{S}_{\beta}}{\partial q} \dot{\mathbf{q}} + \frac{\partial \dot{\mathbf{S}}_{\beta}}{\partial q} \dot{\mathbf{q}} + \frac{\partial \mathbf{S}_{\beta}}{\partial q} \ddot{\mathbf{q}} - \mathrm{ad}_{\mathrm{Ad}_{g_{\beta}^{-1}} \mathbf{G}} \mathbf{S}_{\beta} \quad (24c)$$

The definitions of the adjoint operators and the analytical formulas of $\frac{\partial \mathbf{S}_{\beta}}{\partial q} \dot{\mathbf{q}}$, $\frac{\partial \dot{\mathbf{S}}_{\beta}}{\partial q} \dot{\mathbf{q}}$, and $\frac{\partial \mathbf{S}_{\beta}}{\partial q} \ddot{\mathbf{q}}$ are provided in the Appendices A and C. $\eta_{\beta}^{+}$, $\dot{\eta}_{\beta}^{+}$ are the velocity and acceleration twists of $\beta + 1$ expressed in the frame of $\beta$. $\mathcal{N}_k$ has the same dimensions as $\mathcal{M}_k (\mathbb{R}^{6 \times 6})$, while $\mathbf{R}_{\beta}$ and $\mathbf{Q}_{\beta} \in \mathbb{R}^{6 \times n_{\mathrm{dof}}}$, similar to $\mathbf{S}_{\beta}$. Now we proceed to take the sum over $k$ from $\alpha + 1$ to $n_p$. Substituting (22) and (23) into (21) we get:

$$\frac{\partial ID_{\alpha}}{\partial q} = \frac{\partial S_{\alpha}^{T}}{\partial q} \mathcal{F}_{\alpha}^{C} + S_{\alpha}^{T} (\mathcal{N}_{\alpha}^{C} \mathbf{R}_{\alpha}^{B} + \mathcal{M}_{\alpha}^{C} \mathbf{Q}_{\alpha}^{B} + \mathbf{U}_{\alpha}^{S} + \mathbf{P}_{\alpha}^{S}) \quad (25)$$
where, $R_{\alpha}^{B}$, $Q_{\alpha}^{B}$, $\mathcal{F}_{\alpha}^{C}$, $\mathcal{N}_{\alpha}^{C}$, $\mathcal{M}_{\alpha}^{C}$, $U_{\alpha}^{S}$, and $P_{\alpha}^{S}$ are computed efficiently using recursive formulas according to:

$$R_{\alpha}^{B} = \mathrm{Ad}_{g_{\alpha-1\alpha}}^{-1}(R_{\alpha-1} + R_{\alpha-1}^{B}) \quad (26a)$$

$$Q_{\alpha}^{B} = \mathrm{Ad}_{g_{\alpha-1\alpha}}^{-1}(Q_{\alpha-1} + Q_{\alpha-1}^{B}) \quad (26b)$$

$$\mathcal{F}_{\alpha}^{C} = \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*}(\mathcal{F}_{\alpha+1} + \mathcal{F}_{\alpha+1}^{C}) \quad (26c)$$

$$\mathcal{N}_{\alpha}^{C} = \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*}(\mathcal{N}_{\alpha+1} + \mathcal{N}_{\alpha+1}^{C})\mathrm{Ad}_{g_{\alpha\alpha+1}}^{-1} \quad (26d)$$

$$\mathcal{M}_{\alpha}^{C} = \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*}(\mathcal{M}_{\alpha+1} + \mathcal{M}_{\alpha+1}^{C})\mathrm{Ad}_{g_{\alpha\alpha+1}}^{-1} \quad (26e)$$

$$U_{\alpha}^{S} = \mathcal{N}_{\alpha}^{C} R_{\alpha} + \mathcal{M}_{\alpha}^{C} Q_{\alpha} + \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*} U_{\alpha+1}^{S} \quad (26f)$$

$$P_{\alpha}^{S} = \overline{\mathrm{Ad}}_{\mathcal{F}_{\alpha}^{C}}^{*} S_{\alpha} + \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*} P_{\alpha+1}^{S} \quad (26g)$$

Finally, we obtain:

$$\frac{\partial ID}{\partial \dot{q}} = \sum_{\alpha=1}^{n_p-1} \frac{\partial ID_{\alpha}}{\partial \dot{q}} \quad (27)$$

Note that $R_{\alpha}^{B}$ is the partial derivative of $\eta_{\alpha}$: $\frac{\partial \eta_{\alpha}}{\partial \dot{q}} = R_{\alpha}^{B}$, while $Q_{\alpha}^{B}$ (excluding the gravity term) is associated with the partial derivative of $\dot{\eta}_{\alpha}$: $\frac{\partial \dot{\eta}_{\alpha}}{\partial \dot{q}} = Q_{\alpha}^{*B} - \mathrm{Ad}_{\eta_{\alpha}} R_{\alpha}^{B}$. Quantities with the superscript B are kinematic quantities determined during the forward pass, while those with superscripts C or S are evaluated during the backward pass. Similar to $\mathcal{F}_{\alpha}^{C}$, $\mathcal{N}_{\alpha}^{C}$ and $\mathcal{M}_{\alpha}^{C}$ are the resultant $\mathcal{N}_{k}$ and $\mathcal{M}_{k}$ acting on all points $k > \alpha$, transformed to the point of $\alpha$. All quantities with superscript B and S has a dimension of $\mathbb{R}^{6 \times n_{dof}}$. Hence, the $S_{\alpha}$ projection in (25) returns quantities with dimension $\mathbb{R}^{n_{dof} \times n_{dof}}$. The same rules apply for analytical derivatives with respect to $\dot{q}$ and $\ddot{q}$, discussed in the next sections.

## 4.2 Partial Derivative of ID With Respect to Generalized Velocities

The partial derivative of $ID_{\alpha}$ wrt $\dot{q}$ is given by:

$$\frac{\partial ID_{\alpha}}{\partial \dot{q}} = S_{\alpha}^{T} \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}}^{*} \frac{\partial \mathcal{F}_{k}}{\partial \dot{q}} \quad (28)$$

Only the Coriolis force is a function of $\dot{q}$. We derive (details in Appendix D),

$$\frac{\partial \mathcal{F}_k}{\partial \dot{q}} = \mathcal{N}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathcal{S}_{\beta} + \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathcal{Y}_{\beta} \quad (29)$$

where

$$\mathcal{Y}_{\beta} = \mathcal{R}_{\beta} + \mathrm{Ad}_{\eta_{\beta}} \mathcal{S}_{\beta} + \dot{\mathcal{S}}_{\beta} \quad (30)$$

Substituting (29) into (28) and taking the sum over $k$ from $\alpha + 1$ to $n_p$, we get:

$$\frac{\partial ID_{\alpha}}{\partial \dot{q}} = S_{\alpha}^{T} (\mathcal{N}_{\alpha}^{C} S_{\alpha}^{B} + \mathcal{M}_{\alpha}^{C} \mathcal{Y}_{\alpha}^{B} + V_{\alpha}^{S}) \quad (31)$$

where, $S_{\alpha}^{B}$, $Y_{\alpha}^{B}$, and $V_{\alpha}^{S}$ are recursively computed according to:

$$S_{\alpha}^{B} = \mathrm{Ad}_{g_{\alpha-1\alpha}}^{-1}(S_{\alpha-1} + S_{\alpha-1}^{B}) \quad (32a)$$

$$Y_{\alpha}^{B} = \mathrm{Ad}_{g_{\alpha-1\alpha}}^{-1}(Y_{\alpha-1} + Y_{\alpha-1}^{B}) \quad (32b)$$

$$V_{\alpha}^{S} = \mathcal{N}_{\alpha}^{C} S_{\alpha} + \mathcal{M}_{\alpha}^{C} Y_{\alpha} + \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*} V_{\alpha+1}^{S} \quad (32c)$$

Note that $S_{\alpha}^{B}$ is the partial derivative of $\eta_{\alpha}$: $\frac{\partial \eta_{\alpha}}{\partial \dot{q}} = S_{\alpha}^{B}$, while $Y_{\alpha}^{B}$ is associated with the partial derivative of $\dot{\eta}_{\alpha}$: $\frac{\partial \dot{\eta}_{\alpha}}{\partial \dot{q}} = Y_{\alpha}^{B} - \mathrm{Ad}_{\eta_{\alpha}} S_{\alpha}^{B}$. We have, $J_{\alpha} = S_{\alpha}^{B}$ and $\dot{J}_{\alpha} = Y_{\alpha}^{B} - R_{\alpha}^{B}$.

Summing (31) over $\alpha$ from 1 to $n_p - 1$, we have:

$$\frac{\partial ID}{\partial \dot{q}} = \sum_{\alpha=1}^{n_p-1} \frac{\partial ID_{\alpha}}{\partial \dot{q}} \quad (33)$$

## 4.3 Partial Derivative of ID With Respect to Generalized Accelerations

The partial derivative of $ID_{\alpha}$ with respect to $\ddot{q}$ is given by:

$$\frac{\partial ID_{\alpha}}{\partial \ddot{q}} = S_{\alpha}^{T} \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}}^{*} \frac{\partial \mathcal{F}_{k}}{\partial \ddot{q}} \quad (34)$$

Only the inertial force is a function of $\ddot{q}$. From (20), we get:

$$\frac{\partial \mathcal{F}_k}{\partial \ddot{q}} = \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathcal{S}_{\beta} \quad (35)$$

Substituting this into (34), we get:

$$\frac{\partial ID_{\alpha}}{\partial \ddot{q}} = S_{\alpha}^{T} (\mathcal{M}_{\alpha}^{C} S_{\alpha}^{B} + W_{\alpha}^{S}) \quad (36)$$

where, $W_{\alpha}^{S}$ is computed recursively according to:

$$W_{\alpha}^{S} = \mathcal{M}_{\alpha}^{C} S_{\alpha} + \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*} W_{\alpha+1}^{S} \quad (37)$$

Finally summing over all $\alpha$, we get:

$$\frac{\partial ID}{\partial \ddot{q}} = \sum_{\alpha=1}^{n_p-1} \frac{\partial ID_{\alpha}}{\partial \ddot{q}} = M(q) \quad (38)$$

Note that (38), incidentally, represents a novel way of computing the generalized mass matrix of a soft robot.

## 4.4 Partial Derivatives of Internal Forces

The internal forces include the elastic, damping, and actuation forces. It is given by:

$$\tau = Bu - Kq - D\dot{q} \quad (39)$$

where $D \in \mathbb{R}^{n_{dof} \times n_{dof}}$ is the generalized damping matrix, $K \in \mathbb{R}^{n_{dof} \times n_{dof}}$ is the generalized stiffness matrix, $B(q) \in \mathbb{R}^{n_{dof} \times n_a}$ ($n_a$ being the total number of actuators) is the generalized actuation matrix

The actuation matrix ($B$) for tendon-like actuators was derived in Renda et al. (2020); Renda et al. (2024), while a simple Hooke-like linear elastic and damping model is used for computing $K$ and $D$. The expressions for $B$, $K$, and $D$ can be found in the Appendix D.

The partial derivative of $\tau$ with respect to $q$ and $\dot{q}$ are given by:

$$\frac{\partial \tau}{\partial q} = \frac{\partial B}{\partial q} u - K \quad (40a)$$

$$\frac{\partial \tau}{\partial \dot{q}} = -D \quad (40b)$$

For tendon-like actuators, the partial derivative of $Bu$ with respect to $q$ is provided in the Appendix D.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Schematics of cable-driven soft manipulator with cable routing. $L = 50$ cm, base radius $r_b = 3$ cm, and tip radius $r_t = 1.5$ cm. RGB arrows indicate $x$, $y$, and $z$ axes of the global frame. (b) Actuation input used for the simulation.</figcaption></figure>

## 4.5 Cable-Driven Soft Manipulator

Based on the results so far, we set up the first simulation example. We simulate a cable-driven soft manipulator (CDM) shown in Figure 3(a). The manipulator radius varies linearly from base to tip according to $r(X) = r_b + (r_t - r_b)X/L$. The material properties used in the simulation are as follows: Young's modulus of 1 MPa, Poisson's ratio of 0.5, density of 1000 kg/m³, and a material damping coefficient of $10^4$ Pa·s. The soft manipulator has five tendon-like actuators, as shown in Figure 3(a). Their routing, local coordinates for a given $X$, are provided in Table 2. Gravity acts in the negative $z$-axis direction. We used a Legendre polynomial basis with 4th-order angular strains (torsion about $x$, bending about $y$ and $z$) and 2nd-order linear strains (elongation about $x$, shear in $y$ and $z$). Hence, the total DoF of the system is 24. For the numerical integration, we used 5 Gauss quadrature points. Including the computational points at the base and the tip, the total number of points $n_p = 7$.

<table><caption>**Table 2.** Cable routing for the cable-driven manipulator. $r(X)$ is the manipulator radius and $r_t = r(L)$.</caption><thead><tr><th>Cable Number</th><th>y-Coordinate</th><th>z-Coordinate</th></tr></thead><tbody><tr><td>1</td><td>0</td><td>r(X)</td></tr><tr><td>2</td><td>-√3/r(X)</td><td>-1/2r(X)</td></tr><tr><td>3</td><td>√3/r(X)</td><td>-1/2r(X)</td></tr><tr><td>4</td><td>r<sub>t</sub>sin(4π/X)</td><td>r<sub>t</sub>cos(4π/X)</td></tr><tr><td>5</td><td>r<sub>t</sub>sin(4π/X)</td><td>-r<sub>t</sub>cos(4π/X)</td></tr></tbody></table>

The manipulator is actuated by arbitrary time-varying inputs for 10 s, as shown in Figure 3(b). These inputs were selected to assess the robustness of the approach to different types of actuation. Specifically, the inputs for cables 1 and 3 are continuous sinusoidal functions, the input for cable 2 varies linearly, and the inputs for cables 4 and 5 consist of discontinuous step functions. We used ode15s of MATLAB as the ODE integrator, as it demonstrated the fastest computational speed for this simulation. The $\dot{x}$ for the ODE is computed using (4), while its Jacobian is computed using (10b), (11), (13) and analytical derivatives (27), (33), (38), (40a) and (40b). To validate and compare the impact of analytical derivatives, the ode15s solver was

executed both without (method 1) and with (method 2) their inclusion. For the Newmark-$\beta$ scheme (method 3), a time step $h$ of 0.002 s is used in this example. The time step $h$ for all examples in this paper was chosen as the largest value that ensures convergence of the solution while maintaining a positional mismatch at the millimeter level throughout the entire integration period. Throughout all simulations, the default tolerances in MATLAB were used for ode15s and fsolve: ‘RelTol’ of $10^{-3}$, and ‘AbsTol’ of $10^{-6}$. The PC specifications for all computations presented here are as follows: 13th Gen Intel(R) Core(TM) i9-13900HX, 2.20 GHz, 64.0 GB RAM. The MATLAB version used is R2024a.

Figure 4 shows the simulation results for three cases. The dynamic response is very close for all three methods. Figure 4(a) plots the states of the robot and the tip trajectory. The tip position mismatch computed using $||\Delta r_{tip}||$ is plotted in Figure 4(b). The maximum error is less than 30 $\mu$m between methods 1 and 2, while it is less than 1 mm between methods 2 and 3. Similarly, the state mismatch between the first two methods is on the order of $10^{-4}$, while the mismatch between the second and third methods is higher, as shown in Figure 4(c). The following metric is used to calculate the mismatch between the states of the body, obtained from two different methods, at a given time instant:

$$e_x = \frac{||\Delta x||}{||\text{avg}(x)||} \quad (41)$$

Without analytical derivatives, ode15s took 5.40 s, whereas, with analytical derivatives, it was completed in 1.25 s. This indicates that using analytical derivatives makes the computation more than 4 times faster. Simulation using the Newmark-$\beta$ method was completed in 24.43 s, which can be explained by the fact that this scheme has a fixed time-step in contrast to ode15s whose time-step is adaptive.

The analytical derivatives of the forward dynamics (FD) are validated as follows. Using the finite difference method, we calculated the numerical derivatives of $FD(q, \dot{q}, u)$ with respect to $q$ and $\dot{q}$. The $i^{th}$ column of the numerical Jacobians is computed using the forward finite difference method:

$$\left(\frac{\partial FD}{\partial q_i}\right)_{\text{num}} = \frac{FD(q + \epsilon n_i, \dot{q}, u) - FD(q, \dot{q}, u)}{\epsilon} \quad (42a)$$

$$\left(\frac{\partial FD}{\partial \dot{q}_i}\right)_{\text{num}} = \frac{FD(q, \dot{q} + \epsilon n_i, u) - FD(q, \dot{q}, u)}{\epsilon} \quad (42b)$$

where $n_i$ is the unit vector in the direction of the $i$-th component, and $\epsilon$ is a perturbation value set to $10^{-6}$. We used the forward finite difference method as it requires the least function evaluations to compute numerical derivatives, though this may result in lower accuracy compared to other methods.

Using a similar metric like (41), we computed the mismatch between analytical and numerical Jacobian of FD with respect to $q$ and $\dot{q}$ at every time step.

$$e_{\frac{\partial FD}{\partial q}} = \frac{\left\| \left(\frac{\partial FD}{\partial q}\right)_{\text{ana}} - \left(\frac{\partial FD}{\partial q}\right)_{\text{num}} \right\|}{\left\| \text{avg} \left(\frac{\partial FD}{\partial q}\right) \right\|} \quad (43)$$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Dynamic response of the CDM: (a) Snapshots of the manipulator at different times and the tip trajectory with and without using analytical derivatives (AD) of *FD*. (b) The mismatch between the tip positions ($r_{tip}$). (c) State mismatch. Mismatch between numerical and analytical derivatives of *FD* (d) with respect to *q* and (e) with respect to $\dot{q}$.</figcaption></figure>

A similar equation is used for $\frac{\partial FD}{\partial \dot{q}}$. Equation (43) provides a quantitative assessment of the accuracy of the analytically computed Jacobians. The results are plotted in Figure 4(d) and (e). Numerical values on the order of $10^{-6}$ and $10^{-8}$ indicate that their differences are insignificant. Computing the Jacobians numerically took, on average, 10.5 ms, while the analytical Jacobian required only 1.3 ms, making the analytical method more than 8 times faster.

For the static equilibrium simulation, the residue is computed using (17a), while the analytical Jacobian is obtained from (17b), (27), and (40a) with $\dot{q} = 0$ and $\ddot{q} = 0$. We performed 1000 static equilibrium simulations on the same system. In each simulation, we used random values of cable tensions in the range of 0 to 100 N. Figure 5(a) illustrates 10 equilibrium shapes from these 1000 simulations. Both methods (with and without the analytical Jacobian) yielded identical results. The mismatch in tip positions and static equilibrium solutions ($e_q$) using the two approaches are shown in Figure 5(b) and (c). Without analytical derivatives, the static simulations took an average of 28.2 ms, whereas, with analytical derivatives, it took 3.3 ms, making the computation more than 8 times faster.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Static simulation results: (a) Ten arbitrary equilibrium shapes. (b) The mismatch between tip positions. (c) The mismatch static equilibrium solutions.</figcaption></figure>

## 5 Extension to Hybrid Multi-body System

A hybrid multi-body system may incorporate rigid joints (up to 6 DoF), rigid links, soft links, branched chains, and closed-chain joints (Figure 6). The forward pass of RNEA goes through the whole multibody ($i = 1$ to $N$) and the backward pass from $i = N$ to 1. We propose the following rules to simplify the analysis of hybrid multi-body systems:

1. Each link is defined by a rigid joint and a body (rigid or soft). Accordingly, every rigid link has two computational points: one at $X = 0$ and another at $X = 1$ of the joint. These points are defined in the body's reference frame, typically located at the center

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Schematics of a generic hybrid multi-body system featuring soft and rigid links, rigid joints, branched chains, and closed-chain (CC) joints.</figcaption></figure>
of mass (CM), where the inertia matrix is computed. Soft links have $n_p + 1$ computational points: at $X = 0$ of the rigid joint and $n_p$ along the soft body. $X = 1$ of the joint coincides with $X = 0$ of the soft body.

2. For a rigid body, the inertial and gravitational forces act on the computational frame on the joint at $X = 1$. At $X = 0$ of a rigid joint, $\mathcal{F}_k = 0$, $\mathcal{M}_k = 0$ and $\mathcal{N}_k = 0$.

3. Let $g_f$ denote a fixed transformation between two adjacent frames (Figure 6). The recursive transformation rule, as outlined in equations (26), (32), and (37), applies to both the forward and backward transformations between these two frames when $g_{\alpha-1\alpha}$ or $g_{\alpha\alpha+1}$ is replaced by $g_f$.

4. In a branched chain system, the kinematic term $(\bullet)^B$ denotes the contribution of all joints from the global frame to the computational point $\alpha$, following a serial chain. Meanwhile, $(\bullet)^C$ and $(\bullet)^S$ encompass all the children's joints, including all branches ahead of point $\alpha$. The dashed arrows in Figure 6 illustrate these points.

5. Computation of $(\bullet)^B$ in the forward pass follows the same rule of the geometric Jacobian in (2b):

$$ (\bullet)_{i, \alpha}^{B} = \mathrm{Ad}_{g_{\alpha\alpha-1}} \left( (\bullet)_{i, \alpha-1}^{B} + [\mathbf{0}_{6 \times n_i^{-}} (\bullet)_{i, \alpha-1} \mathbf{0}_{6 \times n_i^{+}}] \right) \quad (44) $$

Hence, while the joint quantities, $S_\alpha$, $R_\alpha$, $Q_\alpha$, and $Y_\alpha$ have dimensions $\mathbb{R}^{6 \times n_{dof_i}}$, their $(\bullet)^B$ counterparts have dimensions $\mathbb{R}^{6 \times n_{dof}}$.

6. Computation of $(\bullet)^S$ in the backward pass follows a similar rule:

$$ (\bullet)_{i, \alpha}^{S} = [\mathbf{0}_{6 \times n_i^{-}} (\bullet)_{i, \alpha} \mathbf{0}_{6 \times n_i^{+}}] + \mathrm{Ad}^*_{g_{\alpha\alpha+1}} (\bullet)_{i, \alpha+1}^{S} \quad (45) $$

Hence, $(\bullet)^S$ also has dimensions $\mathbb{R}^{6 \times n_{dof}}$.

7. For the $i^{th}$ Cosserat rod, equations (25), (31), and (36) are generalized for the multi-body as follows:

$$ \frac{\partial I D_{i, \alpha}}{\partial \boldsymbol{q}} = \left[ \mathbf{0}_{n_{dof_i} \times n_i^{-}} \frac{\partial S_{i, \alpha}^T}{\partial q_i} \boldsymbol{\mathcal{F}}_{i, \alpha}^C \mathbf{0}_{n_{dof_i} \times n_i^{+}} \right] \quad (46a) $$

$$ + S_{i, \alpha}^T (\mathcal{N}_{i, \alpha}^C \mathbf{R}_{i, \alpha}^B + \mathcal{M}_{i, \alpha}^C \mathbf{Q}_{i, \alpha}^B + \mathbf{U}_{i, \alpha}^S + \mathbf{P}_{i, \alpha}^S) $$

$$ \frac{\partial I D_{i, \alpha}}{\partial \dot{\boldsymbol{q}}} = S_{i, \alpha}^T (\mathcal{N}_{i, \alpha}^C S_{i, \alpha}^B + \mathcal{M}_{i, \alpha}^C \mathbf{Y}_{i, \alpha}^B + \mathbf{V}_{i, \alpha}^S) \quad (46b) $$

$$ \frac{\partial I D_{i, \alpha}}{\partial \ddot{\boldsymbol{q}}} = S_{i, \alpha}^T (\mathcal{M}_{i, \alpha}^C S_{i, \alpha}^B + \mathbf{W}_{i, \alpha}^S) \quad (46c) $$

8. Finally, the partial derivatives of $ID$ and $\tau$ are computed by vertically concatenating row-blocks of each Cosserat rod according to:

$$ \frac{\partial I D}{\partial \boldsymbol{q}} = \left\|_{i=1}^{N} \sum_{\alpha=1}^{n_p-1} \frac{\partial I D_{i, \alpha}}{\partial \boldsymbol{q}} \right\| \quad (47a) $$

$$ \frac{\partial I D}{\partial \dot{\boldsymbol{q}}} = \left\|_{i=1}^{N} \sum_{\alpha=1}^{n_p-1} \frac{\partial I D_{i, \alpha}}{\partial \dot{\boldsymbol{q}}} \right\| \quad (47b) $$

$$ \frac{\partial I D}{\partial \ddot{\boldsymbol{q}}} = \left\|_{i=1}^{N} \sum_{\alpha=1}^{n_p-1} \frac{\partial I D_{i, \alpha}}{\partial \ddot{\boldsymbol{q}}} = \mathbf{M} \right\| \quad (47c) $$

$$ \frac{\partial \tau}{\partial \boldsymbol{q}} = \left\|_{i=1}^{N} \frac{\partial \tau_i}{\partial \boldsymbol{q}} \right\| \quad (47d) $$

$$ \frac{\partial \tau}{\partial \dot{\boldsymbol{q}}} = \left\|_{i=1}^{N} \frac{\partial \tau_i}{\partial \dot{\boldsymbol{q}}} \right\| \quad (47e) $$

where $\| \|$ is the vertical concatenation operator. In the backward pass from $N$ to 1, the corresponding derivatives of each Cosserat rod (as block matrices of size $\mathbb{R}^{n_{dof_i} \times n_{dof}}$) are concatenated from bottom to top.

Analytical derivatives of arbitrary hybrid multi-body systems can be evaluated following these rules. In addition to these, the dynamics of a multi-body system is often subject to equality constraints that can be expressed as $c(q, \dot{q}) = 0$. These constraints introduce an equal number of unknown Lagrange multipliers to enforce the equality condition. As a result, the system's ODE is transformed into a DAE. We explore two of the most commonly used types of constraints in robotic systems: joint coordinate controlled rigid joints and closed-chain systems.

4

## 5.1 Joint Coordinate Controlled Rigid Joints

Rigid joints can be actuated by specifying joint wrench $u$ or providing joint coordinates. The latter can be expressed as an equality constraint, with the Lagrange multiplier being $u$. Similar to the soft body actuation (39), the generalized actuation force of a rigid joint is given by $Bu$. The expression of the generalized actuation matrix ($B$) and the derivative of $Bu$ with respect to $q$ are provided in the Appendix D. For 1 DoF joints, $B$ is independent of $q$, and therefore, its derivative is $0$.

The computation of $FD$ and its derivatives remains the same if the joint is actuated by specifying $u$. However, if some of the joints are actuated by specifying joint coordinates, the $FD$ computation is different. One approach is to split $\ddot{q}$ into $\ddot{q}_u$ and $\ddot{q}_k$ and accordingly, $u$ into $u_k$ and $u_u$. The subscripts $k$ and $u$ stand for known and unknown quantities. Note that $\ddot{q}_k$ is the second time derivative of the specified joint coordinates $q_k = f(t)$, making it an index-1 DAE. Bringing the unknown terms to the LHS and known terms to the RHS, the generalized dynamics equation can be rewritten as:

$$ [M_u - B_u] \begin{bmatrix} \ddot{q}_u \\ \dot{u}_u \end{bmatrix} = [B_k - M_k] \begin{bmatrix} \mathbf{u}_k \\ \ddot{q}_k \end{bmatrix} + (\bullet) \quad (48) $$
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Hybrid serial robot consisting of 7 rigid links and 1 soft link with a fixed joint. The rigid links shown in red have revolute joints rotating about their local x-axis, while those in blue rotate about the y-axis. (b) Arbitrary joint angles as a function of time.</figcaption></figure>

where, the subscripts of M and B indicate corresponding columns of known and unknown quantities.

Equation (48) is used to solve $[\ddot{q}_u^T u_u^T]^T$ and FD is found by combining $\ddot{q}_u$ and $\ddot{q}_k$. To ensure numerical stability, during each iteration of the FD, the known joint angles and their derivatives are replaced with the specified $q_k$ and its time derivative $\dot{q}_k$. The derivative of (48) with respect to $q$ gives:

$$ \left[ \frac{\partial \ddot{q}_u}{\partial \dot{q}_u} \frac{\partial \dot{q}_u}{\partial u_u} \frac{\partial \dot{u}_u}{\partial q_u} \right] = [M_u - B_u]^{-1} \left( \frac{\partial \tau}{\partial q_u} - \frac{\partial ID}{\partial q_u} \right) \quad (49) $$

Finally, $\frac{\partial FD}{\partial q}$ is found by combining $\frac{\partial \ddot{q}_u}{\partial q_u}$ with $\frac{\partial \ddot{q}_k}{\partial q} = 0^{n_k \times n_{dof}}$ and $\frac{\partial \ddot{q}_u}{\partial q_k} = 0^{n_u \times n_k}$. The partial derivative of FD with respect to $\dot{q}$ is computed similarly.

The same problem can be solved using index-3 DAE formulation using the Newmark-$\beta$ scheme. In this case, we have $q_k = f(t)$ and the residue (15) and Jacobian (16) are modified to include the unknown $u_{u,n+1}$ and partial derivative with respect to $q_{u,n+1}$ and $u_{u,n+1}$.

$$ R_{dyn} = \tau(q_{n+1}, u_{n+1}) - ID(q_{n+1}) \quad (50a) $$

$$ J_{dyn} = \left[ \frac{\partial \tau}{\partial q_u} - \frac{\partial ID}{\partial q_u} \quad B_u \right] \quad (50b) $$

The statics problem involves finding the unknown generalized coordinates and joint forces ($q_u$ and $u_u$). It can be seen that the residue and Jacobian of the static equilibrium are identical to (50a) and (50b) with $\dot{q}$ and $\ddot{q}$ set to zero.

### 5.1.1 Hybrid Serial Robot
We modeled a hybrid serial robot by combining a rigid serial robotic arm with a soft manipulator, as shown in Figure 7(a). The robot consists of 7 revolute joints that are actuated by arbitrary joint angles according to Figure 7(b). The actuation inputs for the joint coordinate controlled rigid joints were selected as $C^2$-continuous functions since the joint angle, velocity, and acceleration are required for the computation of dynamics using index-1 DAE. Additionally, all joints start with an initial value of zero at $t = 0$. For the chosen arbitrary input, joints 1, 3, and 6 follow sinusoidal functions, joints 2 and 4 are governed by fourth-order polynomials, and joints 5 remain fixed throughout the simulation. The soft link has a

fixed joint, a length of 50 cm, and a radius linearly varying from 1.5 cm to 1 cm. The mechanical properties of the soft body are identical to that of the first example. The soft body is modeled as a Kirchhoff rod (no shear deformation) with a fourth-order strain basis. Hence, in total, the robot has 27 DoF.

The dynamic response of the robot is computed over 10 seconds using three methods: ode15s without providing analytical derivatives, ode15s with analytical derivatives, and the Newmark-$\beta$ method. For the Newmark-$\beta$ approach, a time step of 0.001 s was used, as the solution failed to converge for larger time steps. Snapshots of the simulation results for the first 5 seconds are displayed in Figure 8(a). Without analytical derivatives, the ode15s took 7.08 s, while with analytical derivatives, it was completed in 2.49 s, making it approximately 3 times faster. While for our custom Newmark-$\beta$ code, the simulation time was 3 min and 51 s. Mismatch metrics, similar to those in the first example, are shown in Figures 8(b), (c), (d), and (e). The position and state mismatch between the first two methods is in the order of 0.1 mm and $10^{-3}$, respectively. The differences between the numerical and analytical derivatives of FD indicate that their discrepancies are insignificant. Numerically computing the Jacobians took an average of 32.1 ms, whereas the analytical Jacobian required only 2.5 ms, making the analytical approach nearly 13 times faster. The mismatch metrics in Figure 8 appear noisier than those in Figure 4. This could be due to the more rapid dynamic motion of the serial robot compared to the cable-driven manipulator.

Static equilibrium simulation of the robot involves solving $q_u$ of the soft body and $u_u$ of the rigid joints. Similar to the first example, we performed 1000 static equilibrium simulations on the same system where we used random values of joint angles in the range of $-\pi/4$ to $\pi/4$. Figure 9(a) illustrates 5 equilibrium shapes from these 1000 simulations. Both methods (with and without the analytical Jacobian) yielded very close results. The mismatch in tip positions and static equilibrium solutions using the two approaches are shown in Figure 9(b) and (c). Without analytical derivatives, the static simulations took an average of 49.9 ms, whereas, with analytical derivatives, it took 6.3 ms, making the computation approximately 8 times faster.

## 5.2 Closed Chain Systems

Robots can have a closed-chain (CC) structure, as shown in Figure 6. Closed-chain systems can be modeled using kinematic constraints ($e(q) = 0$) expressed in the Pfaffian form (velocity level): $A(q)\dot{q} = 0$. Accordingly, their generalized dynamic equation is given by Armanini et al. (2022a):

$$ M\ddot{q} = \tau + F + A^T\lambda \quad (51a) $$

$$ A\ddot{q} + \dot{A}\dot{q} + \frac{2}{T_B}A\dot{q} + \frac{1}{T_B^2}e = 0 \quad (51b) $$

where $\lambda = [\lambda_1^T, \lambda_2^T, \dots, \lambda_{n_{CL}}^T]^T$ represents the unknown constraint wrenches, where $n_{CL}$ is the number of closed-chain joints and $T_B$ is the Baumgart stabilization constant.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Dynamic response of the serial robot: (a) Superimposed snapshots of the robot at different times and the tip trajectory with (blue) and without (red) using analytical derivatives. (b) Tip position mismatch. (c) State mismatch. The mismatch between numerical and analytical derivatives of FD (d) with respect to $q$ and (e) with respect to $\dot{q}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Static simulation results: (a) Five arbitrary equilibrium shapes. (b) The mismatch between tip positions. (c) The mismatch static equilibrium solutions ($q_u$ and $u_u$).</figcaption></figure>

$A$, $\dot{A}$, and $e$ are given by:

$$A = \left\|_{k=1}^{n_{CL}} \Phi_{\perp k}^{T} (\mathrm{Ad}_{g_{kBA}} \mathbf{J}_{kA} - \mathbf{J}_{kB}) \right. \quad (52a)$$

$$\dot{A} = \left\|_{k=1}^{n_{CL}} \Phi_{\perp k}^{T} (\mathrm{Ad}_{g_{kBA}} \mathrm{Ad}_{\eta_{kBA}} \mathbf{J}_{kA} + \mathrm{Ad}_{g_{kBA}} \dot{\mathbf{J}}_{kA} - \dot{\mathbf{J}}_{kB}) \right. \quad (52b)$$

$$e = \left\|_{k=1}^{n_{CL}} \Phi_{\perp k}^{T} (\log(g_{kBA}))^{\vee} \right. \quad (52c)$$

where $g_{kBA} = g_{kB}^{-1}g_{kA}$, $\Phi_{\perp k}$ is the basis for the constrained degrees of freedom and $\mathbf{J}_{kA}$ and $\mathbf{J}_{kB}$ are geometric Jacobians at points A and B where the closed-loop joint is defined.

Combining (51a) and (51b) we get:

$$M_A \begin{bmatrix} \ddot{q} \\ \lambda \end{bmatrix} = \begin{bmatrix} \tau + \mathbf{F} \\ -\dot{A}\dot{q} - \frac{2}{T_B} A\dot{q} - \frac{1}{T_B^2} e \end{bmatrix} \quad (53)$$

where,

$$M_A = \begin{bmatrix} M & -A^T \\ A & 0 \end{bmatrix} \quad (54)$$

Equation (53) is used to solve for $\ddot{q}$ and $\lambda$, and $\lambda$ is ignored for the time integration. The partial derivative of (53) with respect to $q$ gives

$$\begin{bmatrix} \frac{\partial FD}{\partial \dot{q}} \\ \frac{\partial \dot{\lambda}}{\partial q} \end{bmatrix} = M_A^{-1} \begin{bmatrix} \frac{\partial \tau}{\partial \dot{q}} - \frac{\partial ID(q, \dot{q}, \ddot{q}, \lambda)}{\partial q} \\ -\frac{\partial A}{\partial q} \ddot{q} - \frac{\partial \dot{A}}{\partial q} \dot{q} - \frac{2}{T_B} \frac{\partial A}{\partial q} \dot{q} - \frac{1}{T_B^2} \frac{\partial e}{\partial q} \end{bmatrix} \quad (55)$$

where $ID$ includes the contributions of the constraint force: a local wrench of $\mathrm{Ad}_{g_{kBA}} \Phi_{\perp k} \lambda_k$ on point A and $-\Phi_{\perp k} \lambda_k$ on point B of the closed-loop joint $k$. $ID$ treats $\lambda$ as an independent variable.

Similarly, the derivative of (53) with respect to $\dot{q}$ gives:

$$\begin{bmatrix} \frac{\partial FD}{\partial \dot{q}} \\ \frac{\partial \dot{\lambda}}{\partial \dot{q}} \end{bmatrix} = M_A^{-1} \begin{bmatrix} \frac{\partial \tau}{\partial \dot{q}} - \frac{\partial ID(q, \dot{q}, \ddot{q}, \lambda)}{\partial \dot{q}} \\ -\frac{\partial \dot{A}}{\partial \dot{q}} \dot{q} - \dot{A} - \frac{2}{T_B} A \end{bmatrix} \quad (56)$$

For the index-3 DAE formulation of the problem, we have additional unknowns $\lambda$ and additional equations given by the kinematic constraint $e(q) = 0$. Hence, the residue (15) and Jacobian (16) of the Newmark-β scheme are modified as follows.

$$\mathbf{R}_{dyn}(\mathbf{q}_{n+1}, \mathbf{\lambda}_{n+1}) = \begin{bmatrix} \tau(\mathbf{q}_{n+1}) - ID(\mathbf{q}_{n+1}, \mathbf{\lambda}_{n+1}) \\ e(\mathbf{q}_{n+1}) \end{bmatrix} \quad (57a)$$

$$\mathbf{J}_{dyn}(\mathbf{q}_{n+1}) = \begin{bmatrix} \frac{\partial \tau}{\partial \dot{q}} - \frac{\partial ID}{\partial \dot{q}} & \mathbf{A}^T \\ \frac{\partial e}{\partial \dot{q}} & 0 \end{bmatrix} \quad (57b)$$

The static equilibrium simulation involves solving $q$ and $\lambda$ where, $\dot{q}$ and $\ddot{q}$ are zero. In this case, the residue and Jacobian are identical to (57a) and (57b). Analytical derivatives of all quantities in (55), (56), and (57b) can be found in the Appendix D.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** (a) A fin-ray finger with 17 soft links and 6 closed-chain joints. The robot is actuated by a prismatic joint along the global z-axis at its base. (b) Displacement input of the prismatic joint in the range -2.5 to 2.5 cm.</figcaption></figure>

### 5.2.1 Fin-Ray Finger
Figure 10 shows a fin-ray finger, a multi-body system with 17 soft links, revolute joints, a prismatic joint, and 6 closed-chain joints. The finger is actuated by a displacement of the prismatic joint, as shown in Figure 10(b). The ribs (horizontal links) are connected to the sidewalls via revolute joints. Each link has a rectangular cross-section with a width of 1 mm. The sidewalls have a height of 2.5 cm and a length of 3 cm, while the ribs have a height of 1 cm. All soft links are modeled using linear bending, constant elongation, and shear strains, resulting in a planar multi-body with 74 DoF. The material used is Acrylonitrile-Butadiene-Styrene (ABS), with a Young's modulus of 2 GPa and a Poisson's ratio of 0.35, commonly used in 3D printing. Five of the six closed-chain joints are revolute, and the top joint is fixed. The problem includes both the joint coordinate constraint and the closed-chain constraint.

In this example, when the analytical Jacobian was not provided, the dynamics simulation using ode15s was practically stuck at 0.04 s. Other MATLAB ODE integrators, such as ode45 and ode113, progressed extremely slowly, completing only ~0.1 s in 12 hours. However, when ode15s was supplied with the analytical derivatives of FD, the simulation was completed in 16 min and 34 s. Meanwhile, the Newmark-β integrator (with h = 0.01 s) completed the same simulation in just 1 min and 9 s, making it the fastest DAE solver for this example. Figure 11(a) illustrates the dynamic response of the finger, which exhibits periodic behavior as expected. Figure 11(b) and (c) compare the tip position and the robot states, respectively. The tip position has a mismatch in the order of mm, while the state mismatch is higher. Using a lower value of time step can reduce this mismatch at the cost of increased computational effort. Figure 11(d) and (e) indicates that the discrepancies between the numerical and analytical derivatives of FD are insignificant. The numerical derivatives took an average of 498 ms to compute, while the analytical derivatives took 23.74 ms, making the analytical approach more than 20 times faster.

# 6 Examples of Common External Forces

In this section, we derive the analytical derivatives of various external loads to which robots are commonly subjected: point wrenches, contact forces, and hydrodynamic forces.

## 6.1 Point Wrench

The point wrench can be expressed in the local frame (as a follower load) or the global frame. In the former case, the partial derivative of $\mathcal{F}_k$ with respect to $q$ is 0. In our formulation, an external wrench expressed in the global frame $\mathcal{F}_{gk}$ must be transformed into the local frame for the projection into the space of generalized coordinates. The transformed wrench in the local frame is given by:

$$ \mathcal{F}_k = \mathrm{Ad}^{-*}_{g_{\theta k}} \mathcal{F}_{gk} \quad (58) $$

where, $g_{\theta k}$ is the rotational component of $g_k$.

The derivative of this term with respect to $q$ is given by:

$$ \frac{\partial \mathcal{F}_k}{\partial q} = -\overline{\mathrm{ad}}_{\mathcal{F}_k}^* I_\theta S_k^B \quad (59) $$

where, $I_\theta = \mathrm{diag}([1 1 1 0 0 0])$.

Accordingly, the partial derivative of $ID$ with respect to $q$ includes an additional term at every computational point $k$:

$$ \frac{\partial ID_k}{\partial q} = (\bullet) + S_k^T \left( \frac{\partial \mathcal{F}_k^C}{\partial q} \right) \quad (60) $$

where, $\frac{\partial \mathcal{F}_k^C}{\partial q}$ is computed similarly to $\mathcal{F}_k^C$:

$$ \frac{\partial \mathcal{F}_k^C}{\partial q} = \mathrm{Ad}^*_{g_{kk+1}} \left( \frac{\partial \mathcal{F}_{k+1}}{\partial q} + \frac{\partial \mathcal{F}_{k+1}}{\partial q}^C \right) \quad (61) $$

### 6.1.1 Hybrid Parallel Robot
Figure 12(a) shows an example of a hybrid parallel robot with three soft pillars and a rigid platform. The properties of the soft body are kept identical to those of CDM. The soft body is modeled as a cuboid with dimensions 1.5 cm × 3 cm × 15 cm. The top rigid platform is an equilateral triangle with a side length of 20 cm and a height of 1 cm. The soft pillars are connected to the platform via spherical joints and are parameterized using cubic angular strains and first-order linear strains. Hence, the total DoF of the robot is 63. The dynamic simulation is solved using FD computed from (53), with analytic Jacobians according to (55) and (56). The platform is subjected to an external point wrench (force and moment) in the global frame according to Figure 12(b), and the robot's dynamic response is calculated over 10 seconds.

Figure 13(a) shows snapshots of the dynamic response. In this example, when the analytical Jacobian was not provided, the simulation using ode15s advanced at an impractically slow pace. Consequently, we switched to MATLAB's ode45 and ode113 (method 1) integrators. The simulation took 1 hr 46 mins with ode45 and 1 hr 15 mins with ode113. In contrast, ode15s with the analytical Jacobian completed the simulation in just 2.46 s, making it about 1800 times faster than ode113. For the Newmark-β scheme, the computational time was 17.62 s for a time step of 0.01 s. Figures 13(b) and (c) demonstrate how closely the dynamic simulation results of all three methods align. The position and state
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Dynamic response of the fin-ray finger: (a) Snapshots of the finger at different times with the trails of the tip position. (b) Average values of mismatch between tip positions. (c) State mismatch. The mismatch between numerical and analytical derivatives of FD (d) with respect to *q* and (e) with respect to *q̇*.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** (a) Hybrid parallel robot consisting of three soft pillars, a rigid platform with three spherical joints, and two closed-chain revolute joints subject to a point wrench (force *f* and moment *m*). (b) Applied external wrenches.</figcaption></figure>

mismatch between the first two methods is in the order of 10⁻⁵ m and 10⁻⁴, respectively. The difference between the analytical and numerical derivatives of FD, displayed in Figures 13(d) and (e), validates the analytical derivatives. Numerically computing the derivatives took an average of 82.9 ms, whereas the analytical derivatives required only 5.1 ms, making the analytical approach nearly 16 times faster. Thus, it is not just the speed of computation; its accuracy also contributed to the significant improvement in the overall computational efficiency of dynamic simulations.

For the static equilibrium simulations, the system is subjected to 1000 randomly applied point wrenches: forces in the range -0.5 to 0.5 N and moments in the range -0.05 to 0.05 Nm. The equilibrium values of *q* and *λ* are computed by solving (57a). When the analytical derivatives (57b) were not provided, the simulation took an average of 166.39 ms. In contrast, the simulation was completed in just 8.97 ms with the analytical derivatives, making it approximately 18 times

faster. Figure 14(a) presents four static equilibrium shapes from the simulation. The average corner position mismatch is depicted in Figure 14(b), while the solution mismatch, including the equilibrium values of *q* and *λ*, is shown in Figure 14(c). These figures suggest that, with and without analytical derivatives, the solutions obtained are identical, except in cases where multiple solutions are present.

## 6.2 Contact Force

Accurate modeling of contact mechanics is paramount in robotics and has been addressed through various models and algorithms Lidec et al. (2024). Here, we examine a simplified internal contact scenario where a cable-driven soft manipulator, such as the one depicted in Figure 3, is actuated inside a hollow cylinder, with the cylinder's centerline aligned to the global z-axis. The following assumptions are made to simplify the contact force:

1. To compute the contact force, the manipulator is discretized into spheres at each computational point ($n_p = 12$ in this example). The radius of each sphere corresponds to the manipulator's radius at that point $k$.

2. Only normal forces are considered; tangential (frictional) forces are assumed to be zero. Hence, the contact force does not cause any local moments.

3. Penetration into the cylinder is allowed, and the contact force is assumed to be a function of the penetration $\delta$.

4. The direction of the contact force is normal to the surface of the cylindrical wall, taking the form $[n_x, n_y, 0]^T$.

To model the contact force, we used the Hertz model Flores (2022), a two-parameter non-linear contact model given by:
<figure><img src="image_16.png" /><figcaption>Figure 13. Dynamic response of the hybrid parallel robot: (a) Snapshots of the robot at different times with the trails of the corner positions of the rigid platform. Solid and dotted arrows indicate the applied point forces and moments, respectively. (b) Average values of mismatch between corner positions. (c) State mismatch. The mismatch between numerical and analytical derivatives of FD (d) with respect to q and (e) with respect to $\dot{q}$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 14. Static simulation results of the hybrid parallel robot: (a) Four arbitrary equilibrium shapes. (b) The mismatch between tip positions. (c) The mismatch static equilibrium solutions (q and λ). Vertical lines in (b) and (c) indicate cases of multiple static solutions obtained from both methods.</figcaption></figure>

$$ \mathcal{F}_{gk} = \begin{cases} \begin{pmatrix} \mathbf{0} \\ k_c \delta_k^p \mathbf{u}_{\perp k} \end{pmatrix} & \text{if } \delta > 0 \\ \mathbf{0} & \text{otherwise} \end{cases} \qquad (62) $$

where $k_c$ represents the contact stiffness coefficient and $p$ is the force exponent, typically dependent on the material properties and geometry of the contact. $\mathbf{u}_{\perp}$ is the unit normal from the center of the sphere to the wall of the cylinder. The penetration $\delta_k$ is given by:

$$ \delta_k = \|\mathbf{n}_{\perp k}\| + r_k - r_{cyl} \qquad (63) $$

where $\mathbf{n}_{\perp k} = C_1 r_{gk}$, $C_1 = \text{diag}([1 \ 1 \ 0])$ and $r_{gk}$ is the vector from the global frame to point $k$. Note that $\mathbf{u}_{\perp k} = \mathbf{n}_{\perp k} / \|\mathbf{n}_{\perp k}\|$. Since $\mathcal{F}_{gk}$ is in the global frame, we need to transform it into the local frame according to (58). The derivative of the local force with respect to $q$ is given by:

$$ \frac{\partial \mathcal{F}_k}{\partial q} = -\overline{\text{ad}}_{\mathcal{F}_k}^* I_\theta S_k^B + \text{Ad}_{g_{\theta k}}^{-1} \frac{\partial \mathcal{F}_{gk}}{\partial q} \qquad (64) $$

The partial derivative of the contact force in the global frame is given by:

$$ \frac{\partial \mathcal{F}_{gk}}{\partial q} = \begin{pmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & k^* C_1 \end{pmatrix} \text{Ad}_{g_{\theta k}} S_k^B \qquad (65) $$

where $k^*$ is given by

$$ k^* = k_c p \delta_k^{p-1} \mathbf{u}_{\perp k} \mathbf{u}_{\perp k}^T + \frac{k_c \delta_k^p}{\|\mathbf{n}_{\perp k}\|} (I_3 - \mathbf{u}_{\perp k} \mathbf{u}_{\perp k}^T) \qquad (66) $$

Similar to the case of point force, the derivative of the local force (64) is projected backward, and their contribution is accounted into the partial derivative of ID. Note that the penetration condition of the contact force (62) also applies to its derivative. Readers interested in detailed derivation may refer to Appendix D.

### 6.2.1 Contact Dynamics
For the dynamic simulation, we used the same manipulator and actuation inputs as shown in Figure 3. We used an inner wall radius of $r_{cyl}$ = 15 cm, $k_c$ = $10^5$ N/m, and $p$ = 1.5. The simulation results are displayed in Figure 15(a). Simulation using ode15s without providing
<figure><img src="image_18.png" /><figcaption>Figure 15. Contact simulation results: (a) Snapshots of dynamics of CDM inside a hollow cylinder with tip trajectories. Inset showing a top view. (b) Mismatch between tip positions. (c) State mismatch. Mismatch between numerical and analytical derivatives of FD: (d) with respect to $q$ and (e) with respect to $\dot{q}$.</figcaption></figure>

analytical derivatives was completed in 30.61 s, while with analytical Jacobian, it was completed in 7.80 s, making it nearly 4 times faster. Integration using the Newmark-$\beta$ approach with a time step of 0.002 s was completed in 46.64 s. The validation between the methods in terms of tip position and robot state are shown in Figure 15(b) and (c), respectively. Between methods 1 and 2, the tip position mismatch is less than 40 $\mu$m and the state mismatch is in the order of $10^{-3}$. Figure 15(d) and (e) compare the derivatives of FD obtained using numerical and analytical methods. While the discrepancies are generally small, the relatively larger mismatches in the partial derivatives of FD with respect to $q$ can be attributed to the penetration condition. The average time to compute the numerical derivatives was 23.7 ms, while for the analytical derivatives, it was 2.3 ms, making the analytical method approximately 10 times faster.

## 6.3 Hydrodynamic Force

Several hybrid soft robots are designed for underwater exploration. In addition to gravity, they are subject to external forces due to buoyancy, drag, lift, and added mass (fluid displacement). Buoyancy acts as an acceleration opposing gravity, while the added mass force can be simplified as additional inertia that the body experiences. To model the effect of buoyancy, we can use the modified acceleration $a'_{G} = (1 - \rho_{w}/\rho_{b})a_{G}$, where $\rho_{w}$ is the density of the surrounding water and $\rho_{b}$ is the density of the robot's body. To tackle the hydrodynamic forces along rods, we follow Boyer et al. (2006) and use a model that combines a simplified version of the reactive theory of Lighthill (1970), with the resistive empirical model of Morison et al. (1950). In this model, the rod inertia matrices are replaced by the modified inertia matrix $M'_{k} = M_{k} + M_{Ak}$, where $M_{Ak}$ represents the added inertia matrix on the k-th computational point on the robot. The effects of viscosity are modeled by a field of drag-lift forces applied along the flagellum, and given by Armanini et al. (2021):

$$ \mathcal{F}_{Dk} = \mathcal{D}_{k} \|v_{k}\| \eta_{k} \quad (67) $$

where, $\mathcal{D}_{k}$ is the drag-lift matrix and $\|v_{k}\| = \sqrt{\eta_{k}^{T} I_{v} \eta_{k}}$, where $I_{v} = \text{diag}([0 \ 0 \ 0 \ 1 \ 1 \ 1])$. Hence, we have:

$$ \frac{\partial \mathcal{F}_{Dk}}{\partial \dot{q}} = \mathcal{D}_{k}^{*} \frac{\partial \eta_{k}}{\partial q} \quad (68) $$

where $\mathcal{D}_{k}^{*} = \mathcal{D}_{k} (\frac{1}{\|v_{k}\|} \eta_{k} \eta_{k}^{T} I_{v} + \|v_{k}\| I_{6})$.

Substituting the partial derivative of the velocity twist (Appendix D), we get the derivative of the drag-lift force:

$$ \frac{\partial \mathcal{F}_{Dk}}{\partial \dot{q}} = \mathcal{D}_{k}^{*} \sum_{\beta<k} \text{Ad}_{g_{\beta k}}^{-1} R_{\beta} \quad (69) $$

Similarly, we can write the partial derivative of the drag-lift force with respect to $\dot{q}$ as follows,

$$ \frac{\partial \mathcal{F}_{Dk}}{\partial \dot{q}} = \mathcal{D}_{k}^{*} \sum_{\beta<k} \text{Ad}_{g_{\beta k}}^{-1} S_{\beta} \quad (70) $$

With the form of (69) and (70), it is easy to see that the contribution of drag-lift force can be accounted for by including an additional term in (24a):

$$ \mathcal{N}_{k} = \overline{\text{ad}}_{\mathcal{M}_{k} \eta_{k}}^{*} + \text{ad}_{\eta_{k}}^{*} \mathcal{M}_{k} - \mathcal{M}_{k} \text{ad}_{\eta_{k}} + \mathcal{D}_{k}^{*} \quad (71) $$

### 6.3.1 Underwater Flagellated Vehicle
Using this, we simulated the dynamics of an underwater (UW) flagellated vehicle. The schematic of the robot is shown in Figure 16. The robot is a hybrid-branched chain system with a mobile (6 DoF) body, two shafts (rigid bodies with revolute joints), six soft bodies (flagella), and six hooks that connect the flagella with the shafts. The three front flagella are 35 cm long, while the rear ones are 50 cm long. Each flagellum has a base radius of 12.5 cm, tapering smoothly to a point at the tips. The material properties of the soft body are kept identical to previous examples. The robot is actuated by joint torques while subject to external hydrodynamic forces. In this example, we assume that all the bodies are neutrally buoyant ($\rho_{w} = \rho_{b}$). Each flagellum is modeled
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Schematic of the underwater hybrid mobile robot. The robot is actuated by joint torque and is subject to drag-lift and added mass forces.</figcaption></figure>

as an inextensible Kirchhoff rod with a cubic strain field, resulting in 80 DoFs for the robot.

We applied a joint torque of 0.25 Nm on the first shaft and -0.25 Nm on the second to cancel out the rotation of the robot body. The dynamics of the robot are computed for 10 s. Similar to the case of the hybrid parallel robot, the simulation using *ode15s* without providing the Jacobian progressed very slowly. As a result, we used *ode45* and *ode113* (method 1), obtaining simulation times of 11 min 58 s and 8 min 16 s, respectively. In contrast, *ode15s* with the analytical Jacobian completed the simulation in just 1.17 s, making it more than 400 times faster. Using the Newmark-β method, the simulation was finished in 60.91s. The dynamic response of the robot is displayed in Figure 17(a). The mismatch metrics comparing the outputs of all three methods are provided in Figure 17(b) and (c). Between methods 1 and 2, the state mismatch is in the order of $10^{-4}$. The validation of the analytical derivatives of *FD* is provided in Figure 17(d) and (e). The numerical method took 216.4 ms on average, while the average time for the analytical approach was 11.9 ms, making it nearly 18 times faster.

# 7 Discussions and Conclusions

In this work, we used the RNEA algorithm to derive the analytical derivatives of the GVS model for the mechanical analysis of hybrid soft-rigid robots. The “pseudo rigid joint” formulation of the GVS model, facilitated by the Magnus expansion of soft rods, allows for the extension and generalization of analytical derivatives of rigid body dynamics algorithms Carpentier and Mansard (2018); Singh et al. (2022). These derivatives were applied to the dynamic and static simulations of various hybrid multi-body systems. We presented six dynamic and three static simulations, demonstrating significant computational improvements. For dynamic simulations, three integration methods were used: method 1 employed either *ode15s* or *ode113* (when *ode15s* was slower or stalled), method 2 utilized *ode15s* with analytical derivatives, and method 3 applied the Newmark-β approach. Method 2 provided exceptional speedup, particularly in multi-body systems with constraints, where speedup factors of over 3 orders of magnitude ($> 10^3$) were observed. The Newmark-β approach was observed to be

faster for the highly constrained fin-ray finger simulation, while the solution using method 1 was impractically slow for the same case. The results for dynamic simulations summarizing the increase in simulation speed and the mismatch between the analytical and numerical Jacobian of *FD* are shown in Table 3. The results emphasize the effectiveness of analytical derivatives and the accuracy of its computation. For static simulations, *fsolve* with analytical Jacobian was at least 8 times faster in all cases due to the improved convergence achieved with the accurate derivative information and efficient recursive implementation.

A limitation of the *ode15s* solver in MATLAB is that it does not allow passing $\ddot{x}$ to the Jacobian function. Since the partial derivative of *FD* (11) requires $\ddot{q}$, we need to compute $\ddot{q}$ twice, slowing down the simulation. Adding the capability to directly pass $\ddot{x}$ could improve simulation speed. Even though the Newmark-β approach is not a built-in MATLAB function, and our implementation is currently far from optimized, it offers several advantages. Its high stability allows large time steps without sacrificing accuracy, especially in handling DAE index-3 problems, as they are commonly considered in nonlinear structural dynamics of hybrid soft-rigid systems, where this formulation is generally preferred to the index-1 formulation of the Baumgart algorithm. In addition, the Newmark-β scheme works particularly well for highly constrained systems. Moreover, preserving the symplectic structure of mechanical systems while maintaining simulation stability and second-order accuracy, it is ideal for long-duration simulations. Here, we used a fixed time step that resulted in a positional mismatch at the millimeter level compared to method 2, however, developing a variable time-step version of the method could further enhance its efficiency and performance. Similarly, replacing the nominal version of the Newmark scheme with the generalized HHT or $\alpha$ schemes should improve its performance, while requiring only minor modifications to the scheme.

We have implemented the GVS model in a GUI-based MATLAB toolbox called SoRoSim Mathew et al. (2022a). The toolbox allows users from various backgrounds to easily simulate hybrid multi-body systems without requiring deep expertise in underlying algorithms. Implementation of analytical derivatives for the analysis of arbitrary robotic systems is a complex task, particularly for new researchers. To further simplify this process, we updated SoRoSim into a fully differentiable simulator, using the theory presented in this work. The latest version of the toolbox can be obtained from MathWorks or GitHub Mathew et al. (2025). The readers are encouraged to explore the toolbox and leverage its differential capabilities for various robotic analyses.

The applications of analytical derivatives extend far beyond the algorithms presented for static and dynamic analysis. They are foundational for a wide range of advanced control and optimization techniques. For example, in design optimization, they enable efficient tuning of robot geometry and material properties to achieve desired performance. In trajectory optimization, they provide efficient gradient information for motion planning, improving computational efficiency and accuracy. They are also crucial for Model Predictive Control (MPC), where real-time system predictions are needed to maintain optimal performance. Moreover,
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Dynamic response of the underwater mobile robot: (a) Snapshots of the robot at different times with the trajectory of two flagella. (b) Average values of mismatch between the tips of all six flagella. (c) State mismatch. The mismatch between numerical and analytical derivatives of FD (d) with respect to *q* and (e) with respect to *dot q*.</figcaption></figure>

<table><caption>**Table 3.** Summary of computational times and the maximum Jacobian mismatch for 10-second dynamic simulations. Computational time is denoted as $t_c$, and the speed-up factor is $S_f$. Method 1 does not use analytical derivatives, while Method 2 refers to ode15s with analytical derivatives, and Method 3 represents Neumark-β with analytical derivatives and a fixed time step $h$.</caption><thead><tr><th rowspan="2">Example</th><th rowspan="2">N</th><th rowspan="2">DoF</th><th rowspan="2">n<sub>p</sub></th><th colspan="2">Method 1</th><th colspan="2">Method 2</th><th colspan="2">Method 3</th><th rowspan="2">max(e<sub>∂FD</sub>/∂q)</th><th rowspan="2">max(e<sub>∂FD</sub>/∂dot q)</th></tr><tr><th>t<sub>c</sub></th><th>Solver</th><th>t<sub>c</sub></th><th>S<sub>f</sub></th><th>h(s)</th><th>t<sub>c</sub></th></tr></thead><tbody><tr><td>CDM</td><td>1</td><td>24</td><td>7</td><td>5.4 s</td><td>ode15s</td><td>1.25 s</td><td>4.32</td><td>0.002</td><td>24.43 s</td><td>2.7 &times; 10<sup>-6</sup></td><td>2.4 &times; 10<sup>-8</sup></td></tr><tr><td>Serial Robot</td><td>8</td><td>27</td><td>22</td><td>7.08 s</td><td>ode15s</td><td>2.49 s</td><td>2.84</td><td>0.001</td><td>3 min 51 s</td><td>8.5 &times; 10<sup>-7</sup></td><td>6.4 &times; 10<sup>-8</sup></td></tr><tr><td>Fin-ray Finger</td><td>23</td><td>74</td><td>136</td><td>Very slow</td><td>N/A</td><td>16 min 34 s</td><td>Very high</td><td>0.01</td><td>1 min 9 s</td><td>1.2 &times; 10<sup>-5</sup></td><td>9.2 &times; 10<sup>-8</sup></td></tr><tr><td>Parallel Robot</td><td>6</td><td>63</td><td>26</td><td>1 hr 15 min</td><td>ode113</td><td>2.46 s</td><td>1829.27</td><td>0.01</td><td>17.62 s</td><td>1.1 &times; 10<sup>-5</sup></td><td>6.7 &times; 10<sup>-8</sup></td></tr><tr><td>Contact</td><td>1</td><td>24</td><td>12</td><td>30.61 s</td><td>ode15s</td><td>7.8 s</td><td>3.92</td><td>0.002</td><td>46.64 s</td><td>8.9 &times; 10<sup>-4</sup></td><td>6.7 &times; 10<sup>-6</sup></td></tr><tr><td>UW Vehicle</td><td>19</td><td>80</td><td>80</td><td>8 min 16 s</td><td>ode113</td><td>1.17 s</td><td>424.27</td><td>0.01</td><td>1 min 9 s</td><td>1.1 &times; 10<sup>-6</sup></td><td>1.1 &times; 10<sup>-5</sup></td></tr></tbody></table>

their use can significantly improve methods like Differential Dynamic Programming (DDP), especially in handling non-linear constraints and complex robot configurations. Analytical derivatives facilitate accurate local model linearization of nonlinear system dynamics, which is essential for control design and stability analysis. They also enable sensitivity analysis by capturing system behavior changes due to parameter variations, a key aspect of robust design and uncertainty quantification. Furthermore, they enhance system identification by improving parameter estimation to minimize the sim-to-real gap. In reinforcement learning (RL), analytical derivatives aid policy optimization, accelerating convergence and improving training stability. Analytical derivatives can make several rigid-body algorithms accessible for soft robotic systems, allowing researchers to simulate and control complex hybrid systems using existing algorithms. The techniques presented in this paper offer broad applicability across the robotics field for enhancing the efficiency and precision of simulations. We believe that this work will play a pivotal role in addressing challenges in soft and hybrid soft-rigid robots, driving innovation, and enabling real-world applications.

## Acknowledgments

By the Khalifa University of Science and Technology under Grants RIG-2023-048 and RC1-2018-KUCARS, as well as the French ANR COSSEROOTS (ANR-20-CE33-0001).

## Appendices

### A Basic SE(3) Formulae

Adjoint operator of $s\epsilon(3)$:

$$ad_v = \begin{pmatrix} \tilde{\mathbf{w}} & \mathbf{0}_{3 \times 3} \\ \tilde{\mathbf{v}} & \tilde{\mathbf{w}} \end{pmatrix} \in \mathbb{R}^{6 \times 6}$$

where, $\mathbf{V} = [\mathbf{w}^T \mathbf{v}^T]^T \in \mathbb{R}^6$ is a screw vector.
Coadjoint operator of $s\epsilon(3)$:

$$ad_{\dot{v}}^* = \begin{pmatrix} \tilde{\mathbf{w}} & \tilde{\mathbf{v}} \\ \mathbf{0}_{3 \times 3} & \tilde{\mathbf{w}} \end{pmatrix} \in \mathbb{R}^{6 \times 6}$$

Coadjointbar operator of $s\epsilon(3)$:

$$\overline{ad_{\dot{v}}^*} = - \begin{pmatrix} \tilde{\mathbf{w}} & \tilde{\mathbf{v}} \\ \tilde{\mathbf{v}} & \mathbf{0}_{3 \times 3} \end{pmatrix} \in \mathbb{R}^{6 \times 6}$$

Adjoint map of $SE(3)$:

$$Ad_{g(X)} = \begin{pmatrix} \mathbf{R} & \mathbf{0}_{3 \times 3} \\ \tilde{\mathbf{r}}\mathbf{R} & \mathbf{R} \end{pmatrix} \in \mathbb{R}^{6 \times 6}$$
CoAdjoint map of SE(3):

$$ \mathrm{Ad}_{g(X)}^{*} = \begin{pmatrix} \mathbf{R} & \tilde{\mathbf{r}}\mathbf{R} \\ \mathbf{0}_{3 \times 3} & \mathbf{R} \end{pmatrix} \in \mathbb{R}^{6 \times 6} $$

Exponential map of SE(3) (replace $\Omega$ with $\xi$ for rigid joints):

$$ \exp(\hat{\Omega}) = I_4 + \hat{\Omega} + \frac{1}{\theta^2}(1 - \cos(\theta))\hat{\Omega}^2 + \frac{1}{\theta^3}(\theta - \sin(\theta))\hat{\Omega}^3 $$

where, $\theta = \sqrt{\Omega^T I_6 \Omega}$, $I_\theta = \mathrm{diag}([1 \ 1 \ 1 \ 0 \ 0 \ 0])$

Tangent operator and its derivative

$$ T(\Omega) = I_6 + \sum_{r=1}^{4} f_r(\theta) \mathrm{ad}_{\Omega}^{r} $$

$$ \dot{T}(\Omega, \dot{\Omega}) = \sum_{r=1}^{4} (f'_{r}(\theta) \dot{\theta} \mathrm{ad}_{\Omega}^{r} + f_{r}(\theta) (\mathrm{ad}_{\Omega}^{r})) $$

where, $\dot{\theta} = \frac{1}{\theta} \Omega^T I_\theta \dot{\Omega}$ and $(\mathrm{ad}_{\Omega}^r) = \sum_{u=1}^{r} \mathrm{ad}_{\Omega}^{r-u} \mathrm{ad}_{\Omega} \mathrm{ad}_{\Omega}^{u-1}$

$$ f_1(\theta) = \frac{1}{2\theta^2} (4 - 4 \cos(\theta) - \theta \sin(\theta)) $$

$$ f_2(\theta) = \frac{1}{2\theta^3} (4\theta - 5 \sin(\theta) + \theta \cos(\theta)) $$

$$ f_3(\theta) = \frac{1}{2\theta^4} (2 - 2 \cos(\theta) - \theta \sin(\theta)) $$

$$ f_4(\theta) = \frac{1}{2\theta^5} (2\theta - 3 \sin(\theta) + \theta \cos(\theta)) $$

$$ f'_1(\theta) = \frac{1}{2\theta^3} (-8 + (8 - \theta^2) \cos(\theta) + 5\theta \sin(\theta)) $$

$$ f'_2(\theta) = \frac{1}{2\theta^4} (-8\theta + (15 - \theta^2) \sin(\theta) - 7\theta \cos(\theta)) $$

$$ f'_3(\theta) = \frac{1}{2\theta^5} (-8 + (8 - \theta^2) \cos(\theta) + 5\theta \sin(\theta)) $$

$$ f'_4(\theta) = \frac{1}{2\theta^6} (-8\theta + (15 - \theta^2) \sin(\theta) - 7\theta \cos(\theta)) $$

# B Identities

## Screw Theory Identities

$$ \mathrm{ad}_{\mathbf{V}}^{*} = -\mathrm{ad}_{\mathbf{V}}^{T} \qquad (\mathrm{B1}) $$

$$ \mathrm{ad}_{\mathbf{V}_1} \mathbf{V}_2 = -\mathrm{ad}_{\mathbf{V}_2} \mathbf{V}_1 \qquad (\mathrm{B2}) $$

$$ \mathrm{ad}_{\mathbf{V}_1}^{*} \mathbf{V}_2 = \overline{\mathrm{ad}_{\mathbf{V}_2}^{*} \mathbf{V}_1} \qquad (\mathrm{B3}) $$

$$ \mathrm{Ad}_{g} \mathrm{ad}_{\mathbf{V}} = \mathrm{ad}_{\mathrm{Ad}_{g} \mathbf{V}} \mathrm{Ad}_{g} \qquad (\mathrm{B4}) $$

$$ \mathrm{ad}_{\mathbf{V}} \mathrm{Ad}_{g} = \mathrm{Ad}_{g} \mathrm{ad}_{\mathrm{Ad}_{g}^{-1} \mathbf{V}} \qquad (\mathrm{B5}) $$

$$ \mathrm{Ad}_{g}^{*} \mathrm{ad}_{\mathbf{V}} = \mathrm{ad}_{\mathrm{Ad}_{g} \mathbf{V}}^{*} \mathrm{Ad}_{g} \qquad (\mathrm{B6}) $$

$$ \mathrm{ad}_{\mathbf{V}}^{*} \mathrm{Ad}_{g}^{*} = \mathrm{Ad}_{g}^{*} \mathrm{ad}_{\mathrm{Ad}_{g}^{-1} \mathbf{V}}^{*} \qquad (\mathrm{B7}) $$

$$ \mathrm{Ad}_{g}^{*} \overline{\mathrm{ad}_{\mathbf{V}}^{*}} = \overline{\mathrm{ad}_{\mathrm{Ad}_{g}^{*} \mathbf{V}}^{*}} \mathrm{Ad}_{g} \qquad (\mathrm{B8}) $$

$$ \overline{\mathrm{ad}_{\mathbf{V}}^{*}} \mathrm{Ad}_{g} = \mathrm{Ad}_{g}^{*} \overline{\mathrm{ad}_{\mathrm{Ad}_{g}^{*} \mathbf{V}}^{*}} \qquad (\mathrm{B9}) $$

$$ \dot{\mathrm{Ad}}_{g} = \mathrm{Ad}_{g} \mathrm{ad}_{\eta} \qquad (\mathrm{B10}) $$

$$ \dot{\mathrm{Ad}}_{g}^{*} = \mathrm{Ad}_{g}^{*} \mathrm{ad}_{\eta}^{*} \qquad (\mathrm{B11}) $$

$$ \dot{\mathrm{Ad}}_{g}^{-1} = -\mathrm{ad}_{\eta} \mathrm{Ad}_{g}^{-1} \qquad (\mathrm{B12}) $$

$$ \dot{\mathrm{Ad}}_{g}^{-*} = -\mathrm{ad}_{\eta}^{*} \mathrm{Ad}_{g}^{-*} \qquad (\mathrm{B13}) $$

$$ \exp(\dot{\hat{\xi}}) = \exp(\hat{\xi}) \left( \mathrm{Ad}_{\exp(\hat{\xi})}^{-1} \mathrm{T}(\hat{\xi}) \dot{\hat{\xi}} \right) \qquad (\mathrm{B14}) $$

$$ \dot{g} = g\hat{\eta} \qquad (\mathrm{B15}) $$

## Summation Identities

For $M$ and $N$ that are matrix functions,

$$ \sum_{i=1}^{n} \sum_{j=i}^{n} M_i N_j = \sum_{j=1}^{n} \left( \sum_{i=1}^{j} M_i \right) N_j \qquad (\mathrm{B16}) $$

$$ \sum_{i=1}^{n} \sum_{j=1}^{i-1} M_i N_j = \sum_{j=1}^{n-1} \left( \sum_{i=j+1}^{n} M_i \right) N_j \qquad (\mathrm{B17}) $$

$$ \sum_{i=k+1}^{n} \sum_{j=k}^{i-1} M_i N_j = \sum_{j=k}^{n-1} \left( \sum_{i=j+1}^{n} M_i \right) N_j \qquad (\mathrm{B18}) $$

$$ \sum_{i=k+1}^{n} \sum_{j=1}^{i-1} M_i N_j = \left( \sum_{i=k+1}^{n} M_i \right) \left( \sum_{j=1}^{k-1} N_j \right) \qquad (\mathrm{B19}) $$

$$ + \sum_{j=k}^{n-1} \left( \sum_{i=j+1}^{n} M_i \right) N_j $$

# C Derivatives of GVS Elements

## Magnus Expansion and its Derivatives

Since a generic Cosserat strain field is non-constant along the material domain, we used the approximate form of Magnus expansion given by the fourth-order Zannah collocation approach Zanna (1999) to compute the kinematic and differential kinematic map from one computational point to the next. According to this, the Magnus expansion of $\xi$ between $X_\alpha$ and $X_{\alpha+1}$ is given by:

$$ \Omega_{\alpha} = \frac{h_{\alpha}}{2} (\xi_{\alpha z1} + \xi_{\alpha z2}) + \frac{\sqrt{3}h_{\alpha}^2}{12} \mathrm{ad}_{\xi_{\alpha z1}} \xi_{\alpha z2} \qquad (\mathrm{C1}) $$

where $h_\alpha = X_{\alpha+1} - X_\alpha$ and subscripts 'z1' and 'z2' refer to quantities evaluated at the first and second Zannah collocation points: $X_\alpha + (1/2 \mp \sqrt{3}/6)h_\alpha$.

First and second time derivatives of $\Omega$ are given by:

$$ \dot{\Omega}_{\alpha} = \frac{h_{\alpha}}{2} (\dot{\xi}_{\alpha z1} + \dot{\xi}_{\alpha z2}) + \frac{\sqrt{3}h_{\alpha}^2}{12} (\mathrm{ad}_{\xi_{\alpha z1}} \dot{\xi}_{\alpha z2} - \mathrm{ad}_{\xi_{\alpha z2}} \dot{\xi}_{\alpha z1}) \qquad (\mathrm{C2}) $$

$$ \ddot{\Omega}_{\alpha} = \frac{h_{\alpha}}{2} (\ddot{\xi}_{\alpha z1} + \dot{\xi}_{\alpha z2}) + \frac{\sqrt{3}h_{\alpha}^2}{12} (\mathrm{ad}_{\xi_{\alpha z1}} \ddot{\xi}_{\alpha z2} - \mathrm{ad}_{\xi_{\alpha z2}} \ddot{\xi}_{\alpha z1}) \qquad (\mathrm{C3}) $$

$$ + \frac{\sqrt{3}h_{\alpha}^2}{6} \mathrm{ad}_{\xi_{\alpha z1}} \dot{\xi}_{\alpha z2} \qquad (\mathrm{C3}) $$

By virtue of the Magnus expansion, the soft body is computationally identical to $n_p - 1$ rigid joints with equivalent constant joint strain ($\Omega_\alpha$) and its time derivatives ($\dot{\Omega}_\alpha$ and $\ddot{\Omega}_\alpha$). By introducing $q$ from (1) and taking partial derivative of (C1) with respect to $q$ we get,

$$ \frac{\partial \Omega_{\alpha}}{\partial q} = \mathcal{Z}_{\alpha} \qquad (\mathrm{C4}) $$

where $\mathcal{Z}_\alpha \in \mathbb{R}^{6 \times n_{\mathrm{dof},\alpha}}$ is given by:

$$ \mathcal{Z}_{\alpha} = \frac{h_{\alpha}}{2} (\Phi_{\xi z1} + \Phi_{\xi z2}) + \frac{\sqrt{3}h_{\alpha}^2}{12} (\mathrm{ad}_{\xi_{z1}} \Phi_{\xi z2} - \mathrm{ad}_{\xi_{z2}} \Phi_{\xi z1}) \qquad (\mathrm{C5}) $$

Note that for a rigid joint, $\Omega_\alpha = \xi_\alpha$ and $\mathcal{Z}_\alpha = \Phi_{\xi\alpha}$.

Using this we can rewrite (C2) and (C3) as:
$$ \dot{\Omega}_{\alpha} = \mathcal{Z}_{\alpha} \dot{q} \qquad (C6) $$

$$ \ddot{\Omega}_{\alpha} = \mathcal{Z}_{\alpha} \ddot{q} + \dot{Z}_{\alpha} \dot{q} \qquad (C7) $$

where,

$$ \dot{\mathcal{Z}}_{\alpha} = \frac{\sqrt{3}h_{\alpha}^{2}}{12} (\mathrm{ad}_{\xi_{z1}} \Phi_{\xi_{z2}} - \mathrm{ad}_{\xi_{z2}} \Phi_{\xi_{z1}}) \qquad (C8) $$

Using (C6) we derive the derivatives of $\dot{\Omega}_{\alpha}$ with respect to $q$ and $\dot{q}$ as follows:

$$ \frac{\partial \dot{\Omega}_{\alpha}}{\partial q} = \frac{\partial \mathcal{Z}_{\alpha}}{\partial q} \dot{q} = \dot{Z}_{\alpha} \qquad (C9) $$

$$ \frac{\partial \dot{\Omega}_{\alpha}}{\partial \dot{q}} = \mathcal{Z}_{\alpha} \qquad (C10) $$

Similarly, the partial derivatives of $\ddot{\Omega}_{\alpha}$ with respect to $q$, $\dot{q}$, and $\dot{q}$ are given by:

$$ \frac{\partial \ddot{\Omega}_{\alpha}}{\partial q} = \frac{\partial \mathcal{Z}_{\alpha}}{\partial q} \ddot{q} = \ddot{Z}_{\alpha}, \qquad (C11) $$

$$ \frac{\partial \ddot{\Omega}_{\alpha}}{\partial \dot{q}} = \frac{\partial \dot{Z}_{\alpha}}{\partial \dot{q}} \dot{q} + \dot{Z}_{\alpha} = 2 \dot{Z}_{\alpha} \qquad (C12) $$

$$ \frac{\partial \ddot{\Omega}_{\alpha}}{\partial \dot{q}} = \mathcal{Z}_{\alpha} \qquad (C13) $$

where,

$$ \ddot{Z}_{\alpha} = \frac{\sqrt{3}h_{\alpha}^{2}}{12} (\mathrm{ad}_{\xi_{z1}} \Phi_{\xi_{z2}} - \mathrm{ad}_{\xi_{z2}} \Phi_{\xi_{z1}}) \qquad (C14) $$

## Joint Motion Subspace and its Derivatives

The joint motion subspace of the virtual rigid joint $\alpha$ is given by:

$$ S_{\alpha} = T_{\alpha} \mathcal{Z}_{\alpha} \qquad (C15) $$

where $T_{\alpha}$ is the tangent operator (Appendix A).

The relative velocity of the joint expressed in the frame at $X = 0$ is given by $S_{\alpha} \dot{q}$. The partial derivative of this term with respect to $q$ is given by:

$$ \begin{aligned} \frac{\partial S_{\alpha}}{\partial \dot{q}} \dot{q} &= \frac{\partial T_{\alpha}}{\partial \dot{q}} \mathcal{Z}_{\alpha} \dot{q} + T_{\alpha} \frac{\partial \mathcal{Z}_{\alpha}}{\partial \dot{q}} \dot{q} \\ &= \frac{\partial T_{\alpha}}{\partial q} \dot{\Omega}_{\alpha} + T_{\alpha} \dot{Z}_{\alpha} \end{aligned} \qquad (C16) $$

To derive partial derivatives of $T_{\alpha}$, we need to compute $\frac{\partial \dot{\theta}_{\alpha}}{\partial q}$ and $\frac{\partial \dot{\Omega}_{\alpha}}{\partial q} \dot{\Omega}_{\alpha}$. From the definition of $\theta$ (Appendix A) we derive,

$$ \frac{\partial \dot{\theta}_{\alpha}}{\partial q} = \frac{1}{\theta} \dot{\Omega}_{\alpha}^{T} I_{\theta} \mathcal{Z}_{\alpha} \qquad (C17) $$

With appropriate mathematical techniques, we can see that

$$ \frac{\partial \dot{\mathrm{ad}}_{\Omega_{\alpha}} \dot{\Omega}_{\alpha}}{\partial q} = - \sum_{u=1}^{r} \mathrm{ad}_{\Omega_{\alpha}}^{u-1} \mathrm{ad}_{\dot{\Omega}_{\alpha}}^{r-u} \dot{\Omega}_{\alpha} \mathcal{Z}_{\alpha} \qquad (C18) $$

Putting these together in (C16) we get:

$$ \begin{aligned} \frac{\partial S_{\alpha}}{\partial \dot{q}} \dot{q} &= \frac{1}{\theta} \sum_{r=1}^{4} f'_{r}(\theta) \mathrm{ad}_{\Omega_{\alpha}}^{r} \dot{\Omega}_{\alpha} \Omega_{\alpha}^{T} I_{\theta} \mathcal{Z}_{\alpha} \\ &\quad - \sum_{r=1}^{4} f_{r}(\theta) \sum_{u=1}^{r} \mathrm{ad}_{\Omega_{\alpha}}^{u-1} \mathrm{ad}_{\dot{\Omega}_{\alpha}}^{r-u} \dot{\Omega}_{\alpha} \mathcal{Z}_{\alpha} + T_{\alpha} \dot{Z}_{\alpha} \end{aligned} \qquad (C19) $$

Note that due to the order of tensor multiplication, here $\frac{\partial S_{\alpha}}{\partial \dot{q}} \dot{q} \neq \dot{S}_{\alpha}$. In the element notation, the former is given by $\frac{\partial S_{\alpha ij}}{\partial \dot{q}_k} \dot{q}_j$, while the latter is given by $\frac{\partial S_{\alpha ij}}{\partial q_k} \dot{q}_k$.

Similarly, we obtain:

$$ \begin{aligned} \frac{\partial S_{\alpha}}{\partial \dot{q}} \ddot{q} = & \frac{1}{\theta} \sum_{r=1}^{4} f'_{r}(\theta) \mathrm{ad}_{\Omega_{\alpha}}^{r} \mathcal{Z}_{\alpha} \ddot{q} \Omega_{\alpha}^{T} I_{\theta} \mathcal{Z}_{\alpha} \\ & - \sum_{r=1}^{4} f_{r}(\theta) \sum_{u=1}^{r} \mathrm{ad}_{\Omega_{\alpha}}^{u-1} \mathrm{ad}_{\dot{\Omega}_{\alpha}}^{r-u} \dot{\Omega}_{\alpha} \dot{q} \mathcal{Z}_{\alpha} + T_{\alpha} \dot{Z}_{\alpha} \end{aligned} \qquad (C20) $$

Moving to the time derivative of joint motion subspace $\dot{S}_{\alpha}$, it is obtained by taking a time derivative of (C15),

$$ \dot{S}_{\alpha} = \dot{T}_{\alpha} \mathcal{Z}_{\alpha} + T_{\alpha} \dot{Z}_{\alpha} \qquad (C21) $$

where $\dot{T}$ is provided in Appendix A. Then, the derivative of $\dot{S}_{\alpha} \dot{q}$ with respect to $q$ is given by:

$$ \begin{aligned} \frac{\partial \dot{S}_{\alpha}}{\partial \dot{q}} \dot{q} &= \frac{\partial \dot{T}_{\alpha}}{\partial \dot{q}} \mathcal{Z}_{\alpha} \dot{q} + \dot{T}_{\alpha} \frac{\partial \mathcal{Z}_{\alpha}}{\partial \dot{q}} \dot{q} + \frac{\partial T_{\alpha}}{\partial \dot{q}} \dot{Z}_{\alpha} \dot{q} \\ &= \frac{\partial \dot{T}_{\alpha}}{\partial \dot{q}} \dot{\Omega}_{\alpha} + \dot{T}_{\alpha} \dot{Z}_{\alpha} + \frac{\partial T_{\alpha}}{\partial \dot{q}} \dot{Z}_{\alpha} \dot{q} \end{aligned} \qquad (C22) $$

To evaluate the partial derivative of $\dot{T}_{\alpha}$, we need to compute $\frac{\partial \dot{\theta}_{\alpha}}{\partial q}$. Using the definition of $\dot{\theta}$ in Appendix A, we get

$$ \frac{\partial \dot{\theta}_{\alpha}}{\partial q} = - \frac{\dot{\theta}}{\theta^2} \Omega_{\alpha}^{T} I_{\theta} \mathcal{Z}_{\alpha} + \frac{1}{\theta} \dot{\Omega}_{\alpha}^{T} I_{\theta} \mathcal{Z}_{\alpha} + \frac{1}{\theta} \Omega_{\alpha}^{T} I_{\theta} \dot{Z}_{\alpha} \qquad (C23) $$

Substituting (C17), (C18), and (C23) in (C22) we get:

$$ \begin{aligned} \frac{\partial \dot{S}_{\alpha}}{\partial \dot{q}} \dot{q} &= \frac{1}{\theta} \sum_{r=1}^{4} \left( f''_{r}(\theta) \dot{\theta} \mathrm{ad}_{\Omega}^{r} + f'_{r}(\theta) (\mathrm{ad}_{\Omega}^{r}) \dot{\Omega}_{\alpha} \right) \dot{\Omega}_{\alpha} \Omega_{\alpha}^{T} I_{\theta} \mathcal{Z}_{\alpha} \\ &\quad + \frac{1}{\theta} \sum_{r=1}^{4} f'_{r}(\theta) \mathrm{ad}_{\Omega}^{r} \dot{\Omega}_{\alpha} \left( (\dot{\Omega}_{\alpha}^{T} - \frac{\dot{\theta}}{\theta} \Omega_{\alpha}^{T}) I_{\theta} \mathcal{Z}_{\alpha} + \Omega_{\alpha}^{T} I_{\theta} \dot{Z}_{\alpha} \right) \\ &\quad - \sum_{r=1}^{4} f'_{r}(\theta) \dot{\theta} \sum_{u=1}^{r} \mathrm{ad}_{\Omega_{\alpha}}^{u-1} \mathrm{ad}_{\dot{\Omega}_{\alpha}}^{r-u} \dot{\Omega}_{\alpha} \mathcal{Z}_{\alpha} \\ &\quad + \mathrm{ad}_{\Omega_{\alpha}}^{u-1} \mathrm{ad}_{\dot{\Omega}_{\alpha}} \sum_{p=1}^{r-u} \mathrm{ad}_{\Omega_{\alpha}}^{p-1} \mathrm{ad}_{\dot{\Omega}_{\alpha}}^{r-u-p} \dot{\Omega}_{\alpha} \mathcal{Z}_{\alpha} \\ &\quad + \mathrm{ad}_{\Omega_{\alpha}}^{u-1} \mathrm{ad}_{\dot{\Omega}_{\alpha}}^{r-u} \dot{\Omega}_{\alpha} \dot{Z}_{\alpha} \\ &\quad + \frac{1}{\theta} \sum_{r=1}^{4} f'_{r}(\theta) \mathrm{ad}_{\Omega}^{r} \dot{Z}_{\alpha} \dot{q} \Omega_{\alpha}^{T} I_{\theta} \mathcal{Z}_{\alpha} \\ &\quad - \sum_{r=1}^{4} f_{r}(\theta) \sum_{u=1}^{r} \mathrm{ad}_{\Omega_{\alpha}}^{u-1} \mathrm{ad}_{\dot{\Omega}_{\alpha}}^{r-u} \dot{\Omega}_{\alpha} q \mathcal{Z}_{\alpha} \\ &\quad + \dot{T}_{\alpha} \dot{Z}_{\alpha} \end{aligned} \qquad (C24) $$

where,

$$ \begin{aligned} f''_{1}(\theta) &= \frac{1}{2\theta^4} (24 - (24 - 6\theta^2) \cos(\theta) - (18\theta - \theta^3) \sin(\theta)) \\ f''_{2}(\theta) &= \frac{1}{2\theta^5} (24\theta - (60 - 9\theta^2) \sin(\theta) + (36\theta - \theta^3) \cos(\theta)) \\ f''_{3}(\theta) &= \frac{1}{2\theta^6} (40 - (40 - 8\theta^2) \cos(\theta) - (28\theta - \theta^3) \sin(\theta)) \\ f''_{4}(\theta) &= \frac{1}{2\theta^7} (40\theta - (90 - 11\theta^2) \sin(\theta) + (50\theta - \theta^3) \cos(\theta)) \end{aligned} $$
Also, noticing that $\frac{\partial \dot{\theta}_{\alpha}}{\partial \dot{q}} = \frac{\partial \theta_{\alpha}}{\partial q}$, and $\frac{\partial \dot{\Omega}_{\alpha}}{\partial \dot{q}} = \frac{\partial \Omega_{\alpha}}{\partial q}$, it can be seen that,

$$ \frac{\partial \dot{S}_{\alpha}}{\partial \dot{q}} \dot{q} = \frac{\partial S_{\alpha}}{\partial q} \dot{q} \qquad (C25) $$

Finally, the first term in the partial derivative of $ID$ with respect to $q$ in (21) is $\frac{\partial S_{\alpha}^{T}}{\partial q} \mathcal{F}^{C}$. From (C15) we get,

$$ \frac{\partial S_{\alpha}^{T}}{\partial q} \mathcal{F}^{C} = \frac{\partial \mathcal{Z}_{\alpha}^{T}}{\partial q} T_{\alpha}^{T} \mathcal{F}^{C} + \mathcal{Z}_{\alpha}^{T} \frac{\partial T_{\alpha}^{T}}{\partial q} \mathcal{F}^{C} \qquad (C26) $$

where,

$$ \mathcal{Z}_{\alpha}^{T} = \frac{h_{\alpha}}{2} \left( \Phi_{\xi_{z1}}^{T} + \Phi_{\xi_{z2}}^{T} \right) + \frac{\sqrt{3}h_{\alpha}^{2}}{12} \left( \Phi_{\xi_{z1}}^{T} \mathrm{ad}_{\xi_{z2}}^{*} - \Phi_{\xi_{z2}}^{T} \mathrm{ad}_{\xi_{z1}}^{*} \right) \qquad (C27) $$

and

$$ T^{T}(\Omega) = I_{6} + \sum_{r=1}^{4} (-1)^{r} f_{r}(\theta) \mathrm{ad}_{\Omega}^{*r} \qquad (C28) $$

Based on previously derived results and the identity (B3) we obtain:

$$ \frac{\partial S_{\alpha}^{T}}{\partial q} \mathcal{F}^{C} = \frac{\sqrt{3}h_{\alpha}^{2}}{12} \left( \Phi_{\xi_{z1}}^{T} \overline{\mathrm{ad}}_{T_{\alpha}^{T} \mathcal{F}^{C}} \Phi_{\xi_{z2}} - \Phi_{\xi_{z2}}^{T} \overline{\mathrm{ad}}_{T_{\alpha}^{T} \mathcal{F}^{C}} \Phi_{\xi_{z1}} \right) \\ + \frac{1}{\theta} \mathcal{Z}_{\alpha}^{T} \sum_{r=1}^{4} (-1)^{r} f_{r}^{\prime}(\theta) \mathrm{ad}_{\Omega_{\alpha}}^{*r} \mathcal{F}^{C} \Omega_{\alpha}^{T} I_{\theta} \mathcal{Z}_{\alpha} \\ + \mathcal{Z}_{\alpha}^{T} \sum_{r=1}^{4} (-1)^{r} f_{r}(\theta) \sum_{u=1}^{r} \mathrm{ad}_{\Omega_{\alpha}}^{*u-1} \overline{\mathrm{ad}}_{\mathrm{ad}_{\Omega_{\alpha}}^{*r-u} \mathcal{F}^{C}} \mathcal{Z}_{\alpha} \qquad (C29) $$

## Derivatives of Adjoint Maps

Using identities (B10) and (B4) we get,

$$ \frac{\partial \mathrm{Ad}_{\exp(\hat{\Omega}_{\alpha})}}{\partial q_p} \mathcal{V} = -\mathrm{Ad}_{\mathrm{Ad}_{\exp(\hat{\Omega}_{\alpha})}} \mathcal{V} \mathcal{S}_{\alpha, p} \qquad (C30) $$

where $q_p$ is the $p^{th}$ component of $q$ and $S_{\alpha,p}$ is the $p^{th}$ column of $S_{\alpha}$. Combining all the columns, we get:

$$ \frac{\partial \mathrm{Ad}_{\exp(\hat{\Omega}_{\alpha})}}{\partial q} \mathcal{V} = -\mathrm{Ad}_{\mathrm{Ad}_{\exp(\hat{\Omega}_{\alpha})}} \mathcal{V} \mathcal{S}_{\alpha} \qquad (C31) $$

Similarly using identities (B11), (B12), and (B13) we derive:

$$ \frac{\partial \mathrm{Ad}^*_{\exp(\hat{\Omega}_{\alpha})}}{\partial q} \mathcal{F} = \overline{\mathrm{Ad}}_{\mathrm{Ad}^*_{\exp(\hat{\Omega}_{\alpha})}} \mathcal{F} \mathcal{S}_{\alpha} \qquad (C32) $$

$$ \frac{\partial \mathrm{Ad}^{-1}_{\exp(\hat{\Omega}_{\alpha})}}{\partial q} \mathcal{V} = \mathrm{Ad}^{-1}_{\exp(\hat{\Omega}_{\alpha})} \mathrm{ad} \mathcal{V} \mathcal{S}_{\alpha} \qquad (C33) $$

$$ \frac{\partial \mathrm{Ad}^{-}_{\exp(\hat{\Omega}_{\alpha})}}{\partial q} \mathcal{F} = -\mathrm{Ad}^{-}_{\exp(\hat{\Omega}_{\alpha})} \overline{\mathrm{ad}}^{*} \mathcal{F} \mathcal{S}_{\alpha} \qquad (C34) $$

# D Section-wise Derivations

## Derivations in 4.1

Inverse dynamics of the soft body, enabled by the “Pseudo rigid joint” formulation of GVS is given by:

$$ ID = \sum_{\alpha=1}^{n_p-1} ID_{\alpha} = \sum_{\alpha=1}^{n_p-1} S_{\alpha}^{T} \mathcal{F}_{\alpha}^{C} \qquad (D1) $$

where $\mathcal{F}_{\alpha}^{C}$ is computed recursively according to:

$$ \mathcal{F}_{\alpha}^{C} = \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}}^{*} \mathcal{F}_{k} = \mathrm{Ad}_{g_{\alpha \alpha+1}}^{*} (\mathcal{F}_{\alpha+1} + \mathcal{F}_{\alpha+1}^{C}) \qquad (D2) $$

The partial derivative of $ID_{\alpha}$ with respect to $q$ is given by:

$$ \frac{\partial ID_{\alpha}}{\partial q} = \frac{\partial S_{\alpha}^{T}}{\partial q} \mathcal{F}_{\alpha}^{C} + S_{\alpha}^{T} \sum_{k=\alpha+1}^{n_p} \frac{\partial \mathrm{Ad}_{g_{\alpha k}}^{*}}{\partial q} \mathcal{F}_{k} \\ + S_{\alpha}^{T} \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}}^{*} \frac{\partial \mathcal{F}_{k}}{\partial q} \\ = \frac{\partial S_{\alpha}^{T}}{\partial q} \mathcal{F}_{\alpha}^{C} + S_{\alpha}^{T} \left( \mathcal{N}_{\alpha}^{C} R_{\alpha}^{B} + \mathcal{M}_{\alpha}^{C} Q_{\alpha}^{B} + U_{\alpha}^{S} + P_{\alpha}^{S} \right) \qquad (D3) $$

The analytical derivative of the first term in (D3) is given by (C29). The term $P_{\alpha}^{S}$ arises from the partial derivative of the coAdjoint maps. The rest of the terms come from the partial derivatives of $\mathcal{F}_{k}$ transformed to the frame of $\alpha$.

## Derivation of $P_{\alpha}^{S}$ in (D3)

Using the identity, $\mathrm{Ad}_{g_{\alpha k}}^{*} = \prod_{\beta=\alpha}^{k-1} \mathrm{Ad}_{g_{\beta \beta+1}}^{*}$ and (C32) we get:

$$ \frac{\partial \mathrm{Ad}_{g_{\alpha k}}^{*}}{\partial q_p} \mathcal{F}_{k} = \sum_{\beta=\alpha}^{k-1} \mathrm{Ad}_{g_{\alpha \beta}}^{*} \overline{\mathrm{Ad}}_{\mathrm{Ad}_{g_{\beta k}}^{*}} \mathcal{F}_{k} S_{\beta, p} \qquad (D4) $$

where $q_p$ is the $p^{th}$ component of $q$ and $S_{\alpha,p}$ is the $p^{th}$ column of $S_{\alpha}$. Combining all the columns, we get:

$$ \frac{\partial \mathrm{Ad}_{g_{\alpha k}}^{*}}{\partial q} \mathcal{F}_{k} = \sum_{\beta=\alpha}^{k-1} \mathrm{Ad}_{g_{\alpha \beta}}^{*} \overline{\mathrm{Ad}}_{\mathrm{Ad}_{g_{\beta k}}^{*}} \mathcal{F}_{k} S_{\beta} \qquad (D5) $$

Using this we derive,

$$ \sum_{k=\alpha+1}^{n_p} \frac{\partial \mathrm{Ad}_{g_{\alpha k}}^{*}}{\partial q} \mathcal{F}_{k} = \sum_{k=\alpha+1}^{n_p} \sum_{\beta=\alpha}^{k-1} \mathrm{Ad}_{g_{\alpha \beta}}^{*} \overline{\mathrm{Ad}}_{\mathrm{Ad}_{g_{\beta k}}^{*}} \mathcal{F}_{k} S_{\beta} \qquad (D6) $$

Using summation identity (B18), and identities (B4) and (B5) we get:

$$ \sum_{k=\alpha+1}^{n_p} \frac{\partial \mathrm{Ad}_{g_{\alpha k}}^{*}}{\partial q} \mathcal{F}_{k} = \sum_{\beta=\alpha}^{n_p-1} \mathrm{Ad}_{g_{\alpha \beta}}^{*} \overline{\mathrm{Ad}}_{\mathcal{F}_{\beta}^{C}} \mathcal{S}_{\beta} = P_{\alpha}^{S} \qquad (D7) $$

Recursive computation of $P_{\alpha}^{S}$ can be derived from (D7):

$$ P_{\alpha}^{S} = \overline{\mathrm{Ad}}_{\mathcal{F}_{\alpha}^{C}} S_{\alpha} + \mathrm{Ad}_{g_{\alpha \alpha+1}}^{*} P_{\alpha+1}^{S} \qquad (D8) $$

## Derivation of $\mathcal{N}_{\alpha}^{C} R_{\alpha}^{B} + \mathcal{M}_{\alpha}^{C} Q_{\alpha}^{B} + U_{\alpha}^{S}$ in (D3)

The resultant point wrench in the local frame is given by,

$$ \mathcal{F}_{k} = \mathcal{M}_{k} \dot{\eta}_{k} + \mathrm{ad}^{*}_{\eta_{k}} \mathcal{M}_{k} \eta_{k} - \mathcal{M}_{k} \mathrm{Ad}_{g_{k}}^{-1} \mathcal{G} \qquad (D9) $$

Using, $\dot{\eta}_{k} = \gamma_{k} + \varphi_{k}$ (inertial and Coriolis components of acceleration), we rewrite $\mathcal{F}_{k}$ as the resultant of inertial, Coriolis, and gravitational forces.

$$ \begin{aligned} \mathcal{F}_{k} &= (\mathcal{M}_{k} \gamma_{k}) + (\mathrm{ad}^{*}_{\eta_{k}} \mathcal{M}_{k} \eta_{k} + \mathcal{M}_{k} \varphi_{k}) - (\mathcal{M}_{k} \mathrm{Ad}_{g_{k}}^{-1} \mathcal{G}) \\ &= \mathcal{F}_{Ik} + \mathcal{F}_{Ck} + \mathcal{F}_{Gk} \end{aligned} \qquad (D10) $$
The partial derivatives of $\mathcal{F}_k$ involves the partial derivatives of $\mathcal{F}_{Ik}$, $\mathcal{F}_{Ck}$, and $\mathcal{F}_{Gk}$.

## Contribution of Gravity

We have,

$$ \mathcal{F}_{Gk} = -\mathcal{M}_k \mathrm{Ad}_{g_k}^{-1} \mathcal{G} \quad (\mathrm{D}11) $$

Using $\mathrm{Ad}_{g_k}^{-1} = \prod_{\beta=k-1}^{1} \mathrm{Ad}_{g_{\beta\beta+1}}^{-1}$ and (C33) we get:

$$ \frac{\partial \mathcal{F}_{Gk}}{\partial \mathbf{q}} = -\mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathrm{ad}_{\mathrm{Ad}_{g_{\beta}}^{-1} \mathbf{g}} \mathcal{S}_{\beta} \quad (\mathrm{D}12) $$

## Contribution of Coriolis Force

The Coriolis component of $\mathcal{F}_k$ is given by:

$$ \mathcal{F}_{Ck} = \mathrm{ad}_{\eta_k}^* \mathcal{M}_k \eta_k + \mathcal{M}_k \varphi_k \quad (\mathrm{D}13) $$

where,

$$ \eta_k = \sum_{l=1}^{k-1} \mathrm{Ad}_{g_{lk}}^{-1} \mathbf{S}_l \dot{\mathbf{q}} \quad (\mathrm{D}14) $$

and

$$ \varphi_k = \sum_{l=1}^{k-1} \mathrm{Ad}_{g_{lk}}^{-1} \left( \mathrm{ad}_{\eta_l} \mathbf{S}_l \dot{\mathbf{q}} + \dot{\mathbf{S}}_l \dot{\mathbf{q}} \right) \quad (\mathrm{D}15) $$

Using (B3) we write,

$$ \frac{\partial \mathcal{F}_{Ck}}{\partial \mathbf{q}} = \left( \overline{\mathrm{ad}}_{\mathcal{M}_k}^* \eta_k + \mathrm{ad}_{\eta_k}^* \mathcal{M}_k \right) \frac{\partial \eta_k}{\partial \mathbf{q}} + \mathcal{M}_k \frac{\partial \varphi_k}{\partial \mathbf{q}} \quad (\mathrm{D}16) $$

To compute this, we need to derive the partial derivatives of $\eta_k$ and $\varphi_k$. From (D14) we get,

$$ \frac{\partial \eta_k}{\partial \mathbf{q}} = \sum_{l=1}^{k-1} \frac{\partial \mathrm{Ad}_{g_{lk}}^{-1}}{\partial \mathbf{q}} \mathbf{S}_l \dot{\mathbf{q}} + \mathrm{Ad}_{g_{lk}}^{-1} \frac{\partial \mathbf{S}_l}{\partial \mathbf{q}} \dot{\mathbf{q}} \quad (\mathrm{D}17) $$

The second term is given by (C19), while the first can be derived similarly to (D12). We get,

$$ \begin{aligned} \sum_{l=1}^{k-1} \frac{\partial \mathrm{Ad}_{g_{lk}}^{-1}}{\partial \mathbf{q}} \mathbf{S}_l \dot{\mathbf{q}} &= \sum_{l=1}^{k-1} \sum_{\beta=l}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathrm{ad}_{\mathrm{Ad}_{g_{l\beta}}^{-1} \mathbf{S}_l \dot{\mathbf{q}}} \mathcal{S}_{\beta} \\ &= \sum_{l=1}^{k-1} \sum_{\beta=l}^{k-1} \mathrm{ad}_{\mathrm{Ad}_{g_{lk}}^{-1} \mathbf{S}_l \dot{\mathbf{q}}} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathcal{S}_{\beta} \end{aligned} \quad (\mathrm{D}18) $$

Based on the summation identity (B16), (B5) we can rewrite this as:

$$ \sum_{l=1}^{k-1} \frac{\partial \mathrm{Ad}_{g_{lk}}^{-1}}{\partial \mathbf{q}} \mathbf{S}_l \dot{\mathbf{q}} = \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathrm{ad}_{\eta_{\beta}^{+}} \mathcal{S}_{\beta} \quad (\mathrm{D}19) $$

where, $\eta_{\beta}^{+} = \sum_{l=1}^{\beta} \mathrm{Ad}_{g_{l\beta}}^{-1} \mathbf{S}_l \dot{\mathbf{q}} = \eta_{\beta} + \mathcal{S}_{\beta} \dot{\mathbf{q}}$.

Substituting this into (D17) we get:

$$ \frac{\partial \eta_k}{\partial \mathbf{q}} = \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{R}_{\beta} \quad (\mathrm{D}20) $$

where,

$$ \mathbf{R}_{\beta} = \mathrm{ad}_{\eta_{\beta}^{+}} \mathcal{S}_{\beta} + \frac{\partial \mathcal{S}_{\beta}}{\partial \mathbf{q}} \dot{\mathbf{q}} \quad (\mathrm{D}21) $$

From (D15) we get,

$$ \frac{\partial \varphi_k}{\partial \mathbf{q}} = \sum_{l=1}^{k-1} \left( \frac{\partial \mathrm{Ad}_{g_{lk}}^{-1}}{\partial \mathbf{q}} \left( \mathrm{ad}_{\eta_l} \mathbf{S}_l \dot{\mathbf{q}} + \dot{\mathbf{S}}_l \dot{\mathbf{q}} \right) - \mathrm{Ad}_{g_{lk}}^{-1} \mathrm{ad}_{\mathrm{S}_l \dot{\mathbf{q}}} \frac{\partial \eta_l}{\partial \mathbf{q}} \right. \\ \left. + \mathrm{Ad}_{g_{lk}}^{-1} \left( \mathrm{ad}_{\eta_l} \frac{\partial \mathbf{S}_l}{\partial \mathbf{q}} \dot{\mathbf{q}} + \frac{\partial \dot{\mathbf{S}}_l}{\partial \mathbf{q}} \dot{\mathbf{q}} \right) \right) \quad (\mathrm{D}22) $$

The first term has the same form of (D19). Hence, we get:

$$ \sum_{l=1}^{k-1} \frac{\partial \mathrm{Ad}_{g_{lk}}^{-1}}{\partial \mathbf{q}} \left( \mathrm{ad}_{\eta_l} \mathbf{S}_l \dot{\mathbf{q}} + \dot{\mathbf{S}}_l \dot{\mathbf{q}} \right) = \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathrm{ad}_{\varphi_{\beta}^{+}} \mathcal{S}_{\beta} \quad (\mathrm{D}23) $$

where, $\varphi_{\beta}^{+} = \varphi_{\beta} + \mathrm{ad}_{\eta_{\beta}} \mathcal{S}_{\beta} \dot{\mathbf{q}} + \dot{\mathbf{S}}_{\beta} \dot{\mathbf{q}}$.

The second term of (D22) can be expanded using equation (D20) and the identities (B4) and (B17) as follows:

$$ \begin{aligned} - \sum_{l=1}^{k-1} \mathrm{Ad}_{g_{lk}}^{-1} \mathrm{ad}_{S_l \dot{\mathbf{q}}} \frac{\partial \eta_l}{\partial \mathbf{q}} &= - \sum_{l=1}^{k-1} \sum_{\beta=1}^{l-1} \mathrm{ad}_{\mathrm{Ad}_{g_{lk}}^{-1} S_l \dot{\mathbf{q}}} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{R}_{\beta} \\ &= - \sum_{\beta=1}^{k-2} \left( \sum_{l=\beta+1}^{k-1} \mathrm{ad}_{\mathrm{Ad}_{g_{lk}}^{-1} S_l \dot{\mathbf{q}}} \right) \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{R}_{\beta} \\ &= - \sum_{\beta=1}^{k-2} \mathrm{ad}_{\eta_k - \mathrm{Ad}_{g_{\beta k}}^{-1} \eta_{\beta}^{+}} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{R}_{\beta} \\ &= \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathrm{ad}_{\eta_{\beta}^{+} - \mathrm{Ad}_{g_{\beta k}} \eta_k} \mathbf{R}_{\beta} \end{aligned} \quad (\mathrm{D}24) $$

Note that we used $\eta_{k-1}^{+} = \mathrm{Ad}_{g_{k-1k}} \eta_k$ for changing the summation limit.

The analytical formulas of the last of (D22) term are provided by (C19) and (C24). Substituting (D23) and (D24) in (D22) we get:

$$ \frac{\partial \varphi_k}{\partial \mathbf{q}} = \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \left( \mathrm{ad}_{\varphi_{\beta}^{+}} \mathcal{S}_{\beta} + \mathrm{ad}_{\eta_{\beta}^{+} - \mathrm{Ad}_{g_{\beta k}} \eta_k} \mathbf{R}_{\beta} \right. \\ \left. + \mathrm{ad}_{\eta_{\beta}} \frac{\partial \mathcal{S}_{\beta}}{\partial \mathbf{q}} \dot{\mathbf{q}} + \frac{\partial \dot{\mathcal{S}}_{\beta}}{\partial \mathbf{q}} \dot{\mathbf{q}} \right) \quad (\mathrm{D}25) $$

Substituting (D25) and (D20) into (D16), we get the final form of the derivative of Coriolis force:

$$ \frac{\partial \mathcal{F}_{Ck}}{\partial \mathbf{q}} = \left( \overline{\mathrm{ad}}_{\mathcal{M}_k}^* \eta_k + \mathrm{ad}_{\eta_k}^* \mathcal{M}_k \right) \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{R}_{\beta} \\ + \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \left( \mathrm{ad}_{\varphi_{\beta}^{+}} \mathcal{S}_{\beta} + \mathrm{ad}_{\eta_{\beta}^{+} - \mathrm{Ad}_{g_{\beta k}} \eta_k} \mathbf{R}_{\beta} \right. \\ \left. + \mathrm{ad}_{\eta_{\beta}} \frac{\partial \mathcal{S}_{\beta}}{\partial \mathbf{q}} \dot{\mathbf{q}} + \frac{\partial \dot{\mathcal{S}}_{\beta}}{\partial \mathbf{q}} \dot{\mathbf{q}} \right) \quad (\mathrm{D}26) $$

## Contribution of Inertial Force

The inertial component of $\mathcal{F}_k$ is given by:

$$ \mathcal{F}_{Ik} = \mathcal{M}_k \gamma_k = \mathcal{M}_k \sum_{l=1}^{k-1} \mathrm{Ad}_{g_{lk}}^{-1} \mathbf{S}_l \ddot{\mathbf{q}} \quad (\mathrm{D}27) $$

It is easy to see that the derivative of this term follows the same form of (D20). We get:

$$ \frac{\partial \mathcal{F}_{Ik}}{\partial \mathbf{q}} = \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \left( \mathrm{ad}_{\gamma_{\beta}^{+}} \mathcal{S}_{\beta} + \frac{\partial \mathcal{S}_{\beta}}{\partial \mathbf{q}} \ddot{\mathbf{q}} \right) \quad (\mathrm{D}28) $$

where, $\gamma_{\beta}^{+} = \sum_{l=1}^{\beta} \mathrm{Ad}_{g_{l\beta}}^{-1} \mathbf{S}_l \ddot{\mathbf{q}} = \gamma_{\beta} + \mathcal{S}_{\beta} \ddot{\mathbf{q}}$.
Combining all Contributions

Combining partial derivatives of gravitational force (D12), Coriolis force (D26), and inertial force (D28), we get:

$$ \frac{\partial \mathcal{F}_k}{\partial \dot{q}} = \mathcal{N}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{R}_{\beta} + \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}}^{-1} \mathbf{Q}_{\beta} \quad (\text{D29}) $$

where,

$$ \mathcal{N}_k = \overline{\mathrm{ad}}_{\mathcal{M}_k \eta_k}^* + \mathrm{ad}_{\eta_k}^* \mathcal{M}_k - \mathcal{M}_k \mathrm{ad}_{\eta_k} \quad (\text{D30}) $$

$$ \begin{aligned} \mathbf{Q}_{\beta} = & \overline{\mathrm{ad}}_{\eta_{\beta}^{+}} \mathbf{S}_{\beta} + \mathrm{ad}_{\eta_{\beta}^{+}} \mathbf{R}_{\beta} \\ & + \mathrm{ad}_{\eta_{\beta}} \frac{\partial \mathbf{S}_{\beta}}{\partial \dot{q}} \dot{\mathbf{q}} + \frac{\partial \dot{\mathbf{S}}_{\beta}}{\partial \dot{q}} \dot{\dot{\mathbf{q}}} + \frac{\partial \mathbf{S}_{\beta}}{\partial \dot{q}} \ddot{\mathbf{q}} - \mathrm{ad}_{\overline{\mathrm{ad}}_{g_{\beta}^{-1}} \mathbf{g}} \mathbf{S}_{\beta} \end{aligned} \quad (\text{D31}) $$

In the backward pass, the partial derivatives of $\mathcal{F}_k$ are transformed to the frames of each “virtual joint” using coAdjoint maps. Using (D29), $\mathrm{Ad}_{g_{\beta k}}^{-1} = \mathrm{Ad}_{g_{\alpha k}}^{-1} \mathrm{Ad}_{g_{\alpha \beta}}$, and summation identity (B19) we get:

$$ \begin{aligned} \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}^*} \frac{\partial \mathcal{F}_k}{\partial \dot{q}} &= \sum_{k=\alpha+1}^{n_p} \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\alpha k}}^* \left( \mathcal{N}_k \mathrm{Ad}_{g_{\alpha k}^{-1}} \mathrm{Ad}_{g_{\alpha \beta}} \mathbf{R}_{\beta} \right. \\ & \qquad \left. + \mathcal{M}_k \mathrm{Ad}_{g_{\alpha k}^{-1}} \mathrm{Ad}_{g_{\alpha \beta}} \mathbf{Q}_{\beta} \right) \\ &= \mathcal{N}_{\alpha}^{C} \sum_{\beta=1}^{\alpha-1} \mathrm{Ad}_{g_{\alpha \beta}} \mathbf{R}_{\beta} + \mathcal{M}_{\alpha}^{C} \sum_{\beta=1}^{\alpha-1} \mathrm{Ad}_{g_{\alpha \beta}} \mathbf{Q}_{\beta} \\ & \qquad + \sum_{\beta=\alpha}^{n_p-1} \mathrm{Ad}_{g_{\alpha \beta}^*} \left( \mathcal{N}_{\beta}^{C} \mathbf{R}_{\beta} + \mathcal{M}_{\beta}^{C} \mathbf{Q}_{\beta} \right) \end{aligned} \quad (\text{D32}) $$

We can simplify (D32) into:

$$ \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}^*} \frac{\partial \mathcal{F}_k}{\partial \dot{q}} = \mathcal{N}_{\alpha}^C \mathbf{R}_{\alpha}^B + \mathcal{M}_{\alpha}^C \mathbf{Q}_{\alpha}^B + U_{\alpha}^S \quad (\text{D33}) $$

where,

$$ \mathbf{R}_{\alpha}^{B} = \sum_{\beta=1}^{\alpha-1} \mathrm{Ad}_{g_{\alpha \beta}} \mathbf{R}_{\beta} \quad (\text{D34}) $$

$$ \mathbf{Q}_{\alpha}^{B} = \sum_{\beta=1}^{\alpha-1} \mathrm{Ad}_{g_{\alpha \beta}} \mathbf{Q}_{\beta} \quad (\text{D35}) $$

$$ \mathcal{N}_{\alpha}^{C} = \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}^{*}} \mathcal{N}_{k} \mathrm{Ad}_{g_{\alpha k}^{-1}} \quad (\text{D36}) $$

$$ \mathcal{M}_{\alpha}^{C} = \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}^{*}} \mathcal{M}_{k} \mathrm{Ad}_{g_{\alpha k}^{-1}} \quad (\text{D37}) $$

$$ U_{\alpha}^{S} = \sum_{\beta=\alpha}^{n_p-1} \mathrm{Ad}_{g_{\alpha \beta}^{*}} \left( \mathcal{N}_{\beta}^{C} \mathbf{R}_{\beta} + \mathcal{M}_{\beta}^{C} \mathbf{Q}_{\beta} \right) \quad (\text{D38}) $$

$\mathbf{R}_{\alpha}^{B}$ and $\mathbf{Q}_{\alpha}^{B}$ can be recursively computed during the forward pass, while $\mathcal{N}_{\alpha}^{C}$, $\mathcal{M}_{\alpha}^{C}$, and $U_{\alpha}^{S}$ can be computed during the backward pass as follows:

$$ \mathbf{R}_{\alpha}^{B} = \mathrm{Ad}_{g_{\alpha-1}^{-1}}^{-1} (\mathbf{R}_{\alpha-1} + \mathbf{R}_{\alpha-1}^{B}) \quad (\text{D39}) $$

$$ \mathbf{Q}_{\alpha}^{B} = \mathrm{Ad}_{g_{\alpha-1}^{-1}}^{-1} (\mathbf{Q}_{\alpha-1} + \mathbf{Q}_{\alpha-1}^{B}) \quad (\text{D40}) $$

$$ \mathcal{N}_{\alpha}^{C} = \mathrm{Ad}_{g_{\alpha+1}^{*}}^{*} (\mathcal{N}_{\alpha+1} + \mathcal{N}_{\alpha+1}^{C}) \mathrm{Ad}_{g_{\alpha+1}^{-1}} \quad (\text{D41}) $$

$$ \mathcal{M}_{\alpha}^{C} = \mathrm{Ad}_{g_{\alpha+1}^{*}}^{*} (\mathcal{M}_{\alpha+1} + \mathcal{M}_{\alpha+1}^{C}) \mathrm{Ad}_{g_{\alpha+1}^{-1}} \quad (\text{D42}) $$

$$ U_{\alpha}^{S} = \mathcal{N}_{\alpha}^{C} \mathbf{R}_{\alpha} + \mathcal{M}_{\alpha}^{C} \mathbf{Q}_{\alpha} + \mathrm{Ad}_{g_{\alpha \alpha+1}^{*}}^{*} U_{\alpha+1}^{S} \quad (\text{D43}) $$

Derivations in 4.2

The partial derivative of $ID_{\alpha}$ with respect to $\dot{q}$ is given by:

$$ \frac{\partial ID_{\alpha}}{\partial \dot{q}} = S_{\alpha}^{T} (\mathcal{N}_{\alpha}^{C} S_{\alpha}^{B} + \mathcal{M}_{\alpha}^{C} Y_{\alpha}^{B} + V_{\alpha}^{S}) \quad (\text{D44}) $$

Only the Coriolis force is a function of $\dot{q}$. Therefore, $\mathcal{N}_{\alpha}^{C} S_{\alpha}^{B} + \mathcal{M}_{\alpha}^{C} Y_{\alpha}^{B} + V_{\alpha}^{S}$ in (D44) originates from the partial derivative of $\mathcal{F}_{Ck}$ with respect to $\dot{q}$ transformed to the frame of $\alpha$. From (D13) we get:

$$ \frac{\partial \mathcal{F}_k}{\partial \dot{q}} = \left( \overline{\mathrm{ad}}_{\mathcal{M}_k \eta_k}^* + \mathrm{ad}_{\eta_k}^* \mathcal{M}_k \right) \frac{\partial \eta_k}{\partial \dot{q}} + \mathcal{M}_k \frac{\partial \varphi_k}{\partial \dot{q}} \quad (\text{D45}) $$

From (D14) we get:

$$ \frac{\partial \eta_k}{\partial \dot{q}} = \sum_{l=1}^{k-1} \mathrm{Ad}_{g_{lk}^{-1}} S_l \quad (\text{D46}) $$

From (D15) we get:

$$ \frac{\partial \varphi_k}{\partial \dot{q}} = \sum_{l=1}^{k-1} \mathrm{Ad}_{g_{lk}^{-1}} \left( -\mathrm{ad} s_l \dot{q} \frac{\partial \eta_l}{\partial \dot{q}} + \mathrm{ad} \eta_l s_l + \frac{\partial \dot{s}_l}{\partial \dot{q}} \dot{q} + \dot{s}_l \right) \quad (\text{D47}) $$

Using similar mathematical operations like that of (D24) we get:

$$ \frac{\partial \varphi_k}{\partial \dot{q}} = \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}^{-1}} \left( \mathrm{ad}_{\eta_{\beta}^{+} - \mathrm{Ad}_{g_{\beta k}^{-1}} \eta_k} \mathbf{S}_{\beta} + \mathrm{ad}_{\eta_{\beta}} \mathbf{S}_{\beta} \right. \\ \left. + \frac{\partial \mathbf{S}_{\beta}}{\partial \dot{q}} \dot{q} + \dot{S}_{\beta} \right) \quad (\text{D48}) $$

Substituting (D48) into (D45) we get:

$$ \begin{aligned} \frac{\partial \mathcal{F}_k}{\partial \dot{q}} = & \left( \overline{\mathrm{ad}}_{\mathcal{M}_k \eta_k}^* + \mathrm{ad}_{\eta_k}^* \mathcal{M}_k \right) \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}^{-1}} \mathbf{S}_{\beta} \\ & + \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}^{-1}} \left( \mathrm{ad}_{\eta_{\beta}^{+} - \mathrm{Ad}_{g_{\beta k}^{-1}} \eta_k} \mathbf{S}_{\beta} + \mathrm{ad}_{\eta_{\beta}} \mathbf{S}_{\beta} \right) \quad (\text{D49}) \\ & + \frac{\partial \mathbf{S}_{\beta}}{\partial \dot{q}} \dot{q} + \dot{S}_{\beta} \right) $$

By rearranging this, we get:

$$ \frac{\partial \mathcal{F}_k}{\partial \dot{q}} = \mathcal{N}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}^{-1}} \mathbf{S}_{\beta} + \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{g_{\beta k}^{-1}} \mathbf{Y}_{\beta} \quad (\text{D50}) $$

where,

$$ \mathbf{Y}_{\beta} = \mathrm{ad}_{\eta_{\beta}^{+}} \mathbf{S}_{\beta} + \mathrm{ad}_{\eta_{\beta}} \mathbf{S}_{\beta} + \frac{\partial \mathbf{S}_{\beta}}{\partial \dot{q}} \dot{q} + \dot{S}_{\beta} \quad (\text{D51}) $$
Transforming the partial derivatives of $\mathcal{F}_k$ at all points $k > \alpha$ to the frame of $\alpha$, similar to (D32), we get:

$$ \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}}^* \frac{\partial \mathcal{F}_k}{\partial \dot{q}} = \mathcal{N}_{\alpha}^C \mathcal{S}_{\alpha}^B + \mathcal{M}_{\alpha}^C \mathcal{Y}_{\alpha}^B + \mathcal{V}_{\alpha}^S \quad (\text{D52}) $$

$\mathcal{S}_{\alpha}^B$ and $\mathcal{Y}_{\alpha}^B$ can be recursively computed during the forward pass, while $\mathcal{V}_{\alpha}^S$ can be computed during the backward pass as follows:

$$ \mathcal{S}_{\alpha}^{B} = \mathrm{Ad}_{g_{\alpha-1\alpha}}^{-1} (\mathcal{S}_{\alpha-1} + \mathcal{S}_{\alpha-1}^{B}) \quad (\text{D53}) $$

$$ \mathcal{Y}_{\alpha}^{B} = \mathrm{Ad}_{g_{\alpha-1\alpha}}^{-1} (\mathcal{Y}_{\alpha-1} + \mathcal{Y}_{\alpha-1}^{B}) \quad (\text{D54}) $$

$$ \mathcal{V}_{\alpha}^{S} = \mathcal{N}_{\alpha}^{C} \mathcal{S}_{\alpha} + \mathcal{M}_{\alpha}^{C} \mathcal{Y}_{\alpha} + \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*} \mathcal{V}_{\alpha+1}^{S} \quad (\text{D55}) $$

### Derivations in 4.3

The partial derivatives of $ID_{\alpha}$ with respect to $\ddot{q}$ is given by:

$$ \frac{\partial ID_{\alpha}}{\partial \ddot{q}} = \mathcal{S}_{\alpha}^{T} (\mathcal{M}_{\alpha}^{C} \mathcal{S}_{\alpha}^{B} + \mathcal{W}_{\alpha}^{S}) \quad (\text{D56}) $$

Since, only the inertial force is a function of $\ddot{q}$, $\mathcal{M}_{\alpha}^{C} \mathcal{S}_{\alpha}^{B} + \mathcal{W}_{\alpha}^{S}$ in (D56) originates from the partial derivative of $\mathcal{F}_{Ik}$ with respect to $\ddot{q}$ transformed to the frame of $\alpha$. From (D27) we get:

$$ \frac{\partial \mathcal{F}_k}{\partial \ddot{q}} = \mathcal{M}_k \sum_{\beta=1}^{k-1} \mathrm{Ad}_{\ddot{g}_{\beta k}}^{-1} \mathcal{S}_{\beta} \quad (\text{D57}) $$

Transforming the partial derivatives of $\mathcal{F}_k$ at all points $k > \alpha$ to the frame of $\alpha$, similar to (D32), we get:

$$ \sum_{k=\alpha+1}^{n_p} \mathrm{Ad}_{g_{\alpha k}}^{*} \frac{\partial \mathcal{F}_k}{\partial \ddot{q}} = \mathcal{M}_{\alpha}^{C} \mathcal{S}_{\alpha}^{B} + \mathcal{W}_{\alpha}^{S} \quad (\text{D58}) $$

where, $\mathcal{W}_{\alpha}^{S}$ is computed recursively during the backward pass.

$$ \mathcal{W}_{\alpha}^{S} = \mathcal{M}_{\alpha}^{C} \mathcal{S}_{\alpha} + \mathrm{Ad}_{g_{\alpha\alpha+1}}^{*} \mathcal{W}_{\alpha+1}^{S} \quad (\text{D59}) $$

### Derivations in 4.4

The generalized actuation force for tendon-like actuators is given by:

$$ \mathcal{B}u = \sum_{i=1}^{n_p} w_i \phi_{\xi_i}^T \sum_{k=1}^{n_a} \begin{pmatrix} \tilde{d}_{ik} t_{ik} \\ t_{ik} \end{pmatrix} u_k \quad (\text{D60}) $$

where, $d_{ik}(X_i) = [0 \ y_{ik} \ z_{ik}]^T$ is the local coordinates of the $k^{th}$ tendon and $t_{ik}$ is the unit tangent vector. We have,

$$ t_{ik} = \frac{\mathcal{T}_{ik}}{\| \mathcal{T}_{ik} \|} \quad (\text{D61}) $$

$$ \mathcal{T}_{ik} = \hat{\xi}_i \begin{pmatrix} d_{ik} \\ 1 \\ 0 \end{pmatrix} + \begin{pmatrix} d'_{ik} \\ 0 \end{pmatrix} \quad (\text{D62}) $$

Using the equation of strain parameterization (1) and (D61) we get:

$$ \frac{\partial \mathcal{B}}{\partial q} u = \sum_{i=1}^{n_p} w_i \phi_{\xi_i}^T \sum_{k=1}^{n_a} \frac{u_k}{\| \mathcal{T}_{ik} \|} \begin{pmatrix} \tilde{d}_{ik} \tilde{t}_{ik}^2 \tilde{d}_{ik} & -\tilde{d}_{ik} \tilde{t}_{ik}^2 \\ \tilde{t}_{ik}^2 \tilde{d}_{ik} & -\tilde{t}_{ik}^2 \end{pmatrix} \phi_{\xi_i} \quad (\text{D63}) $$

Generalized stiffness and damping matrices of a soft body are given by:

$$ \mathbf{K} = \sum_{i=1}^{n_p} w_i \phi_{\xi_i}^T \mathbf{\Sigma}(X_i) \phi_{\xi_i} \quad (\text{D64}) $$

$$ \mathbf{D} = \sum_{i=1}^{n_p} w_i \phi_{\xi_i}^T \mathbf{\Upsilon}(X_i) \phi_{\xi_i} \quad (\text{D65}) $$

where $\mathbf{\Sigma}(X_i) = \mathrm{diag}(GJ_{x_i}, EJ_{y_i}, EJ_{z_i}, EA_i, GA_i, GA_i) \in \mathbb{R}^{6 \times 6}$ is the screw elasticity matrix, $E$ being the Young modulus and $G$ the shear modulus and $\mathbf{\Upsilon}(X_i) = v\mathrm{diag}(J_{x_i}, 3J_{y_i}, 3J_{z_i}, 3A_i, A_i, A_i) \in \mathbb{R}^{6 \times 6}$ is the screw damping matrix, $v$ being the material damping.

### Derivations in 5.1

The generalized actuation force due to rigid joints is given by:

$$ \begin{aligned} \mathcal{B}u &= \sum_k \left[ \mathbf{0}_{1 \times n_k^-} \left( \mathbf{\Phi}_k^T \mathbf{T}_k^T \mathrm{Ad}_{\exp(\hat{\xi}_k)}^* \mathbf{\Phi}_k \mathbf{u}_k \right)^T \mathbf{0}_{1 \times n_k^+} \right]^T \\ &= \sum_k \left[ \mathbf{0}_{1 \times n_k^-} \left( \mathcal{S}_k^T \mathcal{F}_{k0} \right)^T \mathbf{0}_{1 \times n_k^+} \right]^T \end{aligned} \quad (\text{D66}) $$

where $k$ is the index of actuated rigid joints and $\mathcal{F}_{k0} = \mathrm{Ad}_{\exp(\hat{\xi}_k)}^* \mathbf{\Phi}_k \mathbf{u}_k$.

The partial derivative of (D66) term with respect to $q$ can be obtained using:

$$ \begin{aligned} \frac{\partial \mathcal{S}_k^T \mathcal{F}_{k0}}{\partial q_k} &= \frac{1}{\theta} \mathbf{\Phi}_k^T \sum_{r=1}^{4} (-1)^r f'_r(\theta) \mathrm{Ad}_{\hat{\Omega}_k}^{*r} \mathcal{F}_{k0} \Omega_k^T I_\theta \mathbf{\Phi}_k \\ &\quad + \mathbf{\Phi}_k^T \sum_{r=1}^{4} (-1)^r f_r(\theta) \sum_{u=1}^{r} \mathrm{Ad}_{\hat{\Omega}_k}^{*u-1} \overline{\mathrm{Ad}_{\hat{\Omega}_k}^{*r-u}} \mathcal{F}_{k0} \mathbf{\Phi}_k \\ &\quad + \mathcal{S}_k^T \overline{\mathrm{Ad}_{\mathcal{F}_{k0}}^*} \mathcal{S}_k \end{aligned} \quad (\text{D67}) $$

Note that, for 1 DoF joints, $T = \mathrm{Ad}_{\exp(\hat{\xi}_k)}$. Hence, for such joints, $B$ is a constant, and its derivative is $0$.

### Derivations in 5.2

#### Jacobian of Constraint Wrench

Using identity (B13) and (B3) we derive partial derivative of the constraint force acting on $A$ with respect to $q_p$ as follows:

$$ \begin{aligned} \frac{\partial \mathcal{F}_{kA}}{\partial q_p} &= \frac{\partial \mathrm{Ad}_{g_{kBA}}^{-*}}{\partial q_p} \mathbf{\Phi}_{\perp k} \lambda_k \\ &= -\mathrm{Ad}_{J_{kBA,p}}^{*} \mathrm{Ad}_{g_{kBA}}^{-} \mathbf{\Phi}_{\perp k} \lambda_k \\ &= -\overline{\mathrm{Ad}_{\mathrm{Ad}_{g_{kBA}}^{*-} \mathbf{\Phi}_{\perp k} \lambda_k}} (\mathbf{J}_{kA,p} - \mathrm{Ad}_{g_{kBA}}^{-1} \mathbf{J}_{kB,p}) \end{aligned} \quad (\text{D68}) $$

Collecting all columns and using (B9), we get:

$$ \begin{aligned} \frac{\partial \mathcal{F}_{kA}}{\partial q} &= -\overline{\mathrm{Ad}_{\mathrm{Ad}_{g_{kBA}}^{*-} \mathbf{\Phi}_{\perp k} \lambda_k}} (\mathbf{J}_{kA} - \mathrm{Ad}_{g_{kBA}}^{-1} \mathbf{J}_{kB}) \\ &= -\mathrm{Ad}_{g_{kBA}}^{*-} \overline{\mathrm{Ad}_{\mathbf{\Phi}_{\perp k} \lambda_k}} (\mathrm{Ad}_{g_{kBA}} \mathbf{J}_{kA} - \mathbf{J}_{kB}) \end{aligned} \quad (\text{D69}) $$

Constraint Jacobian with respect to $q$

The closed-chain constraint equation is given by:
$$C(\ddot{q}, \dot{q}, q) = A\ddot{q} + \dot{A}\dot{q} + \frac{2}{T_B} A\dot{q} + \frac{1}{T_B^2} e = 0 \quad (\text{D70})$$

The Jacobian of this equation with respect to $q$ is given by:

$$\frac{\partial C}{\partial q} = \frac{\partial A}{\partial q} \ddot{q} + \frac{\partial \dot{A}}{\partial q} \dot{q} + \frac{2}{T_B} \frac{\partial A}{\partial q} \dot{q} + \frac{1}{T_B^2} \frac{\partial e}{\partial q} \quad (\text{D71})$$

From the definition of the Pfaffian matrix and its derivative in equations (52a) and (52b), we can rewrite the first three terms of equation (D71) as follows:

$$\frac{\partial A}{\partial q} \ddot{q} + \frac{\partial \dot{A}}{\partial q} \dot{q} + \frac{2}{T_B} \frac{\partial A}{\partial q} \dot{q} = \left\|_{k=1}^{n_{CL}} \Phi_{\perp k}^{T} \left( \begin{aligned} & \frac{\partial \mathrm{Ad}_{g_{kBA}} \dot{\eta}_{kA}}{\partial q} \\ & - \frac{\partial \dot{\eta}_{kB}}{\partial q} + \frac{\partial \mathrm{Ad}_{g_{kBA}} \mathrm{ad}_{\eta_{kBA}} \eta_{kA}}{\partial q} \\ & + \frac{2}{T_B} \left( \frac{\partial \mathrm{Ad}_{g_{kBA}} \eta_{kA}}{\partial q} - \frac{\partial \eta_{kB}}{\partial q} \right) \end{aligned} \right) \right\| \quad (\text{D72})$$

Equation (D20) provides the Jacobian of velocity twist. From (D28) and (D25) we can infer the expression for the Jacobian of the acceleration twist.

$$\frac{\partial \eta_k}{\partial q} = R_k^B \quad (\text{D73})$$

$$\frac{\partial \dot{\eta}_k}{\partial q} = L_k^B - \mathrm{ad}_{\eta_k} R_k^B \quad (\text{D74})$$

where $L_k^B$ is $Q_k^B$ (26b) without the gravity component.

Using (B10), the derivative of the Adjoint terms in (D72) can be derived as:

$$\frac{\partial \mathrm{Ad}_{g_{kBA}} \dot{\eta}_{kA}}{\partial q} = - \mathrm{ad}_{\mathrm{Ad}_{g_{kBA}}} \dot{\eta}_{kA} (\mathrm{Ad}_{g_{kBA}} \mathbf{J}_{kA} - \mathbf{J}_{kB}) \quad (\text{D75})$$

$$\frac{\partial \mathrm{Ad}_{g_{kBA}} \eta_{kA}}{\partial q} = - \mathrm{ad}_{\mathrm{Ad}_{g_{kBA}}} \eta_{kA} (\mathrm{Ad}_{g_{kBA}} \mathbf{J}_{kA} - \mathbf{J}_{kB}) \quad (\text{D76})$$

The final form of (D72) can be obtained by combining (D73), (D74), (D75), and (D76).

The last term in (D71) includes the Jacobian of the closed-chain kinematic error. Let $\epsilon = (\log(g_{kBA}))^\vee$. We have,

$$\exp(\hat{\epsilon}) = g_{kBA} \quad (\text{D77})$$

Using identities (B14) and (B15) we get:

$$\exp(\hat{\epsilon}) \left( \mathrm{Ad}_{\exp(\hat{\epsilon})} T(\epsilon) \frac{\partial \epsilon}{\partial q_p} \right) = g_{kBA} \hat{J}_{kBA,p} \quad (\text{D78})$$

where $T$ is the tangent operator provided in Appendix A.

Rearranging the terms and considering all columns we get:

$$\frac{\partial \epsilon}{\partial q} = T^{-1} (\log(g_{kBA})^\vee) (\mathrm{Ad}_{g_{kBA}} \mathbf{J}_{kA} - \mathbf{J}_{kB}) \quad (\text{D79})$$

Using this we obtain the Jacobian of kinematic constraints.

$$\frac{\partial e}{\partial q} = \left\|_{k=1}^{n_{CL}} \Phi_{\perp k}^{T} T^{-1} (\log(g_{kBA})^{\vee}) (\mathrm{Ad}_{g_{kBA}} \mathbf{J}_{kA} - \mathbf{J}_{kB}) \right\| \quad (\text{D80})$$

## Constraint Jacobian with respect to $\dot{q}$

The Jacobian of the equation (D70) with respect to $\dot{q}$ is given by:

$$\frac{\partial C}{\partial \dot{q}} = \frac{\partial \dot{A}}{\partial \dot{q}} \dot{q} + \dot{A} + \frac{2}{T_B} A \quad (\text{D81})$$

From the time derivative of the Pfaffian matrix, we can see that:

$$\frac{\partial \dot{A}}{\partial \dot{q}} \dot{q} + \dot{A} = \left\|_{k=1}^{n_{CL}} \Phi_{\perp k} \left( \begin{aligned} & \frac{\partial \mathrm{Ad}_{g_{kBA}} \varphi_{kA}}{\partial \dot{q}} - \frac{\partial \varphi_{kB}}{\partial \dot{q}} \\ & + \frac{\partial \mathrm{Ad}_{g_{kBA}} \mathrm{ad}_{\eta_{kBA}} \eta_{kA}}{\partial \dot{q}} \end{aligned} \right) \right\| \quad (\text{D82})$$

where, the partial derivative of $\varphi_k$ is given by equation (D48). We can rewrite it as:

$$\frac{\partial \varphi_k}{\partial \dot{q}} = Y_k^B - \mathrm{ad}_{\eta_k} S_k^B \quad (\text{D83})$$

## Derivations in 6.2

The contact force in the global frame is given by:

$$f_c = k_c \delta_k^p u_{\perp k} \quad (\text{D84})$$

The partial derivative of $f_c$ with respect to $q$ is given by:

$$\frac{\partial f_c}{\partial q} = k_c u_{\perp k} p \delta_k^{p-1} \frac{\partial \delta}{\partial q} + k_c \delta_k^p \frac{\partial u_{\perp k}}{\partial q} \quad (\text{D85})$$

From the definition of penetration (63) and unit normal, we get:

$$\frac{\partial \delta}{\partial q} = u_{\perp k}^{T} C_{1} \frac{\partial r_{gk}}{\partial q} \quad (\text{D86})$$

and

$$\begin{aligned} \frac{\partial u_{\perp k}}{\partial q} &= \frac{1}{\|n_{\perp k}\|} C_1 \frac{\partial r_{gk}}{\partial q} - n_{\perp k} \frac{1}{\|n_{\perp k}\|^3} n_{\perp k}^T C_1 \frac{\partial r_{gk}}{\partial q} \\ &= \frac{1}{\|n_{\perp k}\|} (I_3 - u_{\perp k} u_{\perp k}^T) C_1 \frac{\partial r_{gk}}{\partial q} \end{aligned} \quad (\text{D87})$$

Linear velocity in the global frame is given by $v_{gk} = \frac{\partial r_{gk}}{\partial q} \dot{q} = [0 \ I_3] \mathrm{Ad}_{g_{\theta k}} S_k^B \dot{q}$. Using this we get,

$$\frac{\partial r_{gk}}{\partial q} = [0 \ I_3] \mathrm{Ad}_{g_{\theta k}} S_k^B \quad (\text{D88})$$

Substituting (D86), (D87), and (D88) in (D85) and using $F_{gk} = [0^{1 \times 3} f_c^T]^T$ we get,

$$\frac{\partial F_{gk}}{\partial q} = \begin{pmatrix} 0 & 0 \\ 0 & k^* C_1 \end{pmatrix} \mathrm{Ad}_{g_{\theta k}} S_k^B \quad (\text{D89})$$

where, $k^* = k_c p \delta_k^{p-1} u_{\perp k} u_{\perp k}^T + \frac{k_c \delta_k^p}{\|n_{\perp k}\|} (I_3 - u_{\perp k} u_{\perp k}^T)$.

## Derivations in 6.3

The drag-lift force is given by:

$$\mathcal{F}_{Dk} = \mathcal{D}_k \|v_k\| \eta_k \quad (\text{D90})$$

The partial derivative of $\mathcal{F}_{Dk}$ with respect to $q$ is,

$$\frac{\partial \mathcal{F}_{Dk}}{\partial q} = \mathcal{D}_k \eta_k \frac{\partial \|v_k\|}{\partial q} + \mathcal{D}_k \|v_k\| \frac{\partial \eta_k}{\partial q} \quad (\text{D91})$$

Using $\|v_k\| = \sqrt{\eta_k^T I_v \eta_k}$ we get,

$$\frac{\partial \|v_k\|}{\partial q} = \frac{1}{\|v_k\|} \eta_k^T I_v \frac{\partial \eta_k}{\partial q} \quad (\text{D92})$$
Substituting (D92) in (D91) we get,

$$ \frac{\partial \mathcal{F}_{Dk}}{\partial \dot{q}} = \mathcal{D}_k^* \frac{\partial \eta_k}{\partial \dot{q}} \qquad (D93) $$

where, $\mathcal{D}_k^* = \mathcal{D}_k (\frac{1}{\|v_k\|} \eta_k \eta_k^T I_v + \|v_k\| I_6)$.

It can be seen that the partial derivative of $\mathcal{F}_{Dk}$ with respect to $\dot{q}$ takes a similar form.

$$ \frac{\partial \mathcal{F}_{Dk}}{\partial \dot{q}} = \mathcal{D}_k^* \frac{\partial \eta_k}{\partial \dot{q}} \qquad (D94) $$

## References

Alkayas, A. Y., Mathew, A. T., Feliu-Talegon, D., Deng, P., Thuruthel, T. G. and Renda, F. (2025), ‘Soft synergies: Model order reduction of hybrid soft-rigid robots via optimal strain parameterization’, *IEEE Transactions on Robotics* **41**, 1118–1137.

Armanini, C., Boyer, F., Mathew, A. T., Duriez, C. and Renda, F. (2023), ‘Soft Robots Modeling: A Structured Overview’, *IEEE Transactions on Robotics* **PP**, 1–21.

Armanini, C., Farman, M., Calisti, M., Giorgio-Serchi, F., Stefanini, C. and Renda, F. (2022a), ‘Flagellate underwater robotics at macroscale: Design, modeling, and characterization’, *IEEE Transactions on Robotics* **38**(2), 731–747.

Armanini, C., Hussain, I., Iqbal, M. Z., Gan, D., Prattichizzo, D. and Renda, F. (2021), ‘Discrete cosserat approach for closed-chain soft robots: Application to the fin-ray finger’, *IEEE Transactions on Robotics* **37**(6), 2083–2098.

Bergou, M., Wardetzky, M., Robinson, S., Audoly, B. and Grinspun, E. (2008), ‘Discrete elastic rods’, *SIGGRAPH’08: International Conference on Computer Graphics and Interactive Techniques, ACM SIGGRAPH 2008 Papers* **2008** **27**(3).

Bjelonic, M., Grandia, R., Geilinger, M., Harley, O., Medeiros, V. S., Pajovic, V., Jelavic, E., Coros, S. and Hutter, M. (2022), ‘Offline motion libraries and online mpc for advanced mobility skills’, *International Journal of Robotics Research* **41**, 903–924.

Boyer, F., Gotelli, A., Tempel, P., Lebastard, V., Renda, F. and Briot, S. (2023), ‘Implicit time-integration simulation of robots with rigid bodies and cosserat rods based on a newton-euler recursive algorithm’, *Trans. Rob.* **40**, 677–696.
URL: https://doi.org/10.1109/TRO.2023.3334647

Boyer, F., Lebastard, V., Candelier, F. and Renda, F. (2020), ‘Dynamics of continuum and soft robots: A strain parameterization based approach’, *IEEE Transactions on Robotics* pp. 1–17.

Boyer, F., Lebastard, V., Candelier, F. and Renda, F. (2022), ‘Statistics and Dynamics of Continuum Robots Based on Cosserat Rods and Optimal Control Theories’, *IEEE Transaction on Robotics* **PP**, 1–19.

Boyer, F., Porez, M. and Khalil, W. (2006), ‘Macro-continuous computed torque algorithm for a three-dimensional eel-like robot’, *Robotics, IEEE Transactions on* **22**(4), 763–775.

Burgess, J. (1992), ‘Bending stiffness in a simulation of undersea cable deployment’.
URL: https://dx.doi.org/

Butcher, J. (2016), *Numerical Differential Equation Methods, John Wiley & Sons, Ltd, chapter 2, pp. 55–142*.

Bächer, M., Knoop, E. and Schumacher, C. (2021), ‘Design and control of soft robots using differentiable simulation’, *Current Robotics Reports* **2**, 211–221.

Cao, D. Q. and Tucker, R. W. (2008), ‘Nonlinear dynamics of elastic rods using the Cosserat theory: Modelling and simulation’, *International Journal of Solids and Structures* **45**(2), 460–477.

Cardona, A. and Géradin, M. (1994), *Numerical Integration of Second Order Differential—Algebraic Systems in Flexible Mechanism Dynamics*, Springer Netherlands, Dordrecht, pp. 501–529.
URL: https://doi.org/10.1007/978-94-011-1166-9_16

Carpentier, J. and Mansard, N. (2018), Analytical derivatives of rigid body dynamics algorithms, in ‘Robotics: Science and systems (RSS 2018)’.

Carpentier, J., Saurel, G., Buondonno, G., Mirabel, J., Lamiraux, F., Stasse, O. and Mansard, N. (2019), The pinocchio c++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives, in ‘2019 IEEE/SICE International Symposium on System Integration (SII)’, pp. 614–619.

Choi, A., Jing, R., Sabelhaus, A. P. and Jawed, M. K. (2024), ‘Dismech: A discrete differential geometry-based physical simulator for soft robots and structures’, *IEEE Robotics and Automation Letters* **9**, 3483–3490.

Chung, J. and Hulbert, G. M. (1993), ‘A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized-α Method’, *Journal of Applied Mechanics* **60**(2), 371–375.

Du, T., Wu, K., Ma, P., Wah, S., Spielberg, A., Rus, D. and Matusik, W. (2022), ‘Diffpd: Differentiable projective dynamics’, *ACM Transactions on Graphics* **41**.

Duriez, C. (2013), Control of elastic soft robots based on real-time finite element method, in ‘Robotics and Automation (ICRA), 2013 IEEE International Conference on’, pp. 3982–3987.

Eugster, S. R. and Harsch, J. (2023), ‘A family of total Lagrangian Petrov-Galerkin Cosserat rod finite element formulations’, *GAMM Mitteilungen* **46**(2), 1–21.

Featherstone, R. (2008), *Rigid Body Dynamics Algorithms*, Springer New York.

Feliu-Talegon, D., Mathew, A. T., Alkayas, A. Y., Adamu, Y. A. and Renda, F. (2024), ‘Dynamic shape estimation of tendon-driven soft manipulators via actuation readings’, *IEEE Robotics and Automation Letters* pp. 1–8.
URL: https://ieeexplore.ieee.org/document/10777508/

Flores, P. (2022), ‘Contact mechanics for dynamical systems: a comprehensive review’, *Multibody System Dynamics* **54**, 127–177.
URL: http://dx.doi.org/10.1007/s11044-021-09803-y

Gazzola, M., Dudte, L. H., McCormick, A. G. and Mahadevan, L. (2018), ‘Forward and inverse problems in the mechanics of soft filaments’, *Royal Society Open Science* **5**(6), 171628.

Geilinger, M., Hahn, D., Zehnder, J., Bächer, M., Thomaszewski, B. and Coros, S. (2020), ‘Add: Analytically differentiable dynamics for multi-body systems with frictional contact’, *ACM Transactions on Graphics* **39**.

Gifthaler, M., Neunert, M., Stäuble, M., Frigerio, M., Semini, C. and Buchli, J. (2017a), ‘Automatic differentiation of rigid body dynamics for optimal control and estimation’, *Advanced Robotics* **31**, 1225–1237.
URL: https://doi.org/10.1080/01691864.2017.1395361
Giftthaler, M., Neunert, M., Stäuble, M., Frigerio, M., Semini, C. and Buchli, J. (2017b), 'Automatic differentiation of rigid body dynamics for optimal control and estimation', Advanced Robotics 31, 1225-1237.

Hafner, C., Schumacher, C., Knoop, E., Auzinger, T., Bickel, B. and Bächer, M. (2019), 'X-cad: Optimizing cad models with extended finite elements', ACM Transactions on Graphics 38.

Hoshyari, S., Xu, H., Knoop, E., Coros, S. and Bächer, M. (2019), 'Vibration-minimizing motion retargeting for robotic characters', ACM Transactions on Graphics 38.

Howell, T. A., Cleac'h, S. L., Brüdigam, J., Kolter, J. Z., Schwager, M. and Manchester, Z. (2023), 'Dojo: A differentiable physics engine for robotics'. URL: https://arxiv.org/abs/2203.00806

Hu, Y., Fang, Y., Ge, Z., Qu, Z., Zhu, Y., Pradhana, A. and Jiang, C. (2018), 'A moving least squares material point method with displacement discontinuity and two-way rigid body coupling', ACM Transactions on Graphics 37.

Huang, Z., Hu, Y., Du, T., Zhou, S., Su, H., Tenenbaum, J. B. and Gan, C. (2021), 'Plasticinelab: A soft-body manipulation benchmark with differentiable physics', ArXiv abs/2104.03311.
URL: https://api.semanticscholar.org/CorpusID:233169183

Ishigaki, T., Ayusawa, K. and Yamamoto, K. (2024), 'Compre-
hensive gradient computation framework of pcs model for
soft robot simulation', IEEE Robotics and Automation Letters
9, 5990-5997.

Jatavallabhula, K. M., Macklin, M., Golemo, F., Voleti, V., Petrini, L., Weiss, M., Considine, B., Parent-Lévesque, J., Xie, K., Erleben, K., Paull, L., Shkurti, F., Nowrouzezahrai, D., Fidler, S., Robotics, M. and Lab, E. A. (2021), Gradsim: Differentiable simulation for system identification and visuomotor control, in 'International Conference on Learning Representations'. URL: https://gradsim.github.io

Kane, C., Marsden, J. E., Ortiz, M. and West, M. (2000), 'Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems', International Journal for Numerical Methods in Engineering 49, 1295-1325.

Kane, T. R. and Levinson, D. A. (1983), 'The Use of Kane's Dynamical Equations in Robotics', The International Journal of Robotics Research 2(3), 3-21.
URL: https://doi.org/10.1177/027836498300200301

Lidec, Q. L., Jallet, W., Montaut, L., Laptev, I., Schmid, C. and Carpentier, J. (2024), 'Contact models in robotics: A comparative analysis', IEEE Transactions on Robotics 40, 3716-3733.

Lighthill, M. J. (1970), 'Aquatic animal propulsion of high hydromechanical efficiency', Journal of Fluid Mechanics 44, 265-301.

Mathew, A. T., Ben Hmida, I., Caradonna, D. and Peringal, A. B. (2025), 'Sorosim: A matlab toolbox for hybrid soft-rigid robots'. URL: https://github.com/SoRoSim/SoRoSim

Mathew, A. T., Feliu-Talegon, D., Alkayas, A. Y., Boyer, F. and Renda, F. (2024), ‘Reduced order modeling of hybrid soft-rigid robots using global, local, and state-dependent strain parameterization’, Original Article The International Journal of Robotics Research 0, 1–26.

Mathew, A. T., Hmida, I. M. B., Armanini, C., Boyer, F. and Renda, F. (2022a), 'SoRoSim: A MATLAB Toolbox for Hybrid

Rigid-Soft Robots Based on the Geometric Variable-Strain Approach', IEEE Robotics and Automation Magazine.

Meier, C., Popp, A. and Wall, W. A. (2017), ‘Geometrically Exact Finite Element Formulations for Slender Beams: Kirchhoff-Love Theory Versus Simo-Reissner Theory’, Archives of Computational Methods in Engineering 26(1), 163–243.

Morison, J., Johnson, J. W. and Schaaf, S. A. (1950), 'The force exerted by surface waves on piles', Journal of Petroleum Technology 2, 149–154.

Murray, R., Li, Z. and Sastry, S. (1994), *A Mathematical Introduction to Robotic Manipulation*, Taylor & Francis, Boca Raton, USA.

Neunert, M., Stauble, M., Giftthaler, M., Bellicoso, C. D., Carius, J., Gehring, C., Hutter, M. and Buchli, J. (2018), ‘Whole-body nonlinear model predictive control through contacts for quadrupeds’, IEEE Robotics and Automation Letters 3, 1458–1465.

Newbury, R., Collins, J., He, K., Pan, J., Posner, I., Howard, D. and Cosgun, A. (2024), ‘A review of differentiable simulators’. URL: http://arxiv.org/abs/2407.05560

Renda, F., Armanini, C., Lebastard, V., Candelier, F. and Boyer, F. (2020), 'A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation', *IEEE Robotics and Automation Letters* 5(3), 4006–4013.

Renda, F., Mathew, A. and Talegon, D. F. (2024), ‘Dynamics and control of soft robots with implicit strain parametrization’, IEEE Robotics and Automation Letters 9, 2782–2789.

Rucker, D. and Webster, R. (2011), 'Statics and dynamics of continuum robots with general tendon routing and external loading', Robotics, IEEE Transactions on 27(6), 1033–1044.

Simo, J. and Vu-Quoc, L. (1988), 'On the dynamics in space of rods undergoing large motions — a geometrically exact approach', Computer Methods in Applied Mechanics and Engineering 66(2), 125–161.

Singh, S., Russell, R. P. and Wensing, P. M. (2022), ‘Efficient analytical derivatives of rigid-body dynamics using spatial vector algebra’, IEEE Robotics and Automation Letters 7, 1776–1783.

Spielberg, A., Du, T., Hu, Y., Rus, D. and Matusik, W. (2023), 'Advanced soft robot modeling in chainqueen', Robotica 41, 74-104.

Sukhija, B., Kohler, N., Zamora, M., Zimmermann, S., Curi, S., Krause, A. and Coros, S. (2023), Gradient-based trajectory optimization with learned dynamics, in 'Proceedings - IEEE International Conference on Robotics and Automation', Vol. 2023-May, Institute of Electrical and Electronics Engineers Inc., pp. 1011–1018.

Tedrake, R. and the Drake Development Team (2019), ‘Drake: Model-based design and verification for robotics’. URL: https://drake.mit.edu

Till, J., Aloi, V. and Rucker, C. (2019), ‘Real-time dynamics of soft and continuum robots based on cosserat rod models’, The International Journal of Robotics Research 38(6), 723–746.

Tjavaras, A. A., Zhu, Q., Liu, Y., Triantafyllou, M. S. and Yue, D. K. P. (1998), ‘The mechanics of highly-extensible cables’, Journal of Sound and Vibration 213, 709–737.

Todorov, E., Erez, T. and Tassa, Y. (2012), Mujoco: A physics engine for model-based control, in 'IEEE International
Conference on Intelligent Robots and Systems', pp. 5026-5033.

Wensing, P. M., Posa, M., Hu, Y., Escande, A., Mansard, N. and Prete, A. D. (2024), 'Optimization-based control for dynamic legged robots', *IEEE Transactions on Robotics* **40**, 43-63.

Zanna, A. (1999), 'Collocation and relaxed collocation for the fer and the magnus expansions', *SIAM J. Numer. Anal.* **36**(4), 1145-1182.