# Towards Real-Time Simulation of Soft Robots with Contacts using a Method of Hybrid Hyper-Reduction

Olivier Goury, Samuel Mohsen Youssef, Simon Le Berre, Christian Duriez

► To cite this version:

Olivier Goury, Samuel Mohsen Youssef, Simon Le Berre, Christian Duriez. Towards Real-Time Simulation of Soft Robots with Contacts using a Method of Hybrid Hyper-Reduction. ROBOSOFT 2025 - 8th IEEE-RAS International Conference on Soft Robotics, Cosimo della Santina; Josie Hughes, Apr 2025, Lausanne, Switzerland. pp.1-6, 10.1109/RoboSoft63089.2025.11020929. hal-04935584

HAL Id: hal-04935584

https://inria.hal.science/hal-04935584v1

Submitted on 7 Feb 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Towards Real-Time Simulation of Soft Robots with Contacts using a Method of Hybrid Hyper-Reduction

Olivier Goury<fnref n="1" />, Samuel M. Youssef<fnref n="2" />, Simon Le Berre<fnref n="1" /> and Christian Duriez<fnref n="2" />

**Abstract**—Soft robotics has emerged as an important part of robotics in recent years. Soft robots have an inherent view of contacts that is dramatically different from traditional rigid robots. Indeed, for rigid robots, contacts are either forbidden to avoid damage to the robot, the environment and humans, or precisely controlled for locomotion or interaction with an object. For soft robots, contacts may happen without damage, and when interacting with an object, local deformations allows for smoother interactions and potentially better performance. These prospects make soft robots attractive for tasks such as grasping. Fast finite element simulation is very useful for control and design. However, simulating collision adds a major numerical cost as it requires first a collision detection algorithm to detect collisions, and most importantly, it requires solving a constrained problem to avoid inter-penetrations and compute contact forces. When the number of contact points is large, this computation slows down the simulation dramatically. In this paper, we apply a hybrid hyper-reduction method to alleviate the FEM cost, the collision detection as well as the contact response computation. The deformations are computed in a low-dimensional subspace computed from offline experiments. The mechanical matrices are reduced through a method of hyper-reduction and the collision model is reduced following a hybrid reduction strategy. We show good agreement between original and reduced simulation while speeding up dramatically the computation. We first apply the method in simulation on a soft bouncing ball to explain the method. We then show an example with a soft gripper. The method is generic and can be used for control, design or learning algorithms.

## I. INTRODUCTION

For the past 10 to 15 years, soft robotics has made great progress within the field of robotics in general. The prospects are better manoeuvrability in constrained environment, safer interactions with humans, performing well at manipulation tasks and less expensive building costs. However, the accurate control of such robots proves to be more challenging than rigid robots, due to their greater amount of degrees of freedom. This is especially true when considering the interaction of the soft robot with the environment and contact handling.

### A. Soft Robot Simulation

Simulation plays an important role in robotics in general, for virtual prototyping, testing or applying control strategies, motion planning or can also be used for training learning procedures such as reinforcement learning. Several platforms

already exists to simulate soft robots dynamics. We can cite VoxCAD [1] which represent soft robots material in Voxels supported by a lattice structure mechanically linked by beams. This model can be evaluated at fast rate and supports multi-material design. Another platform is Chainqueen [2], [3]. It is based on the material point method (MPM) and is fully differentiable, making it nicely applicable for optimisation purposes. We can also refer to SOFA [4]. SOFA is based on the finite element method and can simulate any geometry of soft robots in interactions with its environment, with different kind of actuation (cable, pneumatic, hydraulic, etc...). Computational burden of FEM can be alleviated using GPU acceleration [5] or model order reduction [6].

### B. Contact Handling in Soft Robotics

The problem of treating collisions in the context of soft robotics simulation has already been looked at in the literature. There are two main ways to treat contacts:

* (1) by regularising the contact constraint, applying a penalty coefficient or defining an elastic energy as objects get too close to each other. This is numerically efficient but has the drawback of introducing a contact parameter that is artificial and has to be tuned depending on the application.

* (2) by imposing strict constraints of non inter-penetrations between the robot and its environment and/or himself. This is the most accurate way in principle, but it involves solving an optimisation problem at each timestep (typically a linear complementarity problem (LCP)). Solving this problem can suffer from poor convergence in complex cases.

