Differentiable Simulation of Soft Multi-body Systems
Year: 2022

Abstract:
We present a method for differentiable simulation of soft articulated bodies. Our work enables the integration of differentiable physical dynamics into gradient-based pipelines. We develop a top-down matrix assembly algorithm within Projective Dynamics and derive a generalized dry friction model for soft continuum using a new matrix splitting strategy. We derive a differentiable control framework for soft articulated bodies driven by muscles, joint torques, or pneumatic tubes. The experiments demonstrate that our designs make soft body simulation more stable and realistic compared to other frameworks. Our method accelerates the solution of system identification problems by more than an order of magnitude, and enables efficient gradient-based learning of motion control with soft robots.

heading: 1 Introduction
text: Soft articulated bodies have been studied and utilized in a number of important applications, such as microsurgery [32], underwater robots [37], and adaptive soft grippers [20]. Since the compliance of deformable materials can enable robots to operate more robustly and adaptively, soft biomimetic robots are drawing a lot of attention and have made considerable progress. A snailfish robot dives at a depth of 10,900 meters in the Mariana Trench [43]. Drones equipped with soft manipulators grasp and transmit objects with a 91.7% success rate [17]. Soft hands with pneumatic actuators are able to grasp objects of different shapes, including water bottles, eyeglasses, and sheets of cloth [12]. To enable rapid prototyping of soft robots and efficient design of control algorithms through virtual experiments, we aim to create a realistic deformable multi-body dynamics framework, in which soft articulated robots can be simulated to learn powerful control policies.
Design and control of soft robots are challenging because of their nonlinear dynamics and many degrees of freedom. Differentiable physics has shown great promise to deal with such complex problems [4, 14, 31, 68]. One possibility is to treat soft bodies as volumes that are modeled as sets of particles or finite elements [29, 16]. These methods have made great progress, but the volumetric representations are difficult to scale to large multi-body systems and are poorly suited to modeling internal skeletons. Moreover, contact handling in recent differentiable physics frameworks [51, 30] often does not comply with Coulomb’s Law, which is central to plausible visual realism and correct physical behavior.
In this paper, we design a powerful and accurate differentiable simulator for soft multi-body dynamics. Since our entire framework is differentiable, our method can be embedded with gradient-based optimization and learning algorithms, supporting gradient-based system identification, motion planning, and motor control. Within the simulator, we first use tetrahedral meshes to enable adaptive resolution and more accurate modeling. Next, to couple soft materials with articulated skeletons, we design a
∗Equal contribution.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
top-down matrix assembly algorithm within the local steps of Projective Dynamics [5]. For accurate contact handling, we extend and generalize a dry friction model previously developed for cloth simulation [54] to soft solids and introduce a new matrix splitting strategy to stabilize the solver. In addition, our simulation framework incorporates actuator models widely used in robotics, including muscles [40], joint torques [75], and pneumatic actuators [12]. With the support of the articulated skeleton constraints, dry frictional contact, and versatile actuators, our novel differentiable algorithm can simulate soft articulated robots and compute gradients for a wide range of applications.
The key contributions of this work are as follows.
• A top-down matrix assembly algorithm within Projective Dynamics to make soft-body dynamics compatible with reduced-coordinate articulated systems (Sec. 4).
• An extended and generalized dry friction model for soft solids with a new matrix splitting strategy to stabilize the solver (Sec. 5).
• Analytical models of muscles, joint torques, and pneumatic actuators to enable more realistic and stable simulation results (Sec. 4.3 and Appendix C & D).
• A unified differentiable framework that incorporates skeletons, contact, and actuators to enable gradient computation for learning and optimization (Sec. 6).
• Experimental validation demonstrating that differentiable physics accelerates system identification and motion control with soft articulated bodies up to orders of magnitude (Sec. 6).
In the following paper, Section 2 will discuss related papers in deformable body simulation and differentiable physics. Section 3 is basically a preliminary of projective dynamics to explain the high-level simulation framework and define notations that will be used later. Section 4 and 5 are about how to deal with articulated skeleton and contact in our method. Section 6 will show the ablation studies and comparisons results with other learning methods. Code is available on our project page: https://github.com/YilingQiao/diff_fem

