Year: 2021

Abstract:
Differentiable physics enables efficient gradient-based optimizations of neural network (NN) controllers. However, existing work typically only delivers NN controllers with limited capability and generalizability. We present a practical learning framework that outputs unified NN controllers capable of tasks with significantly improved complexity and diversity. To systematically improve training robustness and efficiency, we investigated a suite of improvements over the baseline approach, including periodic activation functions, and tailored loss functions. In addition, we find our adoption of batching and a modified Adam optimizer effective in training complex locomotion tasks. We evaluate our framework on differentiable mass-spring and material point method (MPM) simulations, with challenging locomotion tasks and multiple robot designs. Experiments show that our learning framework, based on differentiable physics, delivers better results than reinforcement learning and converges much faster. We demonstrate that users can interactively control soft robot locomotion and switch among multiple goals with specified velocity, height, and direction instructions using a unified NN controller trained in our system.

heading: 1 INTRODUCTION
text: Differentiable physical simulators deliver accurate analytical gradients of physical simulations, opening up a promising stage for efficient neural network (NN) controller training via gradient descent. Existing research demonstrates that on simple tasks, learning systems via differentiable physics can effectively leverage simulation gradient information and converge orders of magnitude faster (see, e.g., de Avila Belbute-Peres et al. (2018); Hu et al. (2019)) than reinforcement learning.
However, the capability of existing differentiable physics based learning systems is relatively limited. Typically, optimized controllers can only achieve relatively simple single-goal tasks (e.g., moving in one specific direction, as in ChainQueen Hu et al. (2019)). Those learned controllers often have difficulty generalizing the task to a perturbed version.
In this work, we propose a learning framework for complex locomotion skill learning via differentiable physics. The complexity of our tasks comes from two aspects: first, the degrees of freedom of our soft agents are significantly higher compared with the rigid-body ones; second, the agents are expected to learn multiple skills at the same time. We systematically propose a suite of enhancements (Fig. 4) over existing training approaches (such as Hu et al. (2020); Huang et al. (2021)) to
significantly improve the efficiency, robustness, and generalizability of learning systems based on differentiable physics.
In addition, we investigated the contributions of each enhancements to training in detail by a series of ablation studies. We also evaluated our framework through comparisons against Proximal Policy Optimization (PPO) Schulman et al. (2017), a state-of-the-art reinforcement learning algorithm. Results show that our system is simple yet effective and has a much-improved convergence rate compared to PPO.
We demonstrate the versatility of our framework via various physical simulation environments (mass-spring systems and the moving least squares material point method, MLS-MPM Hu et al. (2018)), robot designs (structured and irregular), locomotion tasks (moving, jumping, turning around, all at a controllable velocity). To the best of our knowledge, this is the first time when a neural network can be trained via differentiable physics to achieve tasks of such complexity (Fig. 1). Our agent is simultaneously trained with multiple goals, the learned skills can be continuously interpolated. For example, we can control both the velocity and height of an agent while it runs. The trained agent can also be controlled in real-time, enabling direct applications in soft robot control and video games. In summary, our key contributions are listed below:
• To the best of our knowledge, we proposed the first differentiable physics based learning framework for multiple locomotion skills learning on soft agents using a single network.
• We developed an end-to-end differentiable physical simulation environment for deformable robot locomotion skills learning, which supports mass-spring system and material point method as dynamic backends.
• We systematically investigated the key factors contributing to the differentiable physics based learning framework. We believe our investigation can inspire further researches to explore the possibilities of differentiable physics for more complex learning tasks.

heading: 2 DIFFERENTIABLE SIMULATION ENVIRONMENTS

heading: 2.1 SIMULATION SETTINGS
text: Our learning framework supports multiple types of differentiable physically-based deformable body simulators. The simulator can be seen as a black box that outputs the states for the next time step given the actuation signal and states from the current time step. The NN controller takes the information given by the simulator as part of its input features and decides the actuation signal for the simulator in the next time step. With our automatic differentiation system, we can evaluate the gradient information to guide the training of our NN controllers. We use the mass-spring systems and the moving least squares material point method (MLS-MPM) as our simulation environments. The design choices for each environment can be found in the appendix.