In the Chainqueen platform [2], collisions are handled directly through the grid of the MPM formulation. Though fast, the accuracy can be questioned since the contact behaviour will be dependent on the grid discretisation. This formulation belongs to the category (1) with an elasticity directly dependent on the stiffness of the material. In [7], in the context of grasping applications, both self-collision and interaction with obstacles are considered. The simulation is geometry-based and the contact law is regularised using springs.

In [8], a contact model including friction is developed for grasping using rigid fingers. The triangulated contact surface is approximated by a quadratic form, but the contact is solved for using complementarity conditions. In [9], contact is also treated using non-smooth mechanics with a complementarity condition. No specific strategy is developed to reduce the contact space that is defined on a triangulated mesh. That is sufficient for coarse discretisation of the surface. In [10], a

<footnote n="1">¹Olivier Goury and Simon Le Berre are member of the team TRIPOP, Université de Grenoble Alpes, Inria, UMR 5224 LJK, F-38330 Montbonnot, France olivier.goury@inria.fr</footnote>

<footnote n="2">²Samuel M. Youssef and Christian Duriez are member of the team DEFROST, Univ. Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000 Lille, France</footnote>
special care is taken regarding accuracy and computational cost in the case of self-collisions. The dual space of Lagrange multipliers is projected onto a reduced non-negative space with positive coefficients. This allows to keep a similar definition of the contact law with complementarity between reduced contact force and self-collision distance. However, the formulation is dependent on the building of the reduced positive space for Lagrange multipliers, which is hard to compute in practice. Also, the method was not able to account for friction.

### C. Contributions

In this paper, we apply a method of hybrid hyper-reduction to soft robot simulation involving contacts. The contributions are:

* The first application of a hybrid hyper-reduction method in soft robots simulation with contact and friction

* We show the generality of the method by using different constitutive laws including multi-material.

* We provide an automatic framework for speeding up FE simulation of soft robots with contacts with both speedup on the computation of the robot deformation as well as contact forces.

The rest of the paper is organised as follows: in section II, the robot dynamics and contact model is explained. In section III, we recall the model order reduction approach applied for the reduction of the primal space and present the hybrid hyper-reduction method which applies to both primal and dual space. In section IV, we apply the method on two cases and show the evolution of the accuracy as well as the computational gain achieved by the method. Conclusions and discussions are drawn in the final section.

## II. ROBOT MOTION AND CONTACT EQUATION

### A. Robot dynamics

Consider a soft robot discretised on a finite element mesh. Its dynamics can be described using the second law of Newton:

$$M(\mathbf{q})\dot{\mathbf{v}} = \mathbf{P}(t) - \mathbb{F}(\mathbf{q}, \mathbf{v}) + \mathbf{H}^T \mathbf{\lambda}, \quad (1)$$

where $\mathbf{q}$ is the vector of positions, $\mathbf{v}$ is the vector of velocities, $M(\mathbf{q})$ is the mass matrix, $\mathbf{P}$ gathers external forces, $\mathbb{F}$ accounts for internal forces (depending on the constitutive law of the material) and $\mathbf{H}^T \mathbf{\lambda}$ is the vector of constraint forces contributions, which gathers contact constraints with friction or actuation constraints.

1) *Time discretization:* Let us now discretise the continuous time space in a discrete time interval $[t_0, t_1, ..., t_{n_t}]$, with $h = t_{n+1} - t_n$.

Using an implicit Euler integration scheme, setting $\mathbf{v} = \dot{\mathbf{q}}$, we obtain:

$$M(\mathbf{v}_{t_{n+1}} - \mathbf{v}_{t_n}) + h\mathbf{F}(\mathbf{q}_{t_{n+1}}, \mathbf{v}_{t_{n+1}}) + h\mathbf{g}_{t_{n+1}} = h\mathbf{H}^T \mathbf{\lambda}, \quad (2)$$

$$\mathbf{q}_{t_{n+1}} = \mathbf{q}_{t_n} + h\mathbf{v}_{t_{n+1}} \quad (3)$$

where $\mathbf{g}_{t_{n+1}}$ is the body force at time $t_{n+1}$ known a priori. The vector of internal forces $\mathbf{F}(\mathbf{q}_{t_{n+1}}, \mathbf{v}_{t_{n+1}})$ is a non-linear function of $\mathbf{q}$ and $\mathbf{v}$. Hence, we perform the following linearisation (which is a decent approximation over a small time interval, involving a small variation of the internal forces):

