A Validated Physical Model For Real-Time Simulation of Soft Robotic Snakes

TL;DR

This paper introduces a validated physics-based model for real-time simulation of soft robotic snakes, enabling accurate and interactive control of complex soft actuators using GPU acceleration.

Contribution

A novel constraint-based dynamics model for pneumatic soft actuators that accounts for internal pressure, latency, and damping, validated on a full soft robotic snake.

Findings

Model accurately predicts static deformation.

Model effectively simulates dynamic locomotion.

GPU acceleration enables real-time interactive control.

Abstract

In this work we present a framework that is capable of accurately representing soft robotic actuators in a multiphysics environment in real-time. We propose a constraint-based dynamics model of a 1-dimensional pneumatic soft actuator that accounts for internal pressure forces, as well as the effect of actuator latency and damping under inflation and deflation and demonstrate its accuracy a full soft robotic snake with the composition of multiple 1D actuators. We verify our model's accuracy in static deformation and dynamic locomotion open-loop control experiments. To achieve real-time performance we leverage the parallel computation power of GPUs to allow interactive control and feedback.

Tables

Table 1. TABLE I: size of the structure for one simulated snake

Type	Quantity
Rigid Bodies	15
Particles	1504
Distance Constraints	1460
Tetrahedral Finite Elements	4536
Rigid Joints	10
Particle Attachments	217

Table 2. TABLE II: Benchmark results for our simulator

