ChainQueen: A Real-Time Differentiable Physical Simulator for Soft Robotics

I Introduction

Robot planning and control algorithms often rely on physical simulators for prediction and optimization . In particular, differentiable physical simulators enable the use of gradient-based optimizers, significantly improving control efficiency and precision. Motivated by this, there has been extensive research on differentiable rigid body simulators, using approximate and exact methods.

Significant challenges remain for deformable objects. First, simulating the motion of deformable objects is slow, because they have much higher degrees of freedom (DoFs). Second, contact detection and resolution is challenging for deformable objects, due to their changing geometries and potential self-collisions. Third, closed-form and efficient computation of gradients is challenging in the presence of contact. As a consequence, current simulation methods for soft objects cannot be effectively used for solving inverse problems such as optimal control and motion planning.

In this paper, we introduce a real-time, differentiable physical simulator for deformable objects, building upon the Moving Least Squares Material Point Method (MLS-MPM) . We name our simulator ChainQueen Or 乾坤 , literally “everything between the sky and the earth.”. The Material Point Method (MPM) is a hybrid Lagrangian-Eulerian method that uses both particles and grid nodes for simulation . MLS-MPM accelerates and simplifies traditional MPM using a moving least squares force discretization. In ChainQueen, we introduce the first fully differentiable MLS-MPM simulator with respect to both state and model parameters, with both forward simulation and back-propagation running efficiently on GPUs. We demonstrate the ability to efficiently calculate gradients with respect to the entire simulation. This enables many novel applications for soft robotics including optimization-based closed-loop controller design, trajectory optimization, and co-design of robot geometry, materials, and control.

As a particle-grid-based hybrid simulator, MPM simulates objects of various states, such as liquid (e.g., water), granular materials (e.g., sand), and elastoplastic materials (e.g., snow and human tissue). ChainQueen focuses on elastic materials for soft robotics. It is fully differentiable and $4-9\times$ faster than the current state-of-the-art. Numerical and experimental validation suggest that ChainQueen achieves high precision in both forward simulation and backward gradient computation.

ChainQueen’s differentiability allows it to support gradient-based optimization for control and system identification. By performing gradient descent on controller parameters, our simulator is capable of solving these inverse problems on a diverse set of complex tasks, such as optimizing a 3D soft walker controller given an objective. Similarly, gradient descent on physical design parameters, enables inference of physical properties (e.g. mass, density and Young’s modulus) of objects and optimizing design for a desired task.

In addition to benchmarking ChainQueen’s performance and demonstrating its capabilities on a diverse set of inverse problems, we have interfaced our simulator with high-level python scripts to make ChainQueen user-friendly. Users at all levels will be able to develop their own soft robotics systems using our simulator, without the need to understand its low-level details. We will open-source our code and data and we hope they can benefit the robotics community.

II Related Work

The material point method has been extensively developed from both a solid mechanics and computer graphics perspective. As a hybrid Eulerian-Langrangian method, MPM has demonstrated its versatility in simulating snow , sand , non-Newtonion fluids , cloth , solid-fluid coupling , rigid body coupling, and cutting . also proposed an adaptive MPM scheme to concentrate computation resources in the regions of interest.

There are many benefits of using MPM for soft robotics. First, MPM is a well-founded and physically-accurate discretization method and can be derived through the weak form of conservation laws. Such a physically-based approach makes it easier to match simulation with real-world experiments. Second, MPM is friendly to parallelization on modern hardware architectures. Closely related to our work is a high-performance GPU implementation by Gao et al., from which we borrow many useful optimization practices. Though efficient when solving forward simulation, their simulator is not differentiable, making it inefficient for inverse problems in robotics and learning. Third, MPM naturally handles large deformation and (self-)collision, which are common in soft robotics, but often not modeled in, e.g., mesh-based approaches due to computational expense. Finally, the continuum dynamics (including soft object collision) are governed by the smooth (and differentiable) potential energy, making the whole system differentiable.

Our simulator, ChainQueen, is fully differentiable and the first simulator that applies MPM to soft robotics.

II-B Differentiable Simulation and Control

Recently, there has been an increasing interest in building differentiable simulators for planning and control. For rigid bodies, , and proposed to approximate object interaction with neural nets; later, explored their usage in control. Approximate analytic differentiable rigid body simulators have also been proposed . Such systems have been deployed for manipulation and planning .