$$\mathbf{F}(\mathbf{q}_{t_{n+1}}, \mathbf{v}_{t_{n+1}})) = \mathbf{F}(\mathbf{q}_{t_n} + d\mathbf{q}, \mathbf{v}_{t_n} + d\mathbf{v}) \quad (4)$$

$$= \mathbf{f}_{t_n} + \frac{\delta \mathbf{F}}{\delta \mathbf{q}} d\mathbf{q} + \frac{\delta \mathbf{F}}{\delta \mathbf{v}} d\mathbf{v} \quad (5)$$

Using $d\mathbf{q} = \mathbf{q}_{t_{n+1}} - \mathbf{q}_{t_n} = h\mathbf{v}_{t_{n+1}}$ and $d\mathbf{v} = \mathbf{v}_{t_{n+1}} - \mathbf{v}_{t_n}$, we obtain:

$$\underbrace{(\mathbf{M} + h \frac{\delta \mathbf{F}}{\delta \mathbf{v}} + h^2 \frac{\delta \mathbf{F}}{\delta \mathbf{q}})}_{\mathbf{A}(\mathbf{q}_{t_n}, \mathbf{v}_{t_n})} d\mathbf{v} \quad (6)$$

$$= -h^2 \frac{\delta \mathbf{F}}{\delta \mathbf{q}} \mathbf{v}_{t_n} - h(\mathbf{g}_{t_{n+1}} + \mathbf{f}_{t_n}) + h \mathbf{H}^T \mathbf{\lambda}. \quad (7)$$

Finally, at each time step $t_n$, the equation can be written:

$$\mathbf{A}_{t_n} d\mathbf{v}_{t_{n+1}} = \mathbf{b}_{t_n} + h \mathbf{H}^T \mathbf{\lambda}. \quad (8)$$

### B. Collision model and contact law

1) *Collision model discretization:* To model collision, a collision discretisation is defined on the surface of the robot and will be used to assess they is no inter-penetration between the robot and its environment. This collision model is typically defined on a 2D triangular mesh covering the surface of the robot. This triangular mesh may differ from the 3D FE-mesh defined for computing the robot deformation but it is often matching the tetrahedral mesh, i.e. the vertices of the triangles match the vertices on the boundary of the 3D mesh. We assume this is always the case in the rest of the paper.<fnref n="1" /> A Lagrange multiplier $\mathbf{\lambda}$ is defined on each vertex, edge and triangle of the collision model. Depending on the contact case, only one or several of these may be active at the same time. Let us define the set:

$$\Lambda = \{\mathbf{\lambda}_p\}_{p \in \mathcal{C}}, \quad (9)$$

the dual space, containing all contact forces, with $\mathcal{C}$ the set of potential contact points (including vertices, edges and/or triangles). The cardinality of the set $\mathcal{C}$ defines the order of the dual space.

Now, we are introducing the Signorini law for contact. Consider two points, each on a surface of two distinct objects. The law states that there is a complementarity relation linking the gap $\delta$ between these two points and the contact force $\mathbf{\lambda}$ they apply on each other :

$$0 \le \delta_n \perp \lambda_n \ge 0, \quad (10)$$

where $\mathbf{\lambda}$ is decomposed into a normal scalar component $\lambda_n$ and a tangential vector component (in 3D) $\mathbf{\lambda}_t$ at contact.

<footnote n="1">¹Note that if that is not true, a mapping can be defined between collision model and deformation model (for example using barycentric coordinates).</footnote>
This law can be expressed in terms of velocity at the moment of contact as:

$$ \text{if } \delta_n = 0, 0 \le \dot{\delta}_n \perp \lambda \ge 0, \qquad (11) $$

which is handy since the motion in equation (8) is expressed in terms of velocity variation.

Friction may be modeled using Coulomb's law with stick and slip conditions: 
$$ \begin{cases} \text{if } \dot{\delta}_t = 0, \|\lambda_t\| < \mu|\lambda_n| \\ \text{if } \dot{\delta}_t \neq 0, \|\lambda_t\| = \mu|\lambda_n|, \exists a \ge 0 \text{ such that } \lambda_t = -a.\dot{\delta}_t, \end{cases} $$
with $\mu$ the friction coefficient.