heading: 2.2 AGENT STRUCTURE AND ACTUATION MODE
text: Our framework supports a variety of differently shaped agents as shown in Fig. 2 and 3. Our actuation signal exerts forces along with the ”muscle” directions of the simulated agents. This signal falls within the interval [−1,+1], where its sign determines whether an agent wants to contract (negative sign) or relax (positive sign) a muscle and its absolute value determines the magnitude of the force. In mass-spring systems, the muscle directions are represented by the spring directions, we change the rest-length of springs to generate forces. The actuation signal is used to scale the rest-length up
to some limit (usually at 20%). In addition, only activated springs, i.e. actuators (springs marked red or blue shown in Fig. 1), are able to generate forces according to control signals. In material point methods, the muscle directions are always along the vertical axis in material space. In our implementation, we simply modify the Cauchy stress in the vertical direction of material space to apply this force. Again, the force is scaled to a user-defined bound which can be adjusted for different simulations.

heading: 3 LEARNING FRAMEWORK
text: In this section, we describe our training framework based on differentiable physics (Fig. 4). In summary, we develop a differentiable simulator with an embedded NN controller. In each time step, the program performs an NN controller inference and a differentiable simulator time integration. After a few hundred time steps, the program back-propagates gradients end-to-end via reverse-mode automatic differentiation, and updates the weights of the neural network using the Adam optimizer Kingma & Ba (2014).

heading: 3.1 TASK REPRESENTATION
text: Given an agent A, its position and corresponding velocity are denoted as x = (x, y, z) and v = (u, v, w). We presents three tasks (running, jumping and crawling) for 2D agents, two (running and rotating) for 3D. In a whole simulation with total steps T , an agent is expected to achieve all goals in a goal sequence G = {G1,G2...Gn|nP = T }, where P is the time period for an agent to achieve a goal, n is the number of goals . For each goal G, there are multiple tasks in it. Each task is encoded using a target value g, which instructs the controller to drive the agent. For example, G = {gv, gh} represent performing running and jumping simultaneously, the target velocity and height are gv and gh. The agent is expected to switch between multiple goals during one simulation.
Running. The running task is defined as a velocity tracking problem. Given a time step t, the center of mass c is used to represent the position of the agent, its velocity is defined as the averaged velocity ṽ over a time period Pr, (Pr < P). Then the agent is expected to achieve a velocity ṽ close to the target velocity gv .
ṽ(t) = 1
Pr (c(t)− c(t− Pr)) (1)
Jumping. The jumping task is defined as reaching a given height gh in a time period P . Due to the gravity, it is not possible for a running robot to stay away from the ground. Therefore, the jumping height h̃h(t) of an agent is defined by the maximum vertical position of its lowest point in a certain period P . To be more formally, given S a set of nodes (or points) that constitute the agent, h is a function that can extract the vertical position given a node (or point), the jumping height over the n-th period is defined as
h̃h(n) = max t∈[0,P] min s∈S
h(s, t+ nP) (2)
where n is the index of the period in a whole simulation.
Crawling. The crawling task is defined as lowering the highest point as much as possible. The target of crawling gc is defined as an indicator function. The agent is controlled to crawl if gc = 1 otherwise if gc = 0. Since crawling is a status that needs to be maintained, the crawling height is defined as the vertical position of its highest point at each time step.
h̃c(t) = max s∈S h(s, t) (3)

heading: 3.2 LOSS FUNCTIONS
text: Consider a locomotion task where an agent is instructed to move at a specified velocity, we think a proper loss function should satisfy the following requirements:
1. Periodicity. The agent is expected to move periodically following a predefined cyclic activation signal. Therefore, the loss function should encourage periodic motions.
2. Delayed evaluation. Due to inertia, it takes time for the agent to start running or adjust running velocity. Therefore, an ideal loss function should take this delay into consideration.
3. Fluctuation tolerance. During a single running cycle, requiring the center of mass to move at a constant velocity at each time step may induce highly fluctuated losses. To avoid that, our loss function should be smooth during one motion cycle.
Task loss. For each task, we defined a tailored loss shown in equation 4 according to the requirements above. The loss for running Lv is defined as the accumulation of the difference between the target and the agent’s velocity from start to current time step. For jumping, we define a sparse loss Lh, which only evaluate once for each period. The crawling loss Lc is defined as the accumulation of the crawling height if the loss is applicable.
L = λv ∑ n∈T P∑ t=Pr
(ṽ(t)− gv(n))2︸ ︷︷ ︸ Lv
+λh ∑ n∈T
(h̃h(n)− gh(n))2︸ ︷︷ ︸ Lh
+λc ∑ n∈T P∑ t=0
gc(n)h̃c(t+ nP)︸ ︷︷ ︸ Lc (4)
The λv, λh, λc are weights for losses of dfferent tasks. Our loss function accumulates the contributions of all steps in a “sliding window”. This prevents the loss function from being too small as well as the vanishing gradient problem. We find that having an accumulated loss evaluated at every step works much better when we want to control the speed of the agent explicitly.
Regularization. Since the agents may keep shaking when no target velocity or height are given, we introduce an actuation loss as a regularization term. The intuition is to penalize comparatively large actuations when small goals value is given:
La = ∑ n∈T ∑ t∈[0,P] (A(t)− µ||gv(n)||)2 (5)
where A(t) is the actuation output of NN, µ is a normalizer. The actuation loss balances the importance of input channels of the NN, i.e., different input channels produce similar contributions to the first layer of the trained NN.

heading: 3.3 NETWORK ARCHITECTURE
text: We use two fully connected (FC) layers with sine activation function as the neural network. The network input vector consists of three parts:
1. Periodic control signal of the same period with different phases, to encourage periodic actuation. In this work, we use sin waves as the periodic control signal. 2. State feature vector extracted from the current state of the agent. Using the positions and velocities of all the vertices of the simulated object as the input feature is not a good idea since these quantities are neither translation-invariant nor rotation-invariant. We use a “centralized pose” to remove global translation and rotation information. To be more specific, we subtract the position of the center of mass (CoM) in both 2D and 3D cases to remove the global translation. In 3D cases, we also compute the rotation around the vertical axis of the agent to remove its global rotation. 3. Targets that encode the task information such as running velocity, jumping height, and orientation for 3D agents. Different tasks use different channels to represents their target values. During the training, we assign random values to these targets. When validating the model, we always test the agent on a fixed set of target values. In interactive settings, we assign various target values to control the agent’s motion. To amplify the importance of the targets, we duplicate these channels multiple times, try to make them comparably important among the other input features.

heading: 3.4 TRAINING
text: In the training process, we run one simulation in one training iteration. The duration of the simulation is divided into several periods P as mentioned in the section 3.1. For each P , a goal G contains several tasks is assigned. For example, G = {gv, gh} represents that the agent is expected to perform both running and jumping in this period. The goal and its target values are generated randomly in a uniform distribution. They are kept fixed inside one period P , but varied between different periods.

heading: 4 EXPERIMENTS AND ANALYSIS
text: In this section, we systematically investigate the key factors that contribute to the training, evaluate our framework through ablation studies and comparisons against PPO Schulman et al. (2017), a state-of-the-art reinforcement learning algorithm. All experiments in the ablation studies are performed at least three times. The results of ablation studies are summarized in Table 1 and Fig. 5.

heading: 4.1 VALIDATION METRIC
text: In validation, we test all the trained agents on a fixed set of goals. These goals are represented by combination of targets targets values, which are uniformly sampled by interpolating between the lower and upper bound of all targets appeared in training.

heading: 4.2 DIFFERENTIABLE PHYSICS GRADIENT ANALYSIS
text: One concern in differentiable physics based learning is the that the gradient may explode during the back-propagation through long simulation steps. Here we visualize the distribution of gradients norm of a whole training process in Fig. 6. It can be observed that most values of the gradients are distributed inside the interval [-10, 10] in log scale, which indicates that our learning approach provides stable gradients during training.
Adam vs SGD. Hu et al. (2020) shows that stochastic gradient descent (SGD) with differentiable physics can achieve satisfying results for tasks with single goals, e.g. running along one direction. However, this no longer holds for tasks with multiple goals. To investigate the differences, we perform a learning task with multiple goals on a 2D mass spring system. In this experiment shown in the left of Fig. 5, SGD (curve Full-OP) makes minor progress while our modified Adam optimizer with reduced momentum drops the validation obviously within 10000 iterations. One key difference between Adam and SGD is that Adam utilises the momentum while vanilla SGD doesn’t. To further
investigate the the effectiveness of the momentum hyperparameters β, we perform a grid search for both β1 and β2 whose default values are 0.9 and 0.999 respectively. The results presented in table 2 show that reducing β2 to 0.9 delivers us the best performance. We therefore use this Adam optimizer with reduced momentum in all our experiments.

heading: 4.3 LEARNING TECHNIQUES
text: Batching We adopt batching into differentiable physics and perform a series of experiments as shown in the right of Fig. 5 to verify its effectiveness, with batch size 1, 8, 32, and 64. It can be observed that without batching, there is a large variance in training loss and it hardly converges. Batching with a proper size helps reduce the variance and effectively improve the training performance. We also find that an overlarge batching size does not help training obviously. For example, batching with size 64 only improve the performance marginally.
Loss design We compare our tailored loss function with a naive loss function. The naive loss function defines the velocity as the position difference of the center of mass between consecutive time steps, i.e., requiring the agent to move at a constant velocity at each time step. The experiments show that the agent under the naive loss design can not be properly trained to run, i.e., there are almost no progress in running loss dropping.
0 2000 4000 6000 8000 10000 Training Iteartions
0.75
0.80
0.85
0.90
0.95
1.00
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - Alpaca sin tanh relu gelu sigmoid
0 2000 4000 6000 8000 10000 Training Iteartions
0.75
0.80
0.85
0.90
0.95
1.00
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - Alpaca sin tanh relu gelu sigmoid
0 2000 4000 6000 8000 10000 Training Iteartions
0.75
0.80
0.85
0.90
0.95
1.00
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - Alpaca sin tanh relu gelu sigmoid
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
0 2000 4000 6000 8000 10000 Training Iteartions
0.70
0.75
0.80
0.85
0.90
0.95
1.00
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - Monster
sin
0 2000 4000 6000 8000 10000 Training Iteartions
0.70
0.75
0.80
0.85
0.90
0.95
1.00
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - Monster
sin
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0 2000 4000 6000 8000 10000 Training Iteartions
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - HugeStool sin
0 2000 4000 6000 8000 10000 Training Iteartions
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - HugeStool sin
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0 2000 4000 6000 8000 10000 Training Iteartions
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - Stool sin
0 2000 4000 6000 8000 10000 Training Iteartions
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ta sk
L os
s
Activation Function Ablation - Stool sin
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0 20 40 60 8000 10000 Training Iteartions
0.6
0.7
0.8
0.9
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
Activation Function Ablation - Alpaca sin tanh relu gelu sigmoid
0 2000 4000 6000 8000 1000 Training Iteartions
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
Activation Function blation - lpaca sin tanh relu gelu sig oid
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
0 20 40 60 8000 10000 Training Iteartions
0.5
0.6
0.7
.8
.
1.0
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
Activation Function Ablation - Monster
sin
0 2000 4000 6000 8000 10000 Training Iteartions
0.5
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
Activation Function Ablation - Monster
sin
0.0 0.2 0.4 . . 0.0
0.2
0.4
. . 0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0 20 40 60 8000 10000 Training Iteartions
0.4
0.5
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
Activation Function Ablation - HugeStool sin
0 2000 4000 6000 8000 10000 Training Iteartions
0.4
0.5
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
Activation Function Ablation - HugeStool sin
0.0 0.2 0.4 . . 0.0
0.2
0.4
0.6
0.8
1.0
. . 0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0 20 40 60 8000 10000 Training Iteartions
0.4
0.5
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
Activation Function Ablation - Stool sin
0 2000 4000 6000 8000 10000 Training Iteartions
0.4
0.5
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ru n.
L os
s
Activation Function Ablation - Stool sin
0.0 0.2 0.4 . . . 0.0
0.2
0.4
0.6
0.8
1.0
. . . . . . .
.
.
.
.
.
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0.80
0.85
0.90
0.95
No rm
liz ed
V al
id at
io n
Ju m
p. L
os s
No rm
liz ed
V al
id at
io n
Ju m
p. L
os s sin
0 2000 4000 6000 r i i It rti s
0.80
0.85
0.90
0.95
1.00
No rm
liz ed
V al
id at
io n
Ju m
p. L
os s
ti ti ti l ti - st r
si
0 2000 4000 6000 8000 10000 Training Iteartions
0.80
0.85
No rm
liz ed
V al
id at
io n
Ju m
p. L
os s
. . . . . . . . 0.6 0.8 1.0 .
.
.
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0 2000 4000 6000 r i i It rti s
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ju m
p. L
os s
ti ti ti l ti - t l si
0 2000 4000 6000 8000 10000 Training Iteartions
0.6
0.7
0.8
0.9
1.0
No rm
liz ed
V al
id at
io n
Ju m
p. L
os s
Activation Function Ablation - HugeStool sin
. . . . .
.
.
.
.
. . . . . 0.6 0.8 1.0 .
.
.
.
.
.
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
0 2000 4000 6000 r i i It rti s
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
No rm
liz ed
V al
id at
io n
Ju m
p. L
os s
ti ti ti l ti - t l si
0 2000 4000 6000 8000 10000 Training Iteartions
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
No rm
liz ed
V al
id at
io n
Ju m
p. L
os s
Activation Function Ablation - Stool sin
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
.
. . . . . . .
. . . 0.6 0.8 1.0 .
.
.
.
.
.
0.0 0.2 0.4 0.6 0.8 1.0 0.0
0.2
0.4
0.6
0.8
1.0
Activation function To validate the effectiveness of our sin activation function, we replaced it for hidden layer with various popular activation functions tanh, relu, gelu and sigmoid. For the output layer, the output actuation is expected to be limited in the range -1 to 1 due to the physical constraints of an actuator, e.g., a spring should not be overly compressed or stretched. Therefore we made some modifications for the activation functi s in the output layer. To be more specific, we clamp the values larger than 1 or smaller than -1. In addition, for activation functions whose values are constantly zero (or very small) when input is negative, we modify the relu and gelu by subtracting the value by 1, for sigmoid we map the its value from [0, 1] to [-1, 1]. The results in Fig.7 show that using sin as activation functions gives the best performance among all three losses. For more results of activation functions ablation studies, please check the appendix
7

heading: 4.4 IMPORTANCE OF INPUT FEATURES
text: Ablation on input features The input features consist of three parts: periodic control signal, state vector, and targets, as mentioned in section 3.3. We investigate each part of the contribution to training by ablation studies as shown in Fig.8. The results show that each of the three parts serves as a necessary component for a converged training. The periodic control signal serves as an important role. The agents can hardly move without this signal as the corresponding task loss does not drop along the training (shown by the red curves in Fig.8). The targets are another essential part, without which the optimizer can hardly work after a minor progress as shown by the yellow curves in Fig.8. Although the state vector is less important in terms of reducing the loss, the overall task performance is degraded without the state vector as shown by the green curves in Fig.8.

heading: 4.5 BENCHMARK AGAINST REINFORCEMENT LEARNING
text: We implement the standard proximal policy optimization (PPO) benchmark in multiple goals settings as shown in Fig. 9. The PPO agent can make progress in single goal learning and learn certain locomotion patterns to move toward the goal. Yet as the task gets complicated, PPO often gets trapped in the local minima. We try our best to fine-tune the PPO hyper-parameters, but still find it tends to overfit to one goal but fails to find a balance between different goals. It can be observed in Fig. 9 where the agents trained by PPO can learn to jump but struggle to run.

heading: 5 RELATED WORK
text: Differentiable simulation Differentiable physical simulation is getting increasingly more attention in the learning community. Two families of methods exist: the first family uses neural networks to approximate physical simulators (Battaglia et al., 2016; Chang et al., 2016; Mrowca et al., 2018; Li et al., 2018). Differentiating the approximating NNs then yields gradients of the (approximate) physical simulation. The second family is more accurate and direct: many of these methods using differentiable programming (specifically, reverse-mode automatic differentiation) to implement physical simulators (Degrave et al., 2016; de Avila Belbute-Peres et al., 2018; Schenck & Fox, 2018;
Heiden et al., 2019; Hu et al., 2019; 2020; Huang et al., 2021). Automatic differentiation works well for explicit time integrators, but when it comes to implicit time integration, people often adopt the adjoint methods (Bern et al., 2019; Geilinger et al., 2020), LCP (de Avila Belbute-Peres et al., 2018) and QR decompositions (Liang et al., 2019; Qiao et al., 2020). In this work, we leverage DiffTaichi (Hu et al., 2020), an automatic differentiation system, to create high-performance parallel differentiable simulators and the built-in NN controllers.
Locomotion skill learning Producing physically plausible characters with various locomotion skills is a challenging problem. One classic approach is to manually design locomotion controllers subject to physics laws of motion in order to generate walking Ye & Liu (2010); Coros et al. (2010) or bicycling Tan et al. (2014) characters. These physically-based controllers often rely on a complex set of control parameters and are difficult to design and time-consuming to optimize. Another line of research suggests using data to control characters kinematically. Given a set of data clips, controllers can be made to select the best fit clip to properly react to certain situations Safonova & Hodgins (2007); Lee et al. (2010); Liu et al. (2010). These kinematic models use the motion clips to build a state machine and add transitions between similar frames in adjacent states. Although being able to produce higher quality motions than most simulation-based methods, the kinematic methods lack the ability to synthesize behaviors for unseen situations. With the recent advances in deep learning techniques, attempts have been explored to incorporate reinforcement learning (RL) into locomotion skill learning Peng et al. (2015); Liu & Hodgins (2017); Peng et al. (2018); Park et al. (2019). These modern controllers gain their effectiveness either from tracking high-quality reference motion clips or from cleverly designed rewards to imitate the reference. However, the exploration space for RL is usually prohibitively large to achieve complicated target motions. Carefully designed early termination strategies Ma et al. (2021); Won et al. (2020), better optimization methods Yang & Yin (2021), and adversarial RL schemes Peng et al. (2021) enable these RL-based methods to achieve richer behaviors. However, it is still time-consuming to get a well-trained RL model in complex tasks such as soft robot control. Our method, on the other hand, leverages the differentiable simulation framework and can achieve much better convergence behavior compared to RL-based methods.

heading: 6 LIMITATIONS AND CONCLUSIONS
text: Currently, difficulty of the learning tasks depends on robot design and physical parameters. The training performance may be degraded with improper robot designs. For instance, in our stool case where unnecessary actuators on its body are allowed, it is more difficult to achieve good results. Offloading the physical parameter tuning and robot designing to an automatic pipeline will be an interesting future research direction.
To summarize, we have presented an effective end-to-end differentiable physics based learning framework for soft robot complex locomotion skills learning. We systematically enhance classical differentiable physics learning systems with a suite of techniques and investigated the key factors contributing to the training in detail. We show that our framework can provide stable gradients without explosion during a simulation with hundreds of steps. Together with batching, the modified Adam dramatically outperforms the stochastic gradient descent (SGD) on complex tasks with multiple goals. Benefited by our tailored loss functions, the network take advantages of the three parts of the input features. We found that the periodic control signal dominates the actuation signal, while state vector and task goals have weaker effects. The over-dominated periodic signal may induce high-frequency noises. To balance the importance between the input features, we additionally pose an actuation loss as a regularization term. It can effectively suppress the noises, which help the agents move more smoothly. In addition, we compare our method with the state-of-the-art reinforcement learning proximal policy optimization. Our method shows advantages in both training performance and convergence efficiency on tasks with multiple goals.
Our framework enables users to flexibly design robots and to teach them locomotion skills. The trained agents can be manipulated to smoothly switch locomotion tasks such as running, jumping, and crawling with different speeds and orientations, interactively. To the best of our knowledge, this is the first time a learning system based on differentiable physics can deliver controllers with such robustness, flexibility, and efficiency.

heading: A REINFORCEMENT LEARNING SETTING
text: Environment. We use the open-source implementation of PPO (Kostrikov (2018)) in our environments. Part of important hyper-parameters are listed in the table 3.
Reward Design. We design a reward function for PPO based on our loss functions. Recall the velocity part of equation 4. We define the total velocity reward as:
Rv = λv ∑ n∈T ∑ t∈[Pr,P] gv(n) 2 − (ṽ(t)− gv(n))2 (6)
If we directly split the reward into each time steps by t, the rewards of first Pr time steps would be zero and it is too difficult for PPO to learn the policy. Therefore, we modified the reward function to tackle this issue. By equation 1, we have
Rv = λv ∑ n∈T ∑ t∈[Pr,P] gv(n) 2 − [ 1 Pr (c(t)− c(t− Pr))− gv(n) ]2 (7)
= λv P2r ∑ n∈T ∑ t∈[Pr,P] [Prgv(n)]2 − [c(t)− (c(t− Pr) + Prgv(n))]2 . (8)
We can define a function fn(t′, t) = [Prgv(n)]2 − [c(t′)− (c(t− Pr) + Prgv(n))]2 and substitute it back to equation 8. The reward can be re-written as
Rv = λv P2r ∑ n∈T ∑ t∈[Pr,P] fn(t, t) (9)
Refer to that for any t ∈ P , we have fn(t− Pr, t) = 0, which means
Rv = λv P2r ∑ n∈T ∑ t∈[Pr,P] ∑ ∆t∈[0,Pr) fn(t−∆t, t)− fn(t−∆t− 1, t) (10)
Let t′ = t−∆t
Rv = λv P2r ∑ n∈T ∑ t′∈[0,P] min(t′,Pr)∑ ∆t=0 fn(t ′, t′ + ∆t)− fn(t′ − 1, t′ + ∆t) (11)
= λv P2r ∑ n∈T ∑ t′∈[0,P] min(t′,Pr)∑ ∆t=0 − [c(t′)− (c(t′ + ∆t− Pr) + Prgv(n))] 2 +
[c(t′ − 1)− (c(t′ + ∆t− Pr) + Prgv(n))] 2 (12)
So we can split the whole reward into each time step t as:
Rv(t) = min(t,Pr)∑
∆t=0
[c(t− 1)− (c(t+ ∆t− Pr) + Prgv(n))]2−
[c(t)− (c(t+ ∆t− Pr) + Prgv(n))]2 (13)
Recall that the goal of a jumping task is to reach a specified height goal gh, the intuition for jumping reward is to encourage the agent to improve its maximum height until reaching the height.
Rh(t) = λh ∑ n∈T (min s∈S h(s, t+ nP)− gh(n))2 − (h̃h(n)− gh(n))2 (14)
where the mins∈S h(s, t+nP) is the agent height at current time step as mentioned in 3.1 and h̃h(n) is the max height record during one task period.

heading: B NETWORK ARCHITECTURE
text: Hidden Layer (64)
Actuation
FC1 with sine
FC2 with sine
Input Features
Our network architecture is shown in Fig. 10. It is a two layers fully connected network with a 64 channels hidden layer. The number of input and output channels varies according to different agent designs.

heading: C SIMULATION SETTING
text: Mass-spring systems We adopt the classic Hookean spring model to represent the elastic force and use dashpot damping Baraff & Witkin (1998) as the damping force to simulate mass-spring system. We found that the drag damping model used by Hu et al. Hu et al. (2020) damps the gradient of the NN controller while the dashpot damping model is able to generate vividly changing gradients. Adopting dashpot damping model makes our agent more flexible.
Material point methods We use MLS-MPM Hu et al. (2018) as our material point method simulator. We further applied the affine particle-in-cell method Jiang et al. (2015) to reduce the artificial damping. Due to the nature of MPM as a hybrid Lagrangian-Eulerian method, it always requires a background grid during the simulation. Naive MPM implementations fix the position of the background grid, hence limit the moving range of the agent. We apply a dynamic background grid that follows the agent all the time to overcome this problem.
D INTERACTIVE CONTROL
Smooth interpolation between different tasks is a key advantage of our approach. In interactive settings, users can smoothly control the agent’s motion via input devices such as a keyboard as shown in Fig. 11. Please refer to our supplemental video for more details.
E 3D RESULTS
In addition to 2D cases, our framework can be applied to 3D cases seamlessly. In a 3D space, an agent is additionally expected to control its orientation on a plane, i.e., rotate. We design several 3D agents (shown in Fig. 2) and train them to run and rotate. These trained agents can achieve specified goals given control signals after trained for several hundred to a few thousand iterations. Please check the supplemental videos for more details.
Loss function in 3D The agents in 3D are trained with the loss function below. The loss for running Lv is inherited from 2D settings. The rotation loss Lr is defined as the accumulation of centralized point-wise distance between current position and target rotated position. L = λv ∑ n∈T P∑ t=Pr
(ṽ(t)− gv(n))2︸ ︷︷ ︸ Lr
+λr ∑ n∈T P∑ t=Pr ∑ s∈S
((s(t)− c(t))−R(s(t− Pr)− c(t− Pr)))2︸ ︷︷ ︸ Lr
(15)
R =
[ cos(Prgω) 0 − sin(Prgω)
0 1 0 sin(Prgω) 0 cos(Prgω)
] , (16)
where gω is the target angular velocity assembled in input features. The format of R is derived by rotation matrix along XZ-plane.
Note, for solving running tasks, the rotation loss also plays an important role. A significant issue of 3D agents running is accumulated orientation error. Utilizing identity rotation matrix R helps adjust agent’s posture. In turn’s of rotation tasks, the rotation center should be maintained still which is equivalent to zero velocity. These two losses complement with each other to achieve the final goals.
Friction model. In real world, agents running is driven by static friction between ”feet” and ground. We configure and train agents with different contact models in our experiments. With zero friction applied, training loss shows no evidence decreasing and the agent is not able to move no matter what actions patterns it learns. With a fully sticky surface, where only upwards vertical velocity is kept and all other components of velocities are projected to zero after contact, the 2D agents work well. However, in 3D, it turns out that the agent tends to stuck. In practice, we adopt the classic slip with friction model from physically-based simulation which involves both kinetic friction and static friction. These two categories of contact are determined by how much pressure force is applied on the surface and can be treated like a branched function numerically which can be handled by our auto-differential system. Please see the supplemental video for a visual comparison between slip and sticky friction models.

heading: F GENERALIZATION TO MANIPULATION TASKS
text: Here we provide results on simple manipulation tasks in additional to locomotion tasks. In these tasks, the agent is expected to manipulate a object (marked in purple) to hit a target point (marked in green). This task has similar periodicity property like locomotion. We designed two scenarios: Juggle and Dribbble and Shot shown in Fig.13. The visual results are also shown in the video in supplemental materials.
Juggle In this scenario, the target point appears in the sky above the the agent. The agent is expected to ’juggle’ the object to the target point.
Dribble and Shot In this scenario, the target point keeps moving toward left. The agent is expected to carry the object and shot it to the target point.
The experiments indicate that our method has the possibility to generalize to tasks beyond locomotion.

heading: G DETAILED RESULTS ON FURTHER EXPERIMENTS
text: Here we show more detailed results on further experiments, including ablation studies on different agents, validation loss on different tasks for agents and more network weights visualization results.

Review summaries:
This work proposes to learn locomotion skills on soft robots made of springs. It proposes to leverage differentiable physics along with a NN controller based on SIREN, which allows to directly learn policies by minimizing loss functions defined on trajectories. The methods allows to learn locomotion and jump behaviors on a various set of 2D and 3D mass-spring systems.

This paper provides a practical guide on how to train an NN-based controller for locomotion tasks with robustness and efficiency.   The tasks are categorized into motion primitives including running, jumping, and crawling, each with a manually designed objective function. Considering the periodic nature of locomotion, the author proposed to use sine activation layers and periodic signals in different phases as input. The experiments show that a larger batch size can lead to faster convergence. Compared to SGD, this paper points out that Adam is better for the tasks. The supplementary video presents that the trained networks enable real-time control of complex skills. 

This paper proposes to use differentiable simulators for learning locomotion skills for deformable characters. The key contribution is the experiments and analysis of a large number of design choices, including activation function, reward design, optimizer, batch size, etc. The paper finds that by combining differentiable simulator with the right design choices, the character can learn general and complex skills with one policy. The paper also shows that this combination can significantly outperform PPO.

The paper proposes to learn neural network controllers in a differentiable simulator for locomotion tasks. A set of technical improvements are made and make the controllers more robust and require fewer samplers to learn. These include: using periodic activation functions, design of loss functions, and using large batch size. The paper shows learned locomotion skills of running, jumping, and crawling for different soft robots, which can be interactively controlled.