heading: 2 Related Work
text: Deformable body simulation using Finite Element Method (FEM) plays an important role in many scientific and engineering problems [32, 20, 57]. Previous works model soft bodies using different representations and methods for specific tasks. There are several kinds of approaches for modeling body actuation. Pneumatic-based methods [9] change the rest shape to produce reaction forces. Rigid bones attached within soft materials are also used to control the motion of deformable bodies [52, 38, 18, 45]. To further simulate biologically realistic motion, it is common to apply joint torques in articulated skeletons [36]. For example, [33, 76] use articulated body dynamics to govern the motion while handling collisions using soft contact. Inspired by animals, different designs of muscle-like actuators for soft-body simulations were also proposed [41, 1, 42].
Regarding contact modeling, spring-based penalty forces are widely used [58, 67, 27] for their simplicity. More advanced algorithms include inelastic projection [6, 24] and barrier-based repulsion [47, 48]. However, these methods do not always conform to Coulomb’s frictional law. We opt for a more realistic dry frictional model [44] to better handle collisions.
Projective Dynamics [5] is widely used for its robustness and efficiency for implicit time integration. It has been extended to model muscles [59], rigid skeletons [46], realistic materials [53], and accurate contact forces [54]. Our method also adopts this framework for faster and more stable time integration. In contrast to the aforementioned methods using Projective Dynamics, our algorithm is the first to enable joint actuation in articulated skeletons together with a generalized dry frictional contact for soft body dynamics.
Differentiable physics has recently been successfully applied to solve control and optimization problems. There are several types of physically-based simulations that are differentiable, including rigid bodies [10, 11], soft bodies [30, 29, 39, 19, 16], cloth [51, 62], articulated bodies [21, 74, 63], and fluids [71, 73, 28, 70]. Differentiable physics simulation can be used for system identification [68, 26], control [69], and design [14, 49]. For differentiable soft-body dynamics, Du et al. [16] propose a system for FEM simulation represented by volume mesh. This system has been applied to robot design [56] and control [15]. Different from this work [16], our approach uses tetrahedral meshes with adaptive resolution to model finer detail and scale better to complex articulated bodies.
Hu et al. [30], Krishna Murthy et al. [39] use source code transformation to differentiate the dynamics, but their contact model does not follow Coulomb’s law. Geilinger et al. [19] simulate soft materials attached to rigid parts with penalty-based contact force, but their use of maximal coordinates makes it difficult to incorporate joint torques. In comparison, our model has realistic contact handling, versatile actuators, and skeletons with joint constraints, thereby enabling our method to simulate a much wider range of soft, multi-body systems not possible before.
There are other works that approximate physical dynamics using neural networks [50, 2, 72, 65]. These methods are inherently differentiable but cannot guarantee physical correctness outside the training distribution.

heading: 3 Soft Body Simulation Using Projective Dynamics
text: We use Projective Dynamics [5] to model the physics of soft, multi-body systems because its efficient implicit time integration can make the simulation more stable. We briefly introduce Projective Dynamics below. The dynamics model can be written as
M(qn+1 − qn − hvn) = h2(∇E(qn+1) + fext), (1)
with M being the mass matrix, qn the vertex locations at frame n, h the time step, vn the velocity, E the potential energy due to deformation, and fext the external forces. We choose implict Euler for a stable time integration, then the state qn+1 can be solved by
qn+1 = arg min q
1
2h2 (q− sn)>M(q− sn) + E(q), (2)
where sn = qn + hvn + h2M−1fext. Projective Dynamics reduces the computational cost by introducing an auxiliary variable p to represent the internal energy as the Euclidean distance between p and q after a projection G:
E(q) = ∑ i ωi 2 ‖Giq− pi‖2F , (3)
where ω is a scalar weight, and E contains internal energy from different sources, such as deformations, actuators, and constraints. The computation of ω, G, and p is dependent on the form of the energy. Combining all energy components into Eq. 2, we have
qn+1 = arg min q
1 2 q> ( M h2 + L ) q + q> ( M h2 sn + Jp ) , (4)
where L = ∑ ωiG > i Gi and J = ∑ ωiG > i Si, and Si is the selector matrix. Since this is a quadratic optimization without constraints, its optimal point is given by the solution of the following linear system: (
M h2 + L
) qn+1 = M
h2 sn + Jp (5)
Note that the estimation of p is based on the current values of q. Therefore we need to alternate between computing p and solving q until convergence. Luckily, both steps are fast and easy to solve: solving q in Eq. 5 is easy because it is a simple linear system, and computing p can be fast because it is local and can be parallelized. We show the generic Projective Dynamics method in Alg. 1.
Algorithm 1 Soft body simulation using Projective Dynamics 1: x1 ← initial condition 2: for t = 1 to n− 1 do 3: while not converged do 4: Compute pi for all energy components i according to q (Local step) 5: Solve q in Eq. 5 according to pi (Global step) 6: end while 7: end for