### III. HYBRID HYPER-REDUCTION METHOD (HHR)

In this section, we present the Hybrid-Hyper-reduction method, first described in [11] and [12]. The method is hybrid in the sense that it reduces both primal unknowns (velocities) and dual unknowns (lagrange multipoles associated with contact constraints) in different ways. We describe the two stages of reduction in the following.

#### A. Primal reduction: Snapshot-POD with Hyper-reduction

a) **Projection onto the reduced space:** In the first part, a standard projection-based model order reduction is performed using proper orthogonal decomposition (POD). The robot displacement is assumed to be contained in a low-dimensional vector space, spanned by few basis vectors:

$$ \mathbf{q}(t) \approx \mathbf{q}(0) + \sum_{i=1}^{N} \phi_i \alpha_i(t) \qquad (12) $$

$$ = \mathbf{q}(0) + \mathbf{\Phi} \boldsymbol{\alpha}(t), \qquad (13) $$

The basis $\mathbf{\Phi} = \{\phi_0, \phi_1, ..., \phi_N\}$ is obtained from the application of a singular value decomposition (SVD) on a snapshot matrix $\mathbf{S}$ containing an exhaustive sample of the robots positions. This snapshot matrix is computed offline and that is a computationally intensive process. The SVD reads:

$$ \mathbf{S} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^T, $$

with $\mathbf{U}$ and $\mathbf{V}$ unitary matrices and $\mathbf{\Sigma}$ a rectangular diagonal matrix containing the singular values of $\mathbf{S}$. The optimal basis $\mathbf{\Phi}$ of order $p$ is then given by the first $p$ left singular vectors (in $\mathbf{U}$) associated to the $p$ largest singular values $\sigma_i$. The truncation error of a POD transform of order $p$ is given by:

$$ \nu^2 = \frac{\sum_{i=p+1}^{n_q} \sigma_i^2}{\sum_{i=1}^{n_q} \sigma_i^2} \qquad (14) $$

This formula provides a way to set the dimension of the reduced basis to reach a user-defined tolerance. Substituting $\mathbf{q}$ by its reduced expression (12) in the motion equation (8) at timestep $n$, and projecting onto the reduced space leads to the reduced motion equation:

$$ \mathbf{\Phi}^T \mathbf{A}_{t_n} \mathbf{\Phi} d\dot{\alpha}_{t_{n+1}} = \mathbf{\Phi}^T \mathbf{b}_{t_n} + h \mathbf{\Phi}^T \mathbf{H}^T \boldsymbol{\lambda}, \qquad (15) $$

with $d\alpha_{t_{n+1}} = \alpha_{t_{n+1}} - \alpha_{t_n}$ is the increment of the reduced displacement between timestep $n$ and $n + 1$, and $\dot{\alpha}$ its derivative with respect to time.

At this stage, it is important to note that the number of degrees of freedom of the robot have been greatly reduced, the dimension of equation (15) being much smaller than that of the original equation (8). However, the dimension of the contact space, materialised by the dual variable $\lambda$ remains identical.

b) **Hyper-reduction:** The reduction performed in the previous section is not enough to provide significant speedup because of the nonlinearities of the material. Indeed, the stiffness matrix is not constant and has to be recomputed at each timestep of the simulation, which is actually dominating the computational cost, since solving the reduced system (15) is fast. We apply a method of hyper-reduction [13], [14] called the energy-conserving weighting and sampling method (ECSW) [15], which allows for fast integration of a surrogate stiffness operator. The integration of the mechanical matrices will be only condensed on a small subset of the elements of the FE model called the sample mesh.

Defining $\mathbf{g} = \mathbf{A}\mathbf{\Phi}d\dot{\alpha} - \mathbf{\Phi}^T\mathbf{b}$ (dropping timestep dependence for clarity), $\mathcal{E}$ the set of all integration elements, $\tilde{\mathcal{E}} \subset \mathcal{E}$ the sample mesh, the ECSW method makes the approximation:

$$ \mathbf{\Phi}^T \mathbf{g} = \sum_{e \in \mathcal{E}} \mathbf{\Phi}_e^T \mathbf{g}_e \approx \sum_{e \in \tilde{\mathcal{E}}} \xi_e \mathbf{\Phi}_e^T \mathbf{g}_e, \qquad (16) $$