Differentiable simulators for deformable objects have been less studied. Recently, proposed SPNets for differentiable simulation of position-based fluids . The particle interactions are coded as neural network operations and differentiability is achieved via automatic differentiation in PyTorch. A hierarchical particle-based object representation using neural networks is also proposed in . Instead of approximating physics using neural networks, ChainQueen differentiates MLS-MPM, a well physically founded discretization scheme derived from continuum mechanics. In summary, our simulator can be used for a more diverse set of objects; it is more physically plausible, and runs faster.

III Forward simulation and back-propagation

We use the moving least squares material point method (MLS-MPM) to discretize continuum mechanics, which is governed by the following two equations:

Particle-to-grid transfer (P2G). Particles transfer mass $m_{p}$, momentum $(m\mathbf{v})_{p}^{n}$, and stress-contributed impulse to their neighbouring grid nodes, using the Affine Particle-in-Cell method (APIC) and moving least squares force discretization , weighted by a compact B-spline kernel $N$:

Grid operations. Grid momentum is normalized into grid velocity after division by grid mass:

Note that neighbouring particles interact with each other through their shared grid nodes, and collisions are handled automatically. Here we omit boundary conditions and gravity for simplicity.

Grid-to-particle transfer (G2P). Particles gather updated velocity $\mathbf{v}_{p}^{n+1}$, local velocity field gradients $\mathbf{C}_{p}^{n+1}$ and position $\mathbf{x}_{p}^{n+1}$. The constitutive model properties (e.g. deformation gradients $\mathbf{F}_{p}^{n+1}$) are updated.

For soft robotics, we additionally introduce an actuation model. Inspired by actuators such as , we designed an actuation model that expands or stretches particle $p$ via an additional Cauchy stress $\mathbf{A}_{p}=\mathbf{F}_{p}\bm{\sigma}_{pa}\mathbf{F}_{p}^{T}$, with $\bm{\sigma}_{pa}=\text{Diag}(a_{x},a_{y},a_{z})$ – the stress in the material space. This framework supports the use of other differentiable actuation models including pneumatic, hydraulic, and cable-driven actuators. Fig. 1 illustrates forward simulation and back propagation.

MLS-MPM is naturally differentiable. Though the forward direction has been extensively used in computer graphics, the backward direction (differentiation or back-propagation) is largely unexplored.

Based on the gradients we have derived analytically, we have designed a high-performance implementation that resembles the traditional forward MPM cycle: backward P2G (scatter particle gradients to grid), grid operations, and backward G2P (gather grid gradients to particles). Please see the supplemental document for the gradient derivations. Gradients of state at the end of a time step with respect to states at the starting of the time step can be computed using the chain rule. With the single-step gradients computed, applying the chain rule at a higher level from the final state all-the-way to the initial state yields gradients of the final state with respect to the initial state, as well as to the controller parameters that are used in each state. We cache all the simulation states in memory, using a “memo” object. Though the underlying differentiation is complicated, we have designed a simple high-level TensorFlow interface on which end-users can build their applications (Fig. 2).

Our high-performance implementationBased the Taichi open source computer graphics library. takes advantage of the computational power of modern GPUs through CUDA. We also implemented a reference implementation in TensorFlow. Note that programming physical simulation as a “computation graph” using high-level frameworks such as TensorFlow is less inefficient. In fact, when all the overheads are gone, our optimized CUDA solver is $132\times$ faster than the TensorFlow reference version. This is because TensorFlow is optimized towards deep learning applications where data granularity is much larger and memory access pattern is much more regular than physical simulation, and limited CPU-GPU bandwidth. In contrast, our CUDA implementation is tailored for MLS-MPM and explicitly optimized for parallelism and locality, thus delivering high-performance.

IV Evaluation

In this section, we conduct a comprehensive study of the efficiency and accuracy of our system, in both 2D and 3D.

Instead of using complex geometries, a simple falling cube is used for performance benchmarking, to ensure easy analysis and reproducibility. We benchmark the performance of our CUDA simulator against NVIDIA Flex , a popular PBD physical simulator capable of simulating deformable objects. Note that both PBD and MLS-MPM needs substepping iterations to ensure high stiffness. To ensure fair comparison, we set a Young’s modulus, Poisson’s ration and density so that visually ChainQueen gives similar results to Flex. We used two steps per frame and four iterations per step in Flex. Note that setting exactly the same parameters is not possible since in PBD there is no explicitly defined physical quantity such as Young’s modulus.

We summarize the quantitative performance in Table II. Our CUDA simulator provides higher speed than Flex, when the number of particles are the same. It is also worth noting that the TensorFlow implementation is much slower, due to excessive runtime overheads.