heading: 4 Articulated Skeletons
text: Skeletons are indispensable for vertebrate animals and nowadays articulated robots. However, adding skeletons into the soft body simulation is challenging. First, the rigid bones cannot be simply replaced by soft materials with large stiffness, since this can make the system unstable and unrealistic. Second, joint connections between bones must be physically valid at all times, and thus also cannot be modeled as soft constraints. Moreover, the formulation should support joint actuation as torques to drive the multi-body system like an articulated robot. Li et al. [45] proposed a method for passive articulated soft-body simulation with ball joint constraints. We extend this method to enable rotational/prismatic joints, torque actuation, and precise joint connections without introducing extra constraints.

heading: 4.1 Rigid Body System
text: When integrated with hard skeletons, vertices on the rigid parts can be expressed as qk = QT r kVk, (6) where Q = (I 0) is the projection from homogeneous coordinates to 3D coordinates, Trk ∈ R4×4 the rigid transformation matrix, and Vk ∈ R4×mk the rest-pose homogeneous coordinates of the kth rigid body.
During the global step in Projective Dynamics, we do not directly solve for Trk, but for the increment ∆zk in its degree-of-freedom (DoF), to avoid nonlinearity. This formulation restricts the changes of the rigid vertices to the tangent space yielded by the current Trk:
qi+1k = q i k + ∆q i k ≈ qik + ∂qik ∂zk ∆zik, (7)
where qik are the vertex locations of the k th body in the ith iteration step, and zk is the variable defining the DoF of the kth rigid body, including the rotation and the translation. The nonrigid part of the vertices can also be integrated with this formulation simply with ∂q i
∂z = I.
Let B = ∂q i
∂z be the Jacobian of the concatenated variables. Eq. 4 can be rewritten as
∆zi = arg min ∆z
1 2 ∆z>B>
( M
h2 + L
) B∆z + ∆z>B> (( M
h2 + L
) qi − ( M
h2 sn + Jp
)) (8)
After solving ∆zi, the new vertex states qi+1 are derived by the new rigid transformation matrix Tr′k using Eq. 6, which is subsequently computed in the local step, discussed below.
Local step. The variable ∆zik is composed of the increment of rotation ωik and translation lik based on the current transformation Trk:
∆qik = [ −[qik]× I ] [ωik lik ] . (9)
Here [qik]× is defined as the vertical stack of the cross product matrices of all vertices in the k th rigid body. During the local step, we compute the SVD of the new transformation matrix after integrating ωik and l i k,
Tk =
[ I + ωi∗k 0
0 1
] Trk + [ 0 lik 0 0 ] = UΣV>, (10)
and restrict it to SO(3) to obtain the new rigid transformation, Tr′k = UV >.
The local step of Projective Dynamics for a single rigid body is the same as in [46]. However, we propose Eq. 7 to generalize the coupling to kinematic trees [55] with precise and actuated articulation.