with $\mathbf{\Phi}_e$ and $\mathbf{g}_e$ the restrictions of $\mathbf{\Phi}$, respectively $\mathbf{g}$ to the element $e$, and the scalars $\{\xi_e\}_{e \in \tilde{\mathcal{E}}}$ positive weights. The ECSW method consider each term $\mathbf{\Phi}_e^T \mathbf{g}_e$ as the work done by the forces $\mathbf{g}_e$ along the displacement $\mathbf{\Phi}_e$. It hence guarantees an energy conservation by incorporating non-negative weights $\xi_e$ to compensate for the energy otherwise lost induced by the sparse integration on $\tilde{\mathcal{E}}$.

The selection of the sample mesh and the corresponding weights is done using a sub-optimal greedy algorithm based on minimising an error from a stored database generated offline just as it is done for the snapshot-POD procedure. See more details in [6], [15].

After applying hyper-reduction, equation (15) is replaced by a surrogate:

$$ \mathbf{\Phi}^T \tilde{\mathbf{A}}_{t_n} \mathbf{\Phi} d\dot{\alpha}_{t_{n+1}} = \mathbf{\Phi}^T \tilde{\mathbf{b}}_{t_n} + h \mathbf{\Phi}^T \mathbf{H}^T \boldsymbol{\lambda}, \qquad (17) $$

which typically preserves the accuracy of the reduced model. We will refer to this reduced model as the HR-ECSW model.

#### B. Dual reduction: HR-ECSW with contacts reduced integration

The reduction applied in the previous section only applies onto the primal variable, i.e. the displacement. The hybrid-hyper-reduction method proposes to extend this reduction to the dual space, i.e. contact forces $\Lambda$. Unlike the reduction on the primal variable, no reduced basis will be defined. Instead, the dual space $\Lambda$ will be reduced by taking its restriction to a sample mesh. Let us define $\Omega$ the original space of integration and $\tilde{\Omega}$ the subspace defined by the sampled mesh:

$$ \Omega = \bigcup_{e \in \mathcal{E}} e, \quad \tilde{\Omega} = \bigcup_{e \in \tilde{\mathcal{E}}} e. \qquad (18) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of primal and dual reduced space. The sampled mesh is displayed in blue as a subset $\tilde{\Omega}$ of the original full integration domain. The reduced collision space $\tilde{C}$ is defined as the intersection of the sampled mesh boundary and the full collision space.</figcaption></figure>

Let us now define the set $\tilde{C}$ as the intersection of the space defined by the sampled mesh and the collision model:

$$ \tilde{C} = \tilde{\Omega} \cap C = \partial\tilde{\Omega} \cap C, \qquad (19) $$

where $\partial\tilde{\Omega}$ is the boundary of the set $\tilde{\Omega}$. The hybrid hyper-reduction method proposes to define the reduced dual space as the restriction of the original dual space $C$ onto $\tilde{C}$:

$$ \tilde{\Lambda} = \Lambda|_{\tilde{C}}. \qquad (20) $$

The reduced dual space is a reduced integration domain (RID) for the contacts.

This definition implies that the reduced order of the space of contact forces is directly dependent on the sampled mesh of the HR-ECSW model. The space is illustrated in Figure 1. In practice, only triangles in full intersection with a tetrahedron from the sampled mesh will be included in the reduced collision space. We neglect the case where only one vertex or one edge is in intersection. In the following, we refer to this fully reduced model (velocities and contacts) as the HR-ECSW-RID model.

# IV. EXAMPLES

## A. Soft Bouncing Ball

As an academic example to illustrate the method, we apply the reduction on a soft bouncing ball onto a plane. A ball is discretised onto a tetrahedral mesh with 716 tetrahedra.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Different models of the ball. Top left and right: full order model and sampled model. Bottom left and right: full collision model and reduced collision model, the ball is turned for a better view.</figcaption></figure>