IV-B Accuracy

We design five test cases to evaluate the accuracy of both forward simulation and backward gradient evaluation:

A1 (analytic, 3D, float32 precision): final position w.r.t. initial velocity (with collision). This case tests conservation of momentum, gradient accuracy and stability of back-propagation.

A2 (analytic, 3D, float32 precision): same as A1 but with one collision to a friction-less wall.

B (numeric, 2D, float64 precision): colliding billiards. This case tests gradient accuracy and stability in more complex cases where analytic solutions do not exist. We used float64 precision for accurate finite difference results.

C (numeric, 2D, float64 precision): finger controller. This case tests gradient accuracy of controller parameters, which are used repeatedly in the whole simulation process.

D1 (experimental, pneumatic actuator, actuation) In order to evaluate our simulator’s real-world accuracy, we compared the deformation of a physical actuator to a virtual one. The physical actuator has four pneumatic chambers which can be inflated with an external pump, arranged in a cross-shape. Inflating the individual chambers bends the actuator away from that chamber. The actuator was cast using Smooth-On Dragon Skin 30.

D2 (experimental, pneumatic actuator, bouncing) In a second test, we dropped the same actuator from a 15 cm height, and compared its dynamic motion to a simulation.

In 3D analytic test cases, where gradients w.r.t. initial velocity can be directly evaluated as in Table III. For the experimental comparisons, the results are shown in Fig. 3. In addition to our simulator’s high performance and accuracy, it is worth noting that that the gradients remain stable in the long term, within up to $1000$ time steps.

V Inference, Control and Co-design

The most attractive feature of our simulator is the existence of quickly computable gradients, which allows the use of much more efficient gradient-based optimization algorithms. In this section, we show the effectiveness of our differentiable simulator on gradient-based optimization tasks, including physical inference, control for soft robotics, and co-design of robotic arms.

ChainQueen can be used to infer physical system properties given its observed motion, e.g. perform gradient descent to infer the relative densities of two colliding elastic balls (see figure above, ball A moving to the right hitting ball B, and ball B arrives the destination C). Gradient-based optimization infers that relative density of ball A is $2.26$, which constitutes to the correct momentum to push B to C. Such capability makes it useful for real-world robotic tasks such as system identification.

V-B Control

We can optimize regression-based controllers for soft robots and efficiently discover stable gaits. The controller takes as input the state vector $\mathbf{z}$, which includes target position, the center of mass position, and velocity of each composed soft component. In our examples, the actuation vector $\mathbf{a}$ for up to $16$ actuators is generated by the controller in each time step. During optimization, we perform gradient descent on variables $\mathbf{W}$ and $\mathbf{b}$, where $\mathbf{a}=\tanh{(\mathbf{W}\mathbf{z}+\mathbf{b})}$ is the actuation-generating controller.

We have designed a series of experiments including the 2D biped runner (Fig. 4) and robotic finger, and 3D quadrupedal runner (Fig. 6), crawler and robotic arm. Gradient-based optimizers successfully compute desired controllers within only tens or hundreds of iterations. Visual results are included in the supplemental video.

To emphasize the merits of gradient-based approaches, we compare our control method with proximal policy optimization (PPO) , a state-of-the-art reinforcement learning algorithm. PPO is an actor-critic method which relies on sampled gradients of a reward function in order to optimize a policy. This sampling-based approach is model-free; it relies on gradients of the rewards with respect to controller parameters, but not with respect to the physical model for updates. For our comparison, we use velocity projected onto the direction toward the goal as the reward. Note that this is functionally extremely similar to a distance loss; the cumulative reward $\int_{t}{=0}^{T}{v_{goal}dt}=D-\|x_{T}-x_{goal}\|$, where $D$ is the initial distance and $x_{T}$ and $x_{goal}$ represent world coordinates of the robot at time $T$ and of the goal, respectively. As velocity toward the goal increases, final distance to the goal decreases. We use a simplified single link version (with only two adjacent actuators) of Fig. 5 and the 2D runner Fig. 4 as a benchmark. Quantitative results for the finger are shown in Fig. 7. We performed a similar comparison on the 2D walker, the controller optimized by ChainQueen for the 2D walker starts functioning well within $20$ minutes; by comparison the policy recovered by PPO still chose nearly-random actions after over $4$ hours of training; demonstrating that for certain soft locomotion tasks our gradient-based method can be more efficient than model-free approaches.

V-C Co-design