heading: 4.2 Top-down Matrix Assembly for Articulated Body Systems
text: The articulated body system formulation is similar to the rigid body one, except that the transformation matrix is now chained:
Trk = ∏ u∈Uk Au, (11)
where Uk contains all ancestor of the kth link (inclusive), and Au is the local transformation matrix defined by joint u.
For rigid bodies, the vertex locations of a rigid body only depend on the body’s own DoF variables. In the articulated system, however, they are also affected by the body’s ancestors. Therefore, B changes from a block diagonal matrix to a block lower triangular matrix if the rigid body vertices are ordered by their kinematic tree depth.
To compute the matrix B, we consider a link u with one of its non-root ancestor v. By the definitions in Sec. 4.1, the corresponding block in matrix B is
Bu,v = ∂TruVu ∂zv = QPv ∂Av ∂zv Sv,uVu, (12)
where Pv is the prefix product of the local transformation matrices of the link chain from root to v (exclusive), and Sv,u is the suffix product from v to u. In boundary cases where u = v, the formulation becomes
Bu,u = QPu ∂Au ∂zu Vu. (13)
When v is the root and thus has the same DoFs as a rigid body, using the results from Sec. 4.1, the formulation can be simplified to
Bu,root = Q [−[qu]× I] . (14)
Computing Eq. 12 requires the matrix products Pv and Sv,u of a link chain (v, u) in the tree. Straightforward approaches here could result in O(N3) complexity, where N is the number of links. However, by utilizing the kinematic tree and conducting the computation in top-down order, the complexity can be reduced to O(N2), which is optimal. The key observation is that the prefix and suffix products can be computed recursively:
Pv = Pv′Av′ (15) Sv′,u = AvSv,u, (16)
assuming v′ is the parent link of v. When we traverse the kinematic tree in a depth-first order, the prefix product can be computed in O(1). The suffix product is also obtained as we iterate along the path back to the root. Algorithm 2 shows the matrix assembly method starting from the root node:
Algorithm 2 Matrix Assembly for the Articulated System 1: Input: tree link u 2: Compute Pu using Eq. 15 3: v ← u 4: while v is not root do 5: Compute Sv,u using Eq. 16 6: Compute Bu,v using Eq. 12 7: v ← parent(v) 8: end while 9: Compute Bu,root using Eq. 14 10: for s in descendants(u) do 11: Solve link s recursively 12: end for
The transformation matrix A and the Jacobian of a joint depend on the joint type. This is further derived in Appendix C.

heading: 4.3 Articulated Joint Actuation
text: Eq. 8 is a quadratic optimization, so the optimal ∆zi is given by the linear system
H∆zi = k, (17)
where H = B> ( M h2 + L ) B and k = −B> (( M h2 + L ) qi − ( M h2 sn + Jp )) . Reordering the vertices into sets of deformable and rigid ones yields the following partitioning of the matrix:[ Hd H > c
Hc Hr ] [ ∆zid ∆zir ] = [ kd kr ] , (18)
where ∗d and ∗r represents the deformable and the rigid parts, respectively. The joint actuation can be directly added to kr since the linear system is analogous to the basic formulation Ma = f where the right hand side represents the sum of forces and/or torques. The formulation of pneumatic and muscle actuators can be found in Appendix D.

heading: 5 Contact Modeling
text: We handle the contact using Coulomb’s frictional law via a Jacobian. To compute the velocities after collisions, we split the left-hand side of Equation 5 into the diagonal mass matrix M and the constraint matrix h2L, and move the latter to the right-hand side:
Mvi+1 = f − h2Lvi + ξi, (19)
where f = Msn − (M + h2L)qn + h2Jp and the contact force ξi is determined according to f − h2Lvi (the current momentum) to enforce non-penetration and static/sliding friction. The idea here is to enforce Coulomb’s law at every iteration, which is ensured by solving vi+1 using the inverse of a diagonal matrix M. As long as the solver converges at the end, the final v and ξ will conform to the frictional law.
This method works well for cloth contacts [54], but cannot be directly applied to soft bodies, because solid continuum is much stiffer than thin sheets, i.e. the elements in h2L on the right-hand side are much larger than those in M on the left-hand side, resulting in severe oscillation or even divergence during the iterative solve.
We show that in order to guarantee the convergence of Equation 19, the time step h has to satisfy a certain condition: Proposition 1. Assuming f and ξ are fixed, Equation 19 converges if the time step h satisfies
h2 < ρ
24 √ 3Tµ ∑3 k=1 ‖qk − q0‖22
(20)
where ρ is the density, µ is the stiffness, T is the number of tetrahedra, and qi are the vertex positions.
Details of the proof can be found in Appendix A. Using the setting in our experiments, where T ≈ 1000, µ ≈ 3× 105, ‖qk−q0‖2 ≈ 10−2, and ρ ≈ 1, we would need to set h < 1/1934 in order to ensure the convergence, which is too strict for the simulation to be useful in general applications.
Splitting scheme. We improve this method to be compatible with soft body dynamics by introducing a new splitting scheme. Eq. 19 is reformulated as
(M + h2D)vi+1 = f − h2(L−D)vi + ξi, (21)
where D are the diagonals of L. Our key observation is that (a) the diagonals of L are necessary and sufficient to stabilize the Jacobian iteration, and (b) adding extra diagonal elements to the left-hand side will not break the Coulomb friction law. We show in Appendix B that under the same assumption as Proposition 1, our method is guaranteed to converge no matter how big h is. This improvement accelerates the simulation since larger time steps mean faster computation.
We also note that the new splitting scheme will not modify the behavior of the collision response because the convergence point of Eq. 19 is the same as that of Eq. 21, and thus Coulomb’s Law is still satisfied at convergence.