A triangular collision model matching the surface of the tetrahedral mesh is defined on the surface of this ball. This triangular model contains 364 triangles. The ball is supposed to deform following an elastic constitutive law defined by a Young Modulus of 1000 kPa and a Poisson Ratio of 0.45 (nearly incompressible). To generate the snapshot space, the ball is simply dropped onto a plane and positions are recorded every 5 timesteps. The storage stops once the ball has bounced. Friction is included with a value $\mu = 0.5$. Applying the reduction method explained in III, we obtain 12 modes with an error tolerance of $1.0e^{-3}$ including 3 translation modes to account for displacement without deformation. The sampled mesh is computed and contains 17 elements. The reduced collision space is then deduced as the intersection of this sampled mesh with the boundary which defines only 6 triangles on the surface of the ball. This is displayed in Figure 2. We can see in Figure 3 that the deformation of the ball onto the plane is preserved when using the reduced collision model. Indeed, even though the contact constraints are only partially applied, the deformation of the ball is constrained by the modes $\Phi$. Concerning computational time, the original ball model without any reduction runs at 25 frames per seconds (fps) on a standard laptop. It moves up to 105 fps applying the primal reduction, and 130 fps applying the full HR-ECSW-RID reduction. Note that this example is of relatively small dimension and greater speedups are expected with more complex examples.

## B. Soft Gripper

In this example, we define a simple gripper with 2 finger actuated by cable. We first apply the reduction onto a single finger and then duplicate it to create the gripper. The finger volumetric mesh contains 3302 vertices and 12398 tetrahedra. The collision mesh is made of 2450 vertices and
<figure><img src="image_4.png" /><figcaption>Fig. 3. Collision of the ball onto the plane using the HHR model: despite the "holes" in the collision model, the ball does not penetrate the plane and deforms correctly.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Multi-material decomposition of the finger. The top part (in blue) is much softer than the bottom part (in green)</figcaption></figure>

4896 triangles. For the constitutive law, we use a bi-material elastic model. The bottom part has a Young Modulus of 600 kPa and the top part is much softer with a Young Modulus of 20 kPa as pictured in Figure 4. The idea is to let the finger surface deform well to conform to a shape to be grabbed thus increasing contact surface and friction force.

After application of the SVD, we obtain 12 modes with a tolerance of $1.0e^{-3}$. The sampled elements are 46 for the integration domain. 17 triangles are saved in the collision model. The reduced model is displayed in Figure 6. We check the accuracy of the reduced model by measuring the position of the tip of the finger while the cube sweeps along the last phalanx. We compare both reduced models: primal reduction with full collision model and complete reduction of both deformation and collision. Results are displayed in Figure 7. Concerning computational time at the time of contact on a standard laptop computer, the full order model runs at 1.8 fps, the HR-ECSW model runs at 16 fps and the HR-ECSW-RID model runs at 40 fps. We can really see the benefit of reducing contact computation in this example since the model runs more than twice as fast with a small loss of accuracy.

We proceed to create a simple gripper robot by adding two instances of the reduced finger within a single simulation scene. The second finger is flipped to face the first finger. Its reduced basis (modes) are rotated to match the new position of the finger. We test the gripper with the task of grasping a rigid cube and a ball, see Figure 8. For the cube case, we compare the timings of the 3 models all along the sequence of grabbing. This is displayed in Figure 9.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Snapshot generation for a single finger with collision. The Finger is actuated with a cable and is in contact with a cube that is moved along the finger last phalanx.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Reduced model of the finger. We can see that most elements of the reduced integration domain are located at the tip. This was expected since the snapshot was generated with the tip in contact with the cube.</figcaption></figure>

## V. CONCLUSIONS

In this paper we have presented a method of Hybrid Hyper-Reduction and have applied it to an academic example of a bouncing ball and gripper. We showed the method can cut off computational time drastically while keeping a good accuracy. This opens up the path for quasi real-time simulation of soft robots with contacts that could be used for design or control tasks. More importantly this could be very beneficial for feeding reinforcement learning (RL) methods which has shown its potential for rigid robotics applications but is harder to apply for soft robots. This could be used in the context of grasping or locomotion. However, some improvements still have to be made. First of all, the method should be validated on more complex scenarios. The HHR
<figure><img src="image_8.png" /><figcaption>Fig. 7. Error of the reduced models along the simulation. We can see that this error fluctuates but remains below a threshold of half a centimeter. The error of the full HR-ECSW-RID model has higher peaks, certainly due to jumps between different parts of the reduced collision model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Soft gripper in action grabbing a cube and a ball. The gripper model is reduced following the HR-ECSW-RID method.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Evolution of the computational time for each timestep along the gripper simulation(log scale). Contacts appear after 10 timesteps and we can see a dramatic increase in computational time, especially for the full order model. The HR-ECSW model mitigates that increase. The HR-ECSW-RID model roughly cuts off computational time by 100 in that example.</figcaption></figure>