Our simulator is capable of not only providing gradients with respect to dynamics and controller parameters, but also with respect to structural design parameters, enabling co-design of soft robots. To demonstrate this, we designed a multi-link robot arm (two links, two joints each with two side-by-side actuators; all parts deformable). Similar to shooting method trajectory optimization, actuation for each time step is solved for, along with the time-invariant Young’s modulus of the system for each particle. In our task, we optimized the end-effector of the arm to reach a goal ball with final arm velocity, and minimized for actuation cost $\sum_{i=0}^{N}u_{i}^{T}u_{i}dt$, where $u_{i}$ is the actuation vector at timestep $i$, and $N$ is the total number of timesteps. This is a dynamic task and the target pose cannot be reached in a static equilibrium. NLOPT’s sequential least squares programming algorithm was used for optimization . We compared our co-design solution to fixed designs. The designed stiffness distribution is shown in Fig. 5, along with controls. The convergence for the different tasks can be seen in Fig. 8. As can be seen, only the co-design arm fully converges to the target goal, and with lower actuation cost. Actuation for each chamber was clamped, and rnges of 30% to 400% of a dimensionless initial Young’s modulus were allowed and chosen large enough such as to require a swing instead of a simple bend.

VI Discussion

We have presented ChainQueen, a differentiable simulator for soft robotics, and demonstrated how it can be deployed for inference, control, and co-design. ChainQueen has the potential to accelerate the development of soft robots. We have also developed a high-performance GPU implementation for ChainQueen, which we plan to open source.

One interesting future direction is to couple our soft object simulation with rigid body simulation, as done in . As derived in , the $\Delta t$ limit for explicit time integration is $C\Delta x\sqrt{\frac{\rho}{E}}$, where $C$ is a constant close to one, $\rho$ is the density, and $E$ is the Young’s modulus. That means for very stiff materials (e.g., rigid bodies), only a very restrictive $\Delta t$ can be used. However, a rigid body simulator should probably be employed in the realm of nearly-rigid objects and coupled with our deformable body simulator. Combining our simulator with existing rigid-body simulators using Compatible Particle-in-Cell can be an interesting direction.

Acknowledgments

We would like to thank Chenfanfu Jiang, Ming Gao and Kui Wu for the insightful discussions.

In this document, we discuss the detailed steps for backward gradient computation in ChainQueen, i.e. the differentiable Moving Least Squares Material Point Method (MLS-MPM) . Again, we summarize the notations in Table IV. We assume fixed particle mass $m_{p}$, volume $V_{p}^{0}$, hyperelastic constitutive model (with potential energy $\psi_{p}$ or Young’s modulus $E_{p}$ and Poisson’s ratio $\nu_{p}$) for simplicity.

VII Variable dependencies

The MLS-MPM time stepping is defined as follows:

The forward variable dependency is as follows:

During back-propagation, we have the following reversed variable dependency:

We reverse swap two sides of the equations for easier differentiation derivation:

In the following sections, we derive detailed gradient relationships, in the order of actual gradient computation. The frictional boundary condition gradients are postponed to the end since it is less central, though during computation it belongs to grid operations. Back-propagation in ChainQueen is essentially a reversed process of forward simulation. The computation has three steps, backward particle to grid (P2G), backward grid operations, and backward grid to particle (G2P).

VIII Backward Particle to Grid (P2G)

(A, P2G) For $\mathbf{v}_{p}^{n+1}$, we have

(B, P2G) For $\mathbf{C}_{p}^{n+1}$, we have

Note, the above two gradients should also include the contributions of $\frac{\partial L}{\partial{{\mathbf{v}_{p}^{n}}}}$ and $\frac{\partial L}{\partial{{\mathbf{C}_{p}^{n}}}}$ respectively, with $n$ being the next time step.

(C, P2G) For $\mathbf{v}_{i}^{n}$, we have

IX Backward Grid Operations

(D, grid) For $\mathbf{p}_{i}^{n}$, we have

X Backward Grid to Particle (G2P)

(F, G2P) For $\mathbf{v}_{p}^{n}$, we have

(G, G2P) For $\mathbf{P}_{p}^{n}$, we have

(H, G2P) For $\mathbf{F}_{p}^{n}$, we have

(I, G2P) For $\mathbf{C}_{p}^{n}$, we have

(J, G2P) For $\mathbf{x}_{p}^{n}$, we have

(K, G2P) For $\bm{\sigma}_{pa}^{n}$, we have

XI Friction Projection Gradients

(L, grid) For $\mathbf{v}_{i}^{n}$, we have

References