heading: 6 Experiments
text: In this section, we first introduce our implementation and then report ablation studies that demonstrate the importance of skeletons and collision contacts in soft-body dynamics. Subsequently, we use the
gradients computed by our method to perform system identification; specifically, we estimate the physical parameters of bridges. Finally, we perform gradient-based learning of grasping and motion planning on robots with various actuators, including a pneumatic gripper, an octopus with muscles, and a skeleton-driven fish. Our method can converge more than an order of magnitude faster than reinforcement learning and derivative-free baselines.

heading: 6.1 Implementation
text: Our simulator is written in C++, the learning algorithms are implemented in PyTorch [60], and Pybind [34] is used as the interface. We run our experiments on two desktops, one with an Intel Xeon W-2123 CPU @ 3.6GHz and the other with an Intel i9-10980XE @ 3.0GHz, respectively. For differentiation, the numerical data structure in our simulator is templatized and integrated into the C++ Eigen library, such that our method can conveniently interoperate with autodiff tools to differentiate the dynamics. Our method can also run in pure C++ to perform forward simulation. We refer to the open-source code from [59] (Apache-2.0), [54] (GNU GPL v3.0), and [45] (MPL2) during our implementation. More details can be found in our code in the supplement.
To further improve the memory efficiency, we introduce a checkpointing scheme [7] into our pipeline. Instead of storing the entire simulation history, we only store the system’s state in each step. During the backward pass, we reload the saved state vector and resume all the intermediate variables before the backpropagation. This strategy can save a major part of the memory, compared to the brute-force implementation. We conduct an experiment to compare the memory consumption with and without this check-
pointing scheme. The results are reported in Table 1. CppAD [3] is used to differentiate the simulation here. In this experiment, we simulate a bridge and estimate its material properties as shown in Figure 3(a). The results show that the memory footprint of the baseline scales linearly with simulation length, while our checkpointing scheme keeps memory consumption nearly constant.

heading: 6.2 Ablation Study
text: Skeleton constraints. Controlling soft characters via skeletons is natural and convenient: vertebrate animals are soft, but are driven by piecewise-rigid skeletons. Our simulator supports skeletons and joint torques within soft bodies. This ablation study compares other designs with ours. In this experiment, a Baymax model [13] in its T-pose is released from above the ground, as shown in Figure 1. We embed 5 bones inside Baymax (4 in arms and legs, and 1 in the torso). When Baymax falls to the ground, we also add torques on its shoulders so it can lift its arms to a target Y-pose. More details of the setting and qualitative results can be found in Appendix E and the supplementary video. Three metrics, summarized in Figure 1, are used to measure realism. The metrics are averaged over 5 repetitions with different initial positions and velocities. For comparison, we simulate a ‘No skeleton’ Baymax without the support of rigid bones. Its bone error is non-zero because of the deformation. The Baymax in a differentiable rigid body simulator [62] is rigid, so the body length error is non-zero. Li et al. [45] simulate the ‘Passive’ skeleton case where there is no joint actuation and joint angles cannot be adjusted to the desired configuration. We also run Difftaichi-MPM [30] by converting the mesh model to the point-based MPM representation. ‘MPM’ does not have skeletons so the errors are high. The arms also detach from the body so the joint error is NaN. Our method attains the highest degree of physical realism and correctness overall.
Contact handling. Good contact handling is critical for simulating multi-body systems that interact with their environment. In this experiment, we throw a 3D soft ball against a 2D thin sheet. Metrics in this experiment are penetration error and indicators of vertical compression and horizontal stretching. Zero penetration error is ideal. ‘Yes’ for compression/stretching indicates that the simulator can model the deformation of the soft ball correctly. The metrics are averaged over 5 experiments with different initial positions and velocities. The dry frictional contact model of Ly et al. [54] does not model the deformation of soft solids, and there could be penetration when the resolutions of the ball and cloth differ a lot due to the nodal collision handling scheme. The rigid differentiable simulator of Qiao et al. [62] can prevent interpenetration, but the ball remains rigid. MPM [30] can model the deformation of both the ball and the cloth, but the cloth is torn apart by the ball and penetration cannot be quantified. In contrast, our method accurately handles collision to avoid interpenetration and correctly simulates the deformation of the ball.

heading: 6.3 Applications
text: System identification. Determining the material parameters of deformable objects can be challenging given their high dimensionality and complex dynamics. In this experiment, we use our differentiable simulator to identify the material property of each finite element cell within the soft body. As shown in Figure 3, there are two bridges with unknown materials: a suspension bridge with both ends fixed and the entire bridge being soft, and an arch bridge that has three piers attached to the ground. Given that the movement of the barycenter under gravity, compared to its rest pose, is ∆x = 8cm, we estimate Young’s modulus and Poisson’s ratio of each finite element cell in the bridge. The loss function is the distance from the actual barycenter to the target. The suspension bridge has n = 668 cells and the arch one has n = 2911 cells. The number of unknowns is 2n. We compare our method with four derivative-free methods (CMA-ES [25], LEAP [8], BOBYQA [61], and Nelder-Mead [66]). Each experiment is repeated 5 times with different random seeds. As shown in the figure, our method converges in ∼10 iterations while others fail to converge even after 100 iterations, indicating that derivative-free methods in this high-dimensional setting become too inefficient to converge to a reasonable solution. By making use of the gradients provided by our method, common gradient-based algorithms can quickly reach the target configuration.
Motion planning. Controlling the motion of deformable bodies is challenging due to their flexible shapes. In this experiment, summarized in Figure 4, the task is to control robots with different actuator types. In general, given an initial state X0, control policy φθ1(·), and material parameters θ2, the simulator can generate a trajectory of the states at all time steps t, {simt(X0, φθ1(·), θ2)} = {X1, X2, ..., Xn}. If we want the system to reach a target state Xtarget at the end of the simulation, we can define an objective function L(X0, θ1, θ2) = ||(simN (X0, φθ1(·), θ2)−Xtarget)||2, where the optimization variable are θ1 and/or θ2.
We compare our method with Reinforcement Learning algorithms (SAC [23], SQL [22], and PPO [64]), and the best derivative-free optimization method from the last experiment, CMA-ES. We also tried MBPO [35], but we found that this method takes too much memory and could not finish any test. All RL methods use the negative of the loss as the reward.
The pneumatic gripper in Figure 4(a) has 56 pneumatic cells in four arms and is attached to an (invisible) drone as in [17]. The pneumatic activation can control the volume of a tetrahedron. When the cells inflate, the arms will move inwards and hold the ball tighter. We control the pneumatic activation as well as the movement of the drone to move the ball from the start (0, 0, 0) to our target (0, 0.3, 0) in 50 steps. The loss is the distance from the actual position to the target position. Our method converges in 10 episodes while CMA-ES and PPO gradually converge in 200 and 500 episodes, respectively.
The muscle-driven octopus in Figure 4(b) has 8 legs, each with 2 muscles inside. It moves forward by actuating the muscles, being pushed by drag and thrust forces induced by the water on the octopus’s surface [59]. The octopus starts at (0, 0, 0) and our target location is (−0.4, 0.8,−0.4). We set the objective to be the distance between the current location and the desired location. The length of the simulation is 400 steps, and the control input in each step is 64-dimensional. In total, there are 64× 400 = 25600 variables to optimize. Our method converges in 50 episodes while other methods fail to converge in 500 episodes.
The fish with an embedded skeleton in Figure 4(c) has 6 bones: 3 in its body, 2 in the fins, and 1 in the tail. The hydrodynamics in this environment is the same as in the octopus experiment. The fish starts at (0, 0, 0) in step 1 and the target location in step 100 is (0, 0, 0.15). The objective function is the distance from the actual location to the target location. For each step, there will be a torque vector of size 5 that represents the joint actuation level. In total, the optimization variable has 500 dimensions. Our method with gradient-based optimization can converge in roughly 50 episodes, while others cannot converge even after 500 episodes.
In summary, gradient-free optimization methods and RL algorithms meet substantial difficulties when tackling problems with high dimensionality, such as soft, multi-body systems. Even when the action space is as small as the one in the gripper case, RL methods still fail to rapidly optimize the policy. By introducing the gradients of the simulation, simple gradient-based optimization outperforms other algorithms. This work hopefully may inspire improvements in RL algorithms that tackle such high-dimensional problems.

heading: 7 Conclusion
text: Our paper has developed a differentiable physics framework for soft, articulated bodies with dry frictional contact. To make the simulation realistic and easy-to-use, we designed a recursive matrix
assembly algorithm and a generalized dry frictional model for soft continuum with a new matrix splitting strategy. Integrated with joint, muscle, and pneumatic actuators, our method can simulate a variety of soft robots. Using our differentiable physics to enable gradient-based optimization, our method converges more than an order of magnitude faster than the baselines and other existing alternatives.
There are some limitations in our contact handling and soft body dynamics. Currently, though our algorithm is more extensive and generalized than existing differentiable physics algorithms and our implementation handles the most commonly found contact configuration, vertex-face collisions, there could still be edge-edge penetration missed in some corner cases. Moreover, the Projective Dynamics pipeline limits the energy to have the form E = ‖Gq− p‖. Some nonlinear material models (e.g., neo-Hookean) are not captured in this framework and new models for differentiable physics will be required for handling nonlinear and heterogeneous materials. For future work, we aim to add edge-edge collision handling in the Projective Dynamics pipeline. The techniques in [53] can be used to incorporate addition material types. GPU or other parallel computing implementation can be used to boost the performance of gradient computation.
Acknowledgements. This research is supported in part by Army Research Office, National Science Foundation, Dr. Barry Mersky and Capital One Endowed Professorship.

Review summaries:
This paper proposes a soft multi-body differentiable physics simulation framework based on projective dynamics. The authors propose a top-down matrix assembly algorithm for rigid body simulation algorithms and a new matrix splitting strategy for a generalized dry friction model for soft continuum. The experiments demonstrate that their designs make soft body simulation more stable and realistic compared to other frameworks, and the gradients help accelerate system identification and motion control of soft robots.

This paper takes on the timely problem of devising differentiable simulators for soft bodies, with a focus on articulated bodies where the soft nature of the material comes with structure in the way different parts are connected together - hence requiring formulation of equations of motion in the sense of multi-body dynamics.  The core contribution is to formulate the simulation problem in a projective dynamics framework, and to exploit the implicit matrix structure in the equations of motion to achieve computational speedups.   This has been implemented in C++ libraries and demonstrated in a few different simulation settings, including a Baymax soft articulated body, a deformable ball falling on cloth and a bridge -  each representing a different form of structure. The main outcome is that this form of differentiable simulation enables faster system identification and control.

The paper presents a new soft-body physics simulator that is end-to-end differentiable, and shows how it supports efficient system identification and control.  Experiments show the value of the simulator for a range of complex realistic situations -- although only in simulation.  The experiments also report ablations of key model components (skeletons, collision contacts). 

The proposed paper introduces a novel method for solving the simulation of soft articulated bodies based on a differentiable solver. The key contribution of the introduced method are a top-down matrix assembly approach for projective dynamics, a generalized friction model, analytic models for muscles, and a unified framework that enables computing gradients in a differentiable manner. The method is evaluated based on a number of meaningful experiments and ablation studies.