method is relying on an exhaustive sampling phase in the offline stage. This sampling can typically be intractable for accounting all contact cases in complex scenarios. Another issue is the treatment of self-collisions for which the method needs to be tuned. Indeed, one needs to make sure the reduced collision space is consistent when facing itself: triangles in the collision space on each side need to face each other. Nevertheless, the method provides a first automatic way of generating reduced models of soft robots accounting for contacts.

ACKNOWLEDGMENT

This work was funded by the French National Research Agency : project ANR-21-CE45-0031 SPECULAR.

REFERENCES

[1] J. Hiller and H. Lipson, "Dynamic Simulation of Soft Multimaterial 3D-Printed Objects," *Soft Robotics*, vol. 1, no. 1, pp. 88-101, 2014.

[2] Y. Hu, J. Liu, A. Spielberg, J. B. Tenenbaum, W. T. Freeman, J. Wu, D. Rus, and W. Matusik, "Chainqueen: A real-time differentiable physical simulator for soft robotics," in *2019 International conference on robotics and automation (ICRA)*. IEEE, 2019, pp. 6265-6271.

[3] A. Spielberg, T. Du, Y. Hu, D. Rus, and W. Matusik, "Advanced soft robot modeling in chainqueen," *Robotica*, vol. 41, no. 1, pp. 74-104, 2023.

[4] E. Coevoet, T. Morales-Bieze, F. Largilliere, Z. Zhang, M. Thieffry, M. Sanz-Lopez, B. Carrez, D. Marchal, O. Goury, J. Dequidt, et al., "Software toolkit for modeling, simulation, and control of soft robots," *Advanced Robotics*, vol. 31, no. 22, pp. 1208-1224, 2017.

[5] H. Courtecuisse, H. Jung, J. Allard, C. Duriez, D. Y. Lee, and S. Cotin, "Gpu-based real-time soft tissue deformation with cutting and haptic feedback," *Progress in biophysics and molecular biology*, vol. 103, no. 2-3, pp. 159-168, 2010.

[6] O. Goury and C. Duriez, "Fast, generic, and reliable control and simulation of soft robots using model order reduction," *IEEE Transactions on Robotics*, vol. 34, no. 6, pp. 1565-1576, 2018.

[7] G. Fang, Y. Tian, A. Weightman, and C. C. Wang, "Collision-aware fast simulation for soft robots by optimization-based geometric computing," in *2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*. IEEE, 2022, pp. 12614-12621.

[8] H. Dang, J. Weisz, and P. K. Allen, "Blind grasping: Stable robotic grasping using tactile feedback and hand kinematics," in *2011 IEEE international conference on robotics and automation*. IEEE, 2011, pp. 5917-5922.

[9] E. Coevoet, A. Escande, and C. Duriez, "Optimization-based inverse model of soft robots with contact handling," *IEEE Robotics and Automation Letters*, vol. 2, no. 3, pp. 1413-1419, 2017.

[10] O. Goury, B. Carrez, and C. Duriez, "Real-time simulation for control of soft robots with self-collisions using model order reduction for contact forces," *IEEE Robotics and Automation Letters*, vol. 6, no. 2, pp. 3752-3759, 2021.

[11] J. Fauque, I. Ramière, and D. Ryckelynck, "Hybrid hyper-reduced modeling for contact mechanics problems," *International Journal for Numerical Methods in Engineering*, vol. 115, no. 1, pp. 117-139, 2018.

[12] S. Le Berre, "Hyper-réduction d'ordre de modèle en mécanique du contact sous formulation mixte potentiellement non symétrique: analyse théorique et numérique," Ph.D. dissertation, Université Paris sciences et lettres, 2023.

[13] D. Ryckelynck, "A priori hyperreduction method: an adaptive approach," *Journal of Computational Physics*, vol. 202(1), pp. 346-366, 2005.

[14] J. A. Hernandez, M. A. Caicedo, and A. Ferrer, "Dimensional hyper-reduction of nonlinear finite element models via empirical cubature," *Computer methods in applied mechanics and engineering*, vol. 313, pp. 687-722, 2017.

[15] C. Farhat, P. Avery, T. Chapman, and J. Cortial, "Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency," *International Journal for Numerical Methods in Engineering*, vol. 98, no. 9, pp. 625-662, 2014.