Snakes (#)	$1/4$	1	2	3	4	5	6	7	8	9	10
Assembly (ms)	0.74	1.11	1.41	2.33	3.41	5.59	6.19	6.45	7.37	8.48	10.31
Solve (ms)	4.96	10.52	25.06	31.90	43.12	50.41	58.63	66.95	79.73	86.89	95.57
Total (ms)	5.70	11.63	26.47	34.23	46.53	56.00	64.82	73.4	87.10	95.37	105.88
Total/Snake (ms)	-	11.63	13.23	11.41	11.63	11.2	10.80	10.48	10.88	10.59	10.58

Equations

$$\displaystyle\mathbf{M}\ddot{\mathbf{q}}-\mathbf{f}(\mathbf{q},\dot{\mathbf{q}})-\mathbf{J}_{b}^{T}\mbox{\boldmath$\lambda$}_{b}-\mathbf{J}_{n}^{T}\mbox{\boldmath$\lambda$}_{n}-\mathbf{J}^{T}_{f}\mbox{\boldmath$\lambda$}_{f}=\mathbf{0}$$

$$\displaystyle\mathbf{c}_{b}(\mathbf{q},\mathbf{p})+\mathbf{E}\mbox{\boldmath$\lambda$}_{b}=\mathbf{0}$$

$$\displaystyle\mathbf{0}\leq\mathbf{c}_{n}(\mathbf{q})\perp\mbox{\boldmath$\lambda$}_{n}\geq\mathbf{0}$$

$$\displaystyle\forall i\in\mathcal{A},\quad\mathbf{D}_{i}^{T}\dot{\mathbf{q}}+\frac{|\mathbf{D}_{i}^{T}\dot{\mathbf{q}}|}{|\mbox{\boldmath$\lambda$}_{f,i}|}\mbox{\boldmath$\lambda$}_{f,i}=\mathbf{0}$$

$$\displaystyle\forall i\in\mathcal{A},\quad 0\leq|\mathbf{D}_{i}^{T}\dot{\mathbf{q}}|~{}~{}\perp~{}~{}\mu_{i}\lambda_{n,i}-|\mbox{\boldmath$\lambda$}_{f,i}|\geq 0$$

$$\displaystyle\forall i\in\mathcal{I},\quad\mbox{\boldmath$\lambda$}_{f,i}=\mathbf{0}.$$

$${\mathbf{q}}_i={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T.$$

$$\displaystyle\mathbf{q}_{i}=\begin{bmatrix}\mathbf{x}_{i}^{T}&\mbox{\boldmath$\theta$}_{i}^{T}\end{bmatrix}^{T}.$$

$$c_{dist}(\mathbf{q},p)=|{\mathbf{q}}_i-{\mathbf{q}}_j|-rϵ(p)=0,$$

$$ϵ(p)=1+p/Y.$$

$$\mathbf{c}_{tetra}(\mathbf{q})+\mathbf{E}_{tetra}\mbox{\boldmath$\lambda$}=\mathbf{0},$$

$${\mathbf{E}}_{tetra}=\frac{1}{V_{e}Y}\left[\begin{array}{cccccc}1& -\nu & -\nu & 0& 0& 0\\ -\nu & 1& -\nu & 0& 0& 0\\ -\nu & -\nu & 1& 0& 0& 0\\ 0& 0& 0& 1+\nu & 0& 0\\ 0& 0& 0& 0& 1+\nu & 0\\ 0& 0& 0& 0& 0& 1+\nu \end{array}\right].$$

$${\mathbf{c}}_{attach}(\mathbf{q})={\mathbf{q}}_x+\mathbf{R}({\mathbf{q}}_{\theta })\mathbf{r}-{\mathbf{q}}_p=0.$$

$$c_n(\mathbf{q})={\mathbf{n}}^T\left[\mathbf{a}(\mathbf{q})-\mathbf{b}(\mathbf{q})\right]\ge 0,$$

$$\displaystyle\left[\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^{T}+\mathbf{E}\right]\Delta\mbox{\boldmath$\lambda$}=\mathbf{b}.$$

$$a_i\equiv \min (1,\max (-1,(\sin (\omega t+\alpha i)+\varphi ))A.$$

$$v^2=\frac{2}{\rho }(p_t-p_s),$$

$$\Delta p_i=\frac{a_i(t+h)-p_i(t)}{p_s},$$

$$p_i(t+h)=\{\begin{array}{ll}{p}_i(t)+p_s{\Delta p_i}^2{k}_i& \text{is inflating}\\ p_i(t)-\min (p_i(t)k_d,T_p)& \text{is deflating}\end{array}$$

Full text

I Introduction

Soft robotic systems can provide many advantages over rigid robots. For example, the compliance of a soft manipulator allows safe interaction with fragile objects, such as uncooked eggs [1, 2]. A soft robotic snake can navigate through narrow passages to perform tasks in regions unreachable with rigid robots [3]. On the other hand, soft robots are harder to model and simulate than rigid robots, due to their infinite degrees of freedom. Thus, one often has to compromise between simulation accuracy and simulation time. In recent work [4, 5], the authors developed a framework of deformable online simulation that can achieve real-time simulation using relatively coarse meshes. Further, the authors proposed a model order reduction method by running expensive offline simulations and applying machine-learning techniques on the generated dataset [6].

However, the compromise between accuracy and simulation time is often hard to make. An accurate model is important to bridge the gap between simulation and reality for control design. On the other hand, online trajectory planning and model predictive control require the prediction of future states with low-latency. In addition, even learning-based controllers, such as those trained with Reinforcement Learning (RL) [7, 8, 9], typically require thousands of samples to converge, making efficient simulation crucial. Approximate models can enable faster simulations; Furthermore, even in rigid-body simulators, the modeling errors, inaccuracies, and lack of actuator dynamics on the simulator create a gap to the reality that often needs to be bridged before applying controllers learned in simulation to the real robot [10].

In this work, we aim to develop a real-time, high-fidelity simulation for soft robotic systems. We present a compliant constraint dynamics model for a soft snake robot [11, 12, 13]. Our model is able to accurately represent the deformations of the snake links while being efficient enough for real-time simulations. To achieve low-latency simulation, we use a GPU-based physics simulator that leverages large-scale parallel iterative solvers to efficiently solve large systems [14]. The resulting simulation is validated against a real robotic snake to verify that the deformation model is accurate and that the dynamics of the simulation match that of the real system. In addition, we measure and validate the actuator and inflation latency to ensure accurate control response. In summary, our main contributions are:

• •

A dynamics and actuation framework for 1D pneumatic soft actuators that accurately represents a large range of deformations;

• •

A model for a modular soft robotic snake that accurately represents its dynamics;

• •

Model validation in static deformation and dynamic locomotion tests;

• •

A simulator framework suitable for performing real-time control of soft robots.

II Related Work

The field of soft robots includes a large and varied range of designs that incorporate compliant materials and actuators in various forms. Soft robots may have a few flexible regions in them and may be driven by tendons [15], while completely soft pneumatic actuators can be driven by exerting a range of pressures within deformable bodies. The use of pressure leads to many intricate designs that exploit the material geometry to achieve the desired actuation [1, 16]. To allow a large range of pressures and material deformations, hybrid materials are often used on the pressure chambers, such as inextensible layers and fiber reinforcements [17]. Our modular soft snake robots use hybrid materials [12, 13] to achieve highly efficient 2D terrain locomotion.

The diversity of soft materials and difficulty in accurate modeling of soft robots lead to an increasing need for building a realistic simulation for deformable robots. Duriez et al. [18] presented a framework for simulating soft robots using a quasi-static approach based on Finite Elements (FEM). In this work, we also use FEM with tetrahedral elements as a fundamental building-block. In addition, we combine these elements in a multiphysics system with spring networks, frictional contact, attachment constraints for soft/hybrid materials and articulated rigid bodies. We perform implicit time-integration to simulate dynamic trajectories. Recent work by Pozzi et al. [19] used a rigid-body model, fitted to an offline FEM simulation augmented with stiff springs to achieve real-time updates. In this work, we simulate the finite-element models and spring networks directly, using the large-scale parallelism of graphics processing units (GPUs) to achieve online update rates, therefore not requiring an expensive offline simulation and data fitting stage.

Tan et al. [20] propose several techniques for bridging the reality gap in simulations. Among them there is a latency introduced for the delay between when a command is sent and when it’s executed. Since soft actuators have an over-damped response, our work introduces system dynamic model response to the command on top of the latency.

III Soft Robotic Snake

The snake is made of soft bending actuation modules with integrated curvature sensing [11], as shown in Fig. 1c. Each soft bending actuator segment is comprised of two soft linear actuators and an inextensible constraint layer between them. The soft links are made of silicone rubber with cavities wrapped in a fiber reinforcement, which limits the radial deformation when pressurized. In the center of the link, between the two chambers, there is a custom integrated curvature sensor, and a plastic film that inhibits linear extension. This set of constraints results in bending the entire soft module when one chamber is pressurized, as seen in Fig. 1b. The curvature sensor utilizes a magnet and a Hall effect sensor, mounted on a flexible circuit board [11]. Caps are attached to both ends of the actuator to seal the chambers and allow modular connections with other segments. The caps are made of two ABS plates sandwiching the rim of the silicone. The actuator is driven by two 3-2 (3-port, 2-state) binary solenoid valves, each connecting one pressure chamber to a common pressure source. The valves can either inflate or deflate a given actuator chamber. The pressure in each chamber is controlled using a 60 Hz Pulse-Width Modulation (PWM) on the valves, which are set to operate in complete antagonism, so when one chamber is inflating the other is always deflating. This reduces the number of required inputs per chamber to one, corresponding to the single (active) DOF of the bending actuator.

IV Simulator

The continuous equations of motion for the multiphysics simulator are derived from Lagrangian mechanics, and are given in general form by the following,

[TABLE]

These equations describe the motion of a generic dynamics system with frictional contact forces. The state of the system is described by a vector of generalized coordinates $\mathbf{q}\in\mathbb{R}^{n_{d}}$ with $n_{d}$ DOFs, determined by the number of particles and rigid bodies on the system. The inertial properties of the system are represented by the mass-matrix $\mathbf{M}\in\mathbb{R}^{n_{d}\times n_{d}}$, with $\mathbf{f}(\mathbf{q},\dot{\mathbf{q}})$ a generalized force function that includes external and gyroscopic forces. The vector $\mathbf{c}_{b}(\mathbf{q})$ is a set of bilateral constraints of length $n_{b}$, with \mbox{\boldmath\lambda}_{b} the associated Lagrange multipliers. Elastic energy potentials are defined in terms of compliant constraints, here $\mathbf{E}\in\mathbf{R}^{n_{b}\times n_{b}}$ is a block-diagonal compliance, or inverse stiffness matrix as described by Servin et al. [21]. The target pressures are grouped into the vector $\mathbf{p}$, which are parameters to the actuation constraints described in section IV-C. The contact and frictional forces are based on Coulomb’s model, which defines an admissible cone of contact forces [22]. Here $\mathbf{c}_{n}(\mathbf{q})$ are unilateral contact constraints, with $n_{c}$ the number of contacts in the system, and $\lambda_{n,i}$ and $\mu_{i}$ the normal force Lagrange multiplier and friction coefficient for the $i$th contact respectively. The frictional forces for a contact are parameterized by \mbox{\boldmath\lambda}_{f,i}, with a corresponding basis $\mathbf{D}_{i}$ that defines the surface tangent plane at the contact point. The active contact set is defined as $\mathcal{A}=\{i\in(1,\cdots,n_{c})\mid\mu_{i}\lambda_{n,i}>0\}$, with inactive contacts $\mathcal{I}$ being its complement. The constraint Jacobians $\mathbf{J}_{b},\mathbf{J}_{n}$ contain the gradient of bilateral and normal constraint functions with respect to $\mathbf{q}$, and we define the set of frictional basis vectors as the matrix $\mathbf{J}_{f}=[\mathbf{D}_{1},\cdots,\mathbf{D}_{n_{c}}]^{T}$.

IV-A Particles

Each deformable link is modeled as a collection of particles connected by constraints. This is a flexible representation that allows fine-grained control over different sections of the soft body, while being efficient enough for real-time simulation. A particle with index $i$ adds three additional DOFs to the system,

[TABLE]

Assuming a lumped mass model, each particle is assigned a fraction of the connected tetrahedral elements mass (Section IV-C2). The mass-block for the particle is then given by $\mathbf{M}_{i}=m\mathbf{1}_{3}$, where $m$ is the particle mass, and $\mathbf{1}_{3}$ is the $3$-dimensional identity matrix.

IV-B Rigid Bodies

We describe the state of a rigid body with index $i$ using a maximal coordinate representation consisting of the position of the body’s center of mass, $\mathbf{x}_{i}\in\mathbb{R}^{3}$, and its orientation expressed as a quaternion \mbox{\boldmath\theta}_{i}=[\theta_{1},\theta_{2},\theta_{3},\theta_{4}]^{T}\in\mathbb{R}^{4}. We group these components together so the state sub-block for a single rigid body is

[TABLE]

IV-C Constraints

IV-C1 Actuation Constraints

To perform actuation we constrain particles together through equations of the form,

[TABLE]

where $\mathbf{q}_{i}$ and $\mathbf{q}_{j}$ are particle positions, and $r$ is a rest length to maintain between them. The target pressure $p$ induces a strain $\epsilon(p)$ that adjusts the rest length and causes contraction or expansion. Assuming that deformation is linear with stress, and that it occurs primarily along the chamber’s main axis, the amount of expansion/contraction is given by the following relation between material stiffness determined by the Young’s Modulus ($Y$) and pressure $p$,

[TABLE]

Furthermore, we use distance constraints with constant rest length to model the structural stiffness in the deformable chamber, as described in Sec. V.

IV-C2 Tetrahedral Finite-Elements

In addition to distance constraints, tetrahedral finite-elements are used to model the solid chamber material. Assuming a constant strain element and a linear isotropic constitutive model, each tetrahedron defines a 6-dimensional constraint vector,

[TABLE]

where $\mathbf{c}_{tetra}(\mathbf{q})=[\epsilon_{xx}~{}\epsilon_{yy}~{}\epsilon_{zz}~{}\epsilon_{yz}~{}\epsilon_{xz}~{}\epsilon_{xy}]^{T}$ is the vector of corotational strains in Voigt notation, and $\mathbf{E}_{tetra}$ is the constant element compliance matrix, given by

[TABLE]

where $V_{e}$ is the element volume, $Y$ and $\nu$ are the material Young’s modulus and Poisson’s ratio, respectively.

IV-C3 Rigid Body to Particle Attachment

In order to connect soft links to rigid bodies, an attachment constraint between a particle and a point on a rigid body is defined as follows,

[TABLE]

This is a vector-valued function that adds three separate constraints, one for each $x,y,z$ axis respectively. Here $\mathbf{q}_{x},\mathbf{q}_{\theta}$ are the rigid body position and orientation respectively, and $\mathbf{q}_{p}$ is the particle position. $\mathbf{R}(\mathbf{q}_{\theta})$ is the rotation matrix obtained from the body’s orientation, and $\mathbf{r}$ is the attachment point expressed in the body’s local frame.

IV-C4 Rigid Body Joints and Contact

Along with the deformable sections, we model the articulated carriage as rigid bodies, with wheels connected to the main frame using hinge joints as described by [23]. Contact between the wheels and the ground is modeled by non-interpenetration constraints of the form:

[TABLE]

where $\mathbf{n}\in\mathbf{R}^{3}$ is the contact plane normal, $\mathbf{a}$ and $\mathbf{b}$ are points on a rigid or deformable body. Frictional forces are included using a Coulomb model derived from a principle of maximum dissipation that limits the contact forces to a cone. We refer the readers to the survey paper by Stewart [24] for more detail.

IV-D Time-Stepping

The simulation is advanced in time with a first-order implicit time-discretization of the equations of motion similar to the method in [25]. An implicit discretization is chosen as it allows taking large time-steps and avoids constraint drift. At each time-step, the nonlinear system of equations resulting from the implicit discretization is solved using Newton’s method. To solve the complementarity conditions associated with contact we use a non-smooth reformulation based on the Fischer-Burmeister function as described in [26]. Each Newton iteration requires the solution of a sparse-matrix equation of the form

[TABLE]

Where $\mathbf{J}=[\mathbf{J}_{b}^{T}~{}\mathbf{J}_{n}^{T}~{}\mathbf{J}_{f}^{T}]^{T}$ is the matrix of constraint Jacobians, $\mathbf{E}$ is a block-diagonal compliance matrix that includes the tetrahedral compliance matrices, and $\mathbf{b}$ includes the constraint function residuals evaluated at the current Newton iterate. This is a positive semi-definite system that we solve using the Preconditioned Conjugate Residual method (PCR) [27]. This is an iterative Krylov method similar to Conjugate Gradient (CG) but with smooth error reduction, making it better suited for real-time applications with a fixed computational budget. Like CG, the primary computation cost of PCR is the performing sparse matrix-vector multiplications. However, these multiplications are highly parallelizable, and can be done efficiently by assembling $\mathbf{J},\mathbf{M},\mathbf{E},\mathbf{b}$ on the GPU in compressed row-storage (CSR) format, and performing the multiplication with optimized kernels [14]. In our simulator we use a simple diagonal Jacobi preconditioner since it is trivial to parallelize.

V Soft robotic snake model

The soft links of the snake robot are made of Ecoflex™ 00-30 silicone rubber which has material parameters $Y=66.243\text{KPa}$, and $\nu=0.4999$ [28]. We construct a triangular mesh of the surface and tetrahedralize it using TetGen [29]. The link mesh was created with evenly distributed particles, we do not explicitly represent the cavity with tetrahedral elements. The mesh was carefully constructed to provide a radially symmetrical tetrahedron structure, as seen in Fig. 3.

Since the inextensible layer in the center of the link has a deformation threshold that is beyond the range of forces to be applied on the soft links, it is acceptable to model it as a non-deformable constraint between particles along the center plane. Similarly, the radial constraint on the chambers are defined as a set of inter-particle constraints over coplanar particles along the link. Although it would be possible to drive each link’s expansion using surface pressure forces directly, the other constraints in the link allow us to simply control the chamber volume using constraints between particles along the primary axis of expansion.

Only one chamber on the link is active (i.e. pressurized) at a time, so this set of actuation constraints only applies to the expanding chamber. Fig. 3 displays the constraints overlaid on the link. The link mesh was subdivided in 13 cross-sections along its length, in order to allow real-time computation, while maintaining good accuracy on the material deformation.

The links are then connected to each other through the rigid bodies that contain the electronics necessary to control the snake robot. In addition, the rigid bodies are attached to the wheels via. hinge joints. The wheels provide contact with the floor and model the anisotropic friction that a real snake has from its scales.

The type and number of all DOFs, and constraints in the simulated snake are displayed in Table I.

V-A Open-Loop Control

The snake assembly consists of four links attached together, as seen in Fig. 1. An undulating motion that propels the snake is given by the control equation 9, which outputs the pressure for the link $i$.

[TABLE]

If $a_{i}$ is positive, the link will inflate one chamber of the link, while if $a_{i}$ is negative, it will inflate the other. The parameters $\omega$, $\alpha$, $\phi$, and $A$ are the base oscillation frequency, measured in $Hz$, the phase shift between links ($radians$), the offset value for the motion ($[-1,1]$, scalar), and the oscillation amplitude ($[0,p_{max}]$, Pa), respectively. The oscillation frequency dictates how fast the actuators will switch direction, and the phase delay between links is what generates the wave pattern that propels the snake forward. These parameters set the base undulating motion. By changing the offset $\phi$, applied to all links, the controller will inflate one side for longer than the other. This results in the snake propelling itself more to the opposite direction than the chambers that are more inflated, making the trajectory of the center of mass moves with a curvature radius determined by this offset and the friction with the ground. Finally, the amplitude $A$ limits the maximum pressure during the oscillation, thereby controlling the snake’s speed. The parameters that make the snake move forward depend on the physical properties of the snake, such as weight, length, friction coefficient with the floor, and were determined experimentally in [11]. This controller is an open-loop method, which generates the forward motion and allows to make turns, but no feedback is given if the trajectory is deviating from the desired trajectory.

The pressure delivery is not instantaneous but limited by the maximum airflow allowed by the valves. Assuming the pressure source can reliably maintain a constant output pressure $p_{s}$, the air flow $v$ to the chamber is given by [30],

[TABLE]

where $\rho$ is the air density and $p_{t}$ is the pressure in the chamber. This means the pressure update is proportional to the square of the difference between the current pressure and the desired pressure. The pressure update for inflation in each step is then defined based on the difference $\Delta p_{i}$, in Eq. 11.

[TABLE]

The deflation releases pressure in the atmosphere while keeping its own pressure relatively constant due to the change of volume. Therefore, the deflation ratio should be close to linear up to a threshold $T_{p}$ when it is proportional to the over-pressure with a damping,

[TABLE]

where $k_{i},k_{d}\in(0,1]$ are the inflation damping parameters, and are tuned according to the experimental data.

VI Experimental Verification

For all tracking experiments, the poses of the rigid links were captured using the motion capture system (MOCAP) by placing four markers on each rigid extremity so that it can collect their full poses. The MOCAP system contains 11 cameras surrounding the observable space of $4\times 3~{}m^{2}$. This redundancy allows the information of links to be collected with high precision and minimal loss of tracking during the experiments. In order to eliminate remaining outliers, every experiment is repeated for 10 times unless otherwise mentioned. To ensure data sanity, every sample that deviates more than an expected maximum displacement from the collected values is pruned and not used in the analysis, as this type of samples results in a data collection with negligible standard deviation. Fig. 1 shows the markers on the top corners of the rigid plates.

A friction coefficient of $\mu_{i}=1$ has been used for all experiments. We implement our simulator in CUDA and run it on a computer with an Intel Core-i7 8520K, 32GB of RAM, and one NVIDIA GTX 1080 Ti GPU. We use a fixed time-step of $h=0.0083s$, each time-step performs 4 Newton iterations, with each linear system solved approximately using 20 PCR iterations to ensure a fixed computational cost.

VI-A Quasi-static Verification

The first experiment on the simulation is to verify whether the pressure actuator follows the same geometrical behavior as the real link. For this experiment, the curvatures of the real links were obtained by subtracting the yaw of the rigid connectors attached to each soft link, for the varying over-pressures (the pressure that exceeds the resting atmospheric pressure) from [math] to $10$ psi, moving up with steps of $1$ psi for both directions on the link, and averaged over $300$ samples. Negative values show the inflation of left chamber and positive values are for the right chamber.

From Fig. 5, it can be seen that the curvature increases linearly with the pressure within a range. Particularly, the spring model is able to closely match the real curvature, and accurately follows the linear model up to 8psi. Ecoflex 00-30 Young modulus is 66 kPa, at 9 psi (62 kPa) it nears 2x expansion, which is the limit at which the material is linear. This behavior is also clearly observed by the relative error plot. When the pressure exceeds 8 psi, the real link starts to bulge over imperfections in the manufacturing, and on the opposite side, it folds in itself, resulting in a deviation from the linear model. Besides, since the extra pressure is forcing over the imperfections of the manufacturing process, there is a potential risk of damaging the links in the long term. For these reasons, it was deemed that the safe pressure threshold shall be 8 psi, and all the remaining tests were restricted to up to that range.

VI-B Dynamic Verification

The dynamic verification starts with the analysis to the step response on the actuators. A single link was used to capture the rate at which it inflates and deflates from rest to $60-100\%$ at $10\%$ steps. These trials were used to tune the gains $k_{i}=k_{d}=0.23$ and $T_{p}=0.68$. The results are seen in Fig. 6.

The simulator trajectory is then tested with an open-loop generator from Eq. 9 and compared with the real snake robot using the following parameters $\omega=2Hz$, $\alpha=\pi/2$, $\psi=0$ and $A=p_{max}=8psi$. Since it’s an open-loop control, it is expected that due to model inaccuracies and unmodeled dynamics the simulator will have error accumulated along the trajectory. As a result, the simulated trajectory may diverge from the actual trajectory over time. The comparison was made using the center of mass of the snake. The divergence can be addressed with the use of closed-loop feedback control for both simulated snake and the real snake robots. However, from Fig. 7, it can be seen that starting from the same initial condition, the simulated trajectory closely matches the real trajectory for the execution time. The inclusion of the latency model for the pressure update makes the trajectory of the center of mass go in the same overall direction and amplitude as the real snake.

VI-C Benchmarking

In order to test scalability of the system, we benchmark simulation times for a single soft link, a full snake with 4 links, and up to 10 snakes, each with 4 links. The linear system for a single Newton iteration corresponding to a single snake is a sparse matrix of roughly $30000\times 30000$. Each frame is sub-divided in two sub-steps for the constraint solver, and the optimization runs with 20 iterations each. The total per-frame simulation times are in Table II. From the results, it can be seen that simulation time increases linearly with the number of links and snakes. Linear scaling is expected, since a single link saturates the GPU. On average, it takes less than $12\text{ ms}$ to simulate each snake, with optimal performance obtained when simulating seven snakes together. One single snake takes $8.9s$ per frame running at 3.3GHz on an Intel Core i7 5820k.

VII Conclusion

We have presented a dynamical model for simulating 1D pneumatic actuators and a framework to simulate it in a multi-physics environment in real time. By comparing the simulation with a real soft robotic snake, it is demonstrated that the simulated snake produces real-time high-fidelity results even in complex scenarios involving a mixture of hybrid soft bodies, rigid bodies, and friction contacts. In the open loop control analysis, it becomes clear that the simulator doesn’t perfectly model the real dynamics; however, it remains in close performance to the real snake. Our next step is now apply a control trained and tested in simulation on the real snake. The transfer from simulation to reality, while still maintaining stability, can also be facilitated by domain randomization techniques [20]. We demonstrated the use of GPU in accelerating high-fidelity simulation for soft robots and present a framework that can generalize to other soft robotic systems. Our future research is to develop learning-based control to train fast and stable snake gaits in a range of terrains, as well as obstacle-aided navigation, and learning specific motion primitives from demonstration [31, 32]. Another step is also study the model different types of soft actuators such as PneuNet [2], and plan to publish data for our snake model to enable an open platform for experimentation and improvement.