On the Numerical Reliability of Differentiable Physics-Based Optimization for Robotic Material Manipulation
Authors: Xintong Yang, Minglun Wei, Yu-Kun Lai, Ze Ji
Organizations: School of Engineering, Cardiff University, Cardiff, United Kingdom · School of Computer Science and Informatics, Cardiff University, Cardiff, United Kingdom · College of Mechanical and Electrical Engineering, Hohai University, China
Differentiable physics is increasingly used in robotic material manipulation for system identification, trajectory or skill optimization, demonstration generation, and robot or end-effector design. These applications depend on gradients propagated through long, contact-rich simulation rollouts. We study the numerical reliability of those gradients using two Material Point Method (MPM) system-identification benchmarks derived from elastoplastic and granular manipulation. The benchmarks provide controlled cases for three effects that also arise in broader differentiable physics-based optimization. GPU many-to-one sums whose order depends on thread scheduling changed long-horizon gradients and reversed the sign of one parameter gradient relative to a deterministic reference. Finite-difference checks became less reliable for longer rollouts because repeated-run loss variation grew much faster than the loss change produced by the tested parameter perturbations. Observation and loss definitions changed optimization behaviour and the solution preferred by an independent metric. These results motivate reproducible accumulation, finite-difference validation that compares perturbation-induced loss changes with repeated-run variation, and explicit reporting of objective construction when differentiable simulation is used for robotic optimization.
Figures & tables
Fig. 1: Robotic material-manipulation settings motivating this study. Top: Real elastoplastic interaction and the corresponding MPM simulation from the DPSI-derived setting [ 8 ] . Bottom: Real and simulated granular digging from the DDBot-derived setting [ 9 ] .
Appendix figures & tables4 assets
Supplementary material from the paper’s appendix.
Appendix
Setting
Configuration
DPSI-derived elastoplastic
Fixed-corotated elasticity with von Mises plasticity; 218, 451, or 876 particles; 94 global steps with 50 MPM substeps per step; PRT-CD, PRT-EMD, PCD-CD, and PCD-EMD.
DDBot-derived granular
Hencky/Saint Venant-Kirchhoff elasticity with Drucker-Prager plasticity; 27,440 particles; 200 or 311 global steps with 20 MPM substeps per step; HMD and PCD-EMD.
Appendix
TABLE I: Benchmark configurations used in the objective matrix.
Global steps
f64 worst error
f32 loss spread
FD time
2
4.4×10−6 to 1.5×10−5
2.2×10−7 to 1.4×10−6
25 s
5
2.75×10−5
2.9×10−5
38 s
10
2.93×10−5
7.1×10−5
59 s
20
3.41×10−1
4.7×10−5
99 s
40
2.70×10−1
2.8×10−5
178 s
Appendix
TABLE II: DPSI-derived horizon experiment. The f64 column reports the worst relative AD/FD discrepancy across the tested material components.
Objective
Descending epochs
Mean HMD (mm)
PRT-CD
78/80 (97.5%)
3738.4
PRT-EMD
54/80 (67.5%)
3698.6
PCD-CD
43/80 (53.8%)
3711.2
PCD-EMD
46/80 (57.5%)
3708.6
Appendix
TABLE III: Reference-density DPSI-derived optimization behaviour. Descent count is the number of epochs whose post-update objective is lower than the entry value for that epoch, aggregated over two seeds and both precisions. HMD is an independent terminal evaluation.
Case
Runtime or outcome
DPSI 94-step SDF contact, CPU loop vs. GPU kernels
260.821 s → 20.488 s (12.73 × )
DPSI FD validation, 2 steps vs. 40 steps
25 s → 178 s
DDBot HMD epoch, 200 steps
about 115 s; about 154 s with the extended line-search set
DDBot PCD-EMD epoch, 311 steps
about 170 s; about 227 s with the extended line-search set
DPSI particle-density change
about +1.5% wall time for denser f32 cases; recorded device memory unchanged
Deterministic repeat cells
5/5 bit-identical across selected f32/f64 and DPSI/DDBot-derived cases
Three properties determine whether a differentiable simulator can drive gradient-based optimization through contact: simulation accuracy, gradient reliability, and per-iteration cost. Tape-based engines such as MJX and Newton Semi-Implicit require timesteps small enough to keep contacts numerically tractable, and their backpropagation memory grows linearly with the number of timesteps T. Surrogate models bound memory by approximating contact away, but the resulting gradients lose the geometry the optimization depends on. We present Ostrich, a GPU-accelerated rigid-body simulator that resolves hard contacts and friction with non-smooth Newton iteration at large timesteps (h ~ 0.1 s), and differentiates the converged residual via the implicit function theorem, reusing the forward Schur complement to compute the adjoint at O(1) memory per timestep. On real-robot trajectories over a pallet obstacle, Ostrich holds MuJoCo's sim-to-real accuracy up to a 50x larger timestep. Its gradients converge from random initializations where MJX descends slowly and Newton Semi-Implicit stalls; a warm iteration runs 211x faster than MJX's and 4.7x faster than Semi-Implicit's. On the same scene Ostrich differentiates 8,192 parallel worlds on a single 24 GB GPU, sustaining 29x checkpointed MJX's optimization throughput; without checkpointing both baselines exhaust memory at far fewer worlds. We close with a gradient-based trajectory optimization demonstration over triangle-mesh terrain across a 10 s horizon, a setting where prior engines either restrict to primitive geometry or face the convergence and memory limits shown above.
Aleš Kučera, Karel Zimmermann
Department of Cybernetics, Faculty of Electrical Engineering, Czech Technical University in Prague, Prague, Czech Republic
Differentiable simulation of soft bodies is a foundation for system identification, trajectory optimization, and Real2Sim transfer. Yet, existing methods such as the differentiable Projective Dynamics (DiffPD) struggle when faced with heterogeneous materials with extreme stiffness contrasts, hyperelasticity under large deformations, and contact-rich interactions, which are common scenarios in the real world. We present DiffPhD, a unified GPU-accelerated differentiable Projective Dynamics framework for heterogeneous materials that tackles these intertwined challenges simultaneously. Our key insight is a careful integration of: (i) stiffness-aware projective weights to embed heterogeneity into the global system; (ii) trust-region eigenvalue filtering lifted to the backward pass for stable hyperelastic gradients and a type-II Anderson Acceleration scheme with dual-gate convergence to stabilize forward iteration under large stiffness contrasts; and (iii) a unified GPU pipeline that reuses a single sparse factor across forward, backward, and contact computations, with stiffness-amplified Rayleigh damping folded into the same factor for heterogeneity-aware dissipation at zero recurring cost. DiffPhD achieves strict gradient accuracy while delivering up to an order-of-magnitude speedup over prior differentiable solvers on heterogeneous, hyperelastic, contact-rich benchmarks. Crucially, this speedup does not come at the cost of stability: DiffPhD remains convergent on stiffness contrasts up to 100x where prior PD solvers degrade. This unlocks end-to-end gradient-based optimization on regimes previously bottlenecked by either solver fragility or per-iteration cost -- shell--joint composite creatures, soft characters wielding stiff weapons, and soft-gripper robotic manipulation -- all handled within a single forward--backward pass.
Shih-Yu Lai, Sung-Han Tien, Jui-I Huang +9
National Taiwan University · MoonShine Animation Studio · National University of Singapore +3
Robotics simulators have improved significantly in computational speed and scalability, enabling them to generate years of simulated data for complex systems in minutes or hours. Despite these advances, efficiently and accurately computing simulation derivatives remains an open challenge. Addressing this would accelerate the convergence of reinforcement learning and trajectory optimization algorithms, particularly for contact-rich problems. This paper introduces a unifying framework for robotic simulation that accounts for all factors, including dynamics, collisions, and friction. The resulting algorithm computes analytical derivatives of the simulation by implicit differentiation, explicitly handling the intrinsic non-smoothness of the collision and frictional stages while exploiting the sparsity induced by the multi-body structure. Benchmark results demonstrate state-of-the-art performance, with timings ranging from 5μs for a 7-dof manipulator to 95μs for a 36-dof humanoid, an improvement of at least two orders of magnitude over alternative methods. Implemented in C++, the code will be open-sourced after the review process to support applications such as simulation-driven learning and real-time control.
Quentin Le Lidec, Louis Montaut, Yann de Mont-Marin +3
Inria, Ecole normale sup´erieure CNRS, PSL Research University 75005 Paris, France