Differentiable Simulation of Hard Contacts with Soft Gradients for Learning and Control
Authors: Anselm Paulus
Andreas René Geist
Pierre Schumacher
Vít Musil
Simon Rappenecker
Georg Martius
Venue: ICLR 2026 Poster (2026)
Keywords: differentiable simulation
contact forces
musculoskeletal control
Primary area: applications to robotics, autonomy, planning
TL;DR: Improve gradients of penalty-based simulators to enable efficient gradient-based optimization with contact-rich robotic systems.

Abstract: Contact forces introduce discontinuities into robot dynamics that severely limit the use of simulators for gradient-based optimization. Penalty-based simulators such as MuJoCo, soften contact resolution to enable gradient computation. However, realistically simulating hard contacts requires stiff solver settings, which leads to incorrect simulator gradients when using automatic differentiation. Contrarily, using non-stiff settings strongly increases the sim-to-real gap. We analyze penalty-based simulators to pinpoint why gradients degrade under hard contacts. Building on these insights, we propose DiffMJX, which couples adaptive time integration with penalty-based simulation to substantially improve gradient accuracy. A second challenge is that contact gradients vanish when bodies separate. To address this, we introduce contacts from distance (CFD) which combines penalty-based simulation with straight-through estimation. By applying CFD exclusively in the backward pass, we obtain informative pre-contact gradients while retaining physical realism.

Published as a conference paper at ICLR 2026
DIFFERENTIABLE SIMULATION OF HARD CONTACTS WITH SOFT GRADIENTS FOR LEARNING AND CONTROL
Anselm Paulus ; A. René Geist"
Pierre Schumacher", Vít Musil Simon Rappenecker',
Georg Martius
1 University of Tübingen, Tübingen, Germany
2 MyoLab.AI. New York City, USA
3 Masaryk University, Brno, Czechia
(anselm, watssmm.ppuulsss.. paulus, rene. geist, georg martius] @un bunriuubenentttteennnnneee de
ABSTRACT
Contact forces introduce discontinuities into robot dynamics that severely limit the use simulators for gradient- optimization. Penalty-based simulators such MuJoCo, soften contact resolution to enable gradient computation. However, realistically simulating hard contacts requires stiff solver settings, which leads incorrect simulator gradients when using automatic differentiation. Contrarily, using non-stiff settings strongly increases the sim-to-real gap. We analyze penalty- based simulators pinpoint why gradients degrade under hard contacts. Building on these insights, we propose DifiMJX, which couples adaptive time integration with penalty-based simulation substantially improve gradient accuracy. A second challenge is that contact gradients vanish when bodies separate. To address this, we introduce contacts from distance (CFD) which combines penalty- -based simulation with ttraighttthrough estimation. By applying CFD) exclusivel in the backward pass, we obtain informative pre-contact gradients while retaining physical realism. Project page: https: / ggtttb....../gggggg com mart ius- Lab/ / dii fm
INTRODUCTION
Gradients have powered major advances in machine learning ranging from video language models to robot control. In robotics, imitation and reinforcement learning widely rely on gradient- based optimization. Yet, despite the dominance of sim-to-real techniques that leverage robot simulators for policy learning (Tan et al., 2018; Lee et al., 2020; Andrychowicz et al., 2020; Radosavovic et al., 2024; Li et al., 2022), most methods conspicuously avoid using simulator gradients. This is missed opportunity considering that simulator gradients offer direct route to updating actions and learning model parameters. If these gradients were accurate, could fit simulators to real data significantly narrowing the sim-to-real gap and accelerate policy learning - enabling policies for new tasks to be trained in seconds rather than hours. Given the utility of simulator gradients, what prevents their use in robot learning? ? In practice, two fundamental issues hinder the use of simulator gradients: (i) discontinuities arising from contacts yield erroneous gradients, and (ii), if objects do not touch, contact gradients are zero. In this work, we chart path that tackles both challenges.
Differentiating through contacts. The choice of contact model has significant implications simulator's ' differentiability. Complementary-basee solvers, such as Taylor et al. (2022), compute contact forces exactly, such that sudden jumps in the dynamics aggravate gradient computation. Hence, recent literature on differentiating through complementary-bast bassed solvers proposes analytical reformulations of the dynamics, e.g. based on the implicit function theorem (Werling et al., 2021; Taylor et al., 2022), resorts to randomized smoothing (Tassa & Todorov, 2010; Duchi et al., 2012; Suh al., 2022; Bouyarmane al., 2009; Xu al., 2010). While differing in computational cost, both approaches have similar empirical performance Pang et al. (2023); Schwarke et al. (2024). To improve computational efficiency, MuJoCo (Todorov, 2014) reformulates complementary constraints
'These authors contributed equally.

---

Published as a conference paper at ICLR 2026
convex optimization problem, where a constraint's ability to generate force grows proportionally to the constraint violation. While such a penalty-based simulator can be made smooth, gradient inaceuracies increase with the contact stiffness, the relative contact velocity, and the integration step size. While common advice points to reducing the step size to improve gradient accuracy, using sufficiently small step sizes results in prohibitively slow simulation.
Computing gradients of stiff' simulations. In the sequel, we show that gradient errors in penalty- based simulators arise from computing a discrete approximation of continuous physics put simply, from numerically integrating stiff differential equations. Consequently, the gradients obtained via automatic differentiation are "incorrect" insofar as they do not align with those of the underlying continuous system, Although adaptive integration well studied, it's rarely used robotics and to the best of our knowledge its utility for differentiable simulation has not been explored. This perspective complements prior proposals for time-of impactt correction (Hu ct al., 2020; Schwarke et al., 2024) in impulse-baser simulators. We incorporate adaptive timestep integration into MuJoCo XLA, incurring a slight computational overhead while obtaining correct gradients in the presence of hard contacts, and remaining compatible with existing MuloCo libraries.
Computing gradients between non-colliding objects. Another obstacle for gradient-based optimization for policy generation and system identification is the non- informativeness of gradients about unrealized contacts, For example, when a robot's hand is not in contact with an object, then there no gradient directing it make contact for task facilitation. Therefore, drawing inspiration from prior works on contact-invariant optimization (Mordatch et al., 201200222 2012b;a), we propose Contacts From Distance (CFD) to address this challenge. However, naively intro- ducing artificial contact forces considerably alters the simulation, resulting in too large sim-to-real gap. In order to preserve the simulation realism, we propose using the introduce CFD solely in the gradient computation. We implement our changes in Mujoco XLA, and ad- ditionally fix some low-level collision routines to be truly differentiable. CFD complements position- based costs often used in robot learning by removing the need specify exact contact locations.
LoV
Fx
Loss
Simulator
TE-
L
P
collision detection & contact solver
Target
Figure I: Common computational graph for robot control synthesis.
to
In summary, DiffMJX and CFD provide additional tuning knobs to enable practitioners to set gradient correctness and inform optimization on how non-colliding objects would need to move to influence each other's s state. In turn, this work provides a new perspective on how to obtain useful gradients in penalty- based differentiable simulators thereby enabling parameter estimation and policy synthesis for collision-rich and high-dimensional systems.
2 ROBOT SIMULATION
As illustrated in Fig. we want use automatic differentiation obtain the correct gradients of a loss functional L( Fx, P) p) where the next state of the robotic system is governed by the discrete-ting dynamics 1++1 step( at. p) with the state IyK == xxxx) [gk: vyl time & consisting of the system's 's generalized position 9& E R" and velocity E R" control actions as E R"=,, and model parameters p E R.r. Typically, multi-body dynamics simulators consist of forward dynumics model and numerical integration method. The forward dynamics govern the system's 's acceleration via the equations of motion
MM (
(1)
with the joint-space inertia matrix (() RR. applied forces T (ar, a) E R. bias force c(x) E R"* constraint space Jacobian J(q) E Xn., and constraint forces f(x) fe + fr + fc E R consisting of the equality constraint, the generalized friction, and the contact constraint forces. The contact-free dynamics are typically derived via recursive multi-body algorithms (Featherstone, 2014) while computing contacts are resolved through an intricate interplay of collision detection and contact force optimization. As we will see throughout this work, numerical integration plays pivotal role in understanding how contact forces may hinder correct gradient computation. For that, we will resort to MuJoCo XLA concrete example of penalty-based simulator.
2

---

Published as a conference paper at ICLR 2026
2.1 CONSTRAINT RESOLUTION IN MUJoCo XLA
MuJoCo XLA (MJX) is a reimplementation of MuJoCo using the Python library JAX (Bradbury al., 2018), which enables GPU-paralleiizabl gradient computation via automatic differentiation. MuloCo has become the de facto standard in robotics, alongside other widely used simulators such as Bullet (Coumans & Bai, 2016), Drake (Russ Tedrake, 2019), Diff'aichi (Hu et al., 2020) and Nvidia PhysX (Liang et al., 2018) 8). MuloCo 's importance is underpinned by NVIDIA and Google recently releasing MuJoCo Warp (Howell, 2025). what follows, we provide a brief overview of how MuJoCo resolves contacts, which serves as a hands-on example of a penalty-based simulator.
Collision detection. Given the state x and geometry parameterizations of two bodies, a collision detector returns the signed distance r between potential contact points alongside with their surface normals. A contact is only considered active - that is, a contact point can exert contact force - < O. While collision detection is not the primary focus of this work, we later show how to remove discontinuities from its functions to ensure differentiability.
Contact force solver. A detailed description of MuJoCo is given in its documentation, here we provide a summary of key equations. Given Eq. (I), MuJoCo resolves constraints via a relaxation of Gauss' principle (GauB, 1829) which in its primal formulation reads
h(r) kr
d(
(b, = arg minl|z M (T -c) <) 1 + l1y aref Il (#.,)
d(rr)
w
(2)
dw
subject to Jgæ - ye 0, JFI - - 0, Je x yc E EK
do
with the friction cone (dual) KC" the regularizer R > O, and Huber norm || - || The reference acceleration Oret denotes the solver's target for the constraint space acceleration w. Drawing inspiration from Baumgarte (1972), aret follows damped harmonic oscillator
contact 0 contact active inactive
Figure 2: The position-level reference acceleration h(()) and impedance ddr) deter- mine the contact force magni- tudes that the solver can apply.
2 Oref,5 = (Jv); ker, = dw Le (Jv), - -2 Ts (3)
d(r.)
whose dynamics are determined by the impedance d(r), and the solref parameters consisting of the time constant Le and the damping
ratio da. The impedance d(r)) is the central function for determining constraint forces (such contact forces). As illustrated in Fig. 2, the impedance E [0, 1] is a function of the constraint violation r and is specified by the solimp parameters (do, the, w, midpoint, power) which define its shape nomial spline. MuloCo's documentation refers to the impedance as a "constraints ability to generate force" it determines aref and weights the cost for applying constraint forces in Eq. (2) via the diago- nal matrix R, which is computed as R2g+ Áa((1 d d,, where À denotes an approximation A diag (JM- In turn, the constraint is hard if R - O, and approaches an infinitely soft (i.c. non- existent) constraint in the limit R + co. If 1C" solely contains pyramidal or elliptic cone constraints, then problem (2) reduces to convex problem that is solved efficiently via an exact Newton method.
poly-
3 CORRECTING THE CONTACT GRADIENTS OF PENALTY-BASED SIMULATION
To evaluate the correctness of gradients in the presence of contacts, as illustrated in the top row of Fig. 3, we unroll the trajectory of several primitives bouncing against plane observe the final position and its gradient with respect to the initial velocity. In Fig. 3 (MJX row), at the default integration step of 0.002 s. we observe that the loss is oscillating at high frequency, which results in large fluctuations of the gradient. Similar oscillations have been discussed in previous works (Hu al., 2020; Schwarke et al., 2024). In particular. Hu al. (2020) emphasizes that the "time-of-impate (TOI) causes gradient oacillations when using ideal elastic or complementary collisions. Yet, TOI gradient errors differ in nature from errors observed in penalty- -basd simulators; as illustrated below.
Example Point collision: As illustrated in Fig. (left), a point mass starts at height go with velocity vo -1 and is integrated for N steps with semi-implicit Euler in the absence of gravity. It collides with a flat surface and bounces back up. Contact collision is either resolved via a minimal version of a penalty-based contact model found in MuJoCo or ideal elastic collision
3

---

Published as a conference paper at ICLR 2026
box
capsule
Pensity-based simulation
Ideal-elastic simulation
1.0
Re00
Large stepstze
Central differencess
MJK
Dring Gradient
90
0:55
ON
Grodlont Gridaart
0u00
643
Loan Too starte
-0.50
Cricken scaled 100
-50
1:00
INQ H-20003
and 00cr
Smalll sdepaalie
-0.2
surface
-0.4
Silucly chengs
-0.6
CorrEcIaed 15
scaied by
2 velccity %
postor
veiscity 42
Figure 4: Toy simulation of point mass colliding with a surface that either resorts to an ideal-elastic contact model used in DiffTaichi
Figure 3: Simulation of primitives thrown onto a surface. For stiff con- tacts, MJX's gradients of the toss dis- tance deviate from central differences, while DifIMIX maintains agreement.
penalty- -ased contact model similar to MuJoCo. Decreasing the integration stepsize h solely reduces errors in the penalty-based simulation, in the ideal- elastic simulation the gradient sign remains wrong.
as in (Hu et al., 2020). The corresponding JAX code is shown in Fig. 16 in the Appendix. The loss -= IgN Qr| is the distance between the point's final state gN and the target height @r 1. For ideal-elastic collisions, we observe sawtooth-like loss oscillations in L, resulting in the gradient V. a0 L with the wrong sign, independently of the stepsize as observed by Hu et al. (2020). For the penalty-based simulation, gradient oscillations notably reduce when lowering the stepsize.
We observe that state oscillations and the corresponding gradient artifacts are not a characteristic of the specific simulator. Instead, more fundamentally, they arise from time-discretization errors the ODE integration. While this affects both genalty-based and ideal-elastic simulators, we will see that the issue has to be addressed differently.
TOI correction does not fix gradients for penalty-based simulators, but small stepsizes do. For an ideal elastic collision, the ODE of our minimal example is piecewise linear: The dynamics are linear (due to the absence of gravity) before and after the contact, at which the velocity is inverted (see also Fig. 16 in the Appendix). The TOI approach proposed by Hu et al. (2020) exploits this structure by dynamically splitting the ODE into two linear segments at the time of contact. Integrating those separately thereby eliminates the discretization error and yields correct gradients. In penalty- based simulation, the ODE is linear before and after the collision, but is non-linear with variable stiffness over the time of the collision. Therefore, cannot be easily divided into large linear segments.
Luckily, for penalty-based simulation we can instead rely a different technique for reducing the gradient error. The simple solution is to reduce the stepsize, which does not work in the ideal clastic case, seen in Fig. 4. From ODE perspective, this works because the penalty- based ODE is smooth, allowing us to continuously control the integration error by reducing the stepsize. Unfortunately, simply reducing the step size is not practical solution, as it necessitates extremely small steps that substantially increase the computational and memory demands of gradient computation. This trade-off raises a critical question: Can we retain correct gradients of realistic contacts while maintaining practical simulation speeds?
3.1
ADAPTIVE STEPSIZE INTEGRATION: NUMERICAL PRECISION ON DEMAND
standard method for integrating ODEs with variable stiffness is adaptive integration. The idea behind adaptive integration elegant: Two numerical integrators of different orders compute the next state. Their difference provides an estimate of the error. If the error is smaller than a given threshold, the step is accepted; otherwise, the step is rejected and the procedure is repeated with different stepsize chosen by a feedback controller. For further details on the rich history of adaptive stepsize integration, see e.g. Hairer al. (2008); Hairer & Wanner (2002); Söderlind (2002; 2003).
We use Diffrrax (Kidger, 2021) for efficient numerical integration in JAX, taking advantage of its solver flexibility and multiple backpropagation modes. Notably, detailed in Appendix C.2, we
4

---

Published as a conference paper at ICLR 2026
hrr)) position-leve
goal
reference acoclcration
loss
F
Fz
d(r).
impedance
Vanilla MIX CFD
Gradient
0.040
0.2
Wec
Staight-though
Gradient
LoSs
mc
0.0
do
0.025
de
0.2
0
-4
-2
0
4
F2
Figure 6: Contacts from dis- tance (CFD): let MuJoCo create small contact forces be- tween non-colliding objects, reference acceleration h(r) and impedance d(r) are ad-
Applying contact forces for in the forward pass of the simulation causes robot to hover. Bottom: The straight- through-trick is used to replace the original MJX derivative with the derivative of MJX CFD, evalu-
Figure 7: Top:
Figure 8: Billiard simulation. Top: Force acts on the white ball affecting the loss. Bottom: Despite the loss if the balls do not collide, DiffMJX with CFD provides informa- tive non-zero gradients.
derivative being
justed to be nonzero for pos- itive signed distances > O.
ated at the unaltered trajectory.
devoted substantial effort seamlessly integrating quaternions and stateful actuators, which enables seamless compatibility between Diffrax and MJX while ensuring efficient adaptive integration and backwards compatibility with other MuJoCo libraries.
Testing our implementation on the cube bounce toy example, we observe in Fig. 5 that adaptive integration reduces the error in the loss and gradient by multiple orders of magnitude given the same computational budget;
Emor 10
10
detailed analysis is in Fig. 17 in the Appendix. Importantly, the adaptive stepsize selection in Diffrax independent across parallelized environments. Therefore, when vmapping across environments, contacts in one simulation do not slow down other simulations. For experimental verification, see Fig. 20 in Appendix D.
10-3
10-5
10-7
10-1
109
10!
102
Resolving problems of Collision Detection. Using an adaptive integrator with MJX eliminates oscillations in the bounce example. However, gradients for some object primitives (capsule, cylinder, box) experience gradient artifacts due non-differentiable opera- tions in the collision detector arising from discrete case distinctions. We smoothed them with standard proxies, finally leading the sults in the bottom row of Fig. 3, where analytical gradients nearly match central differences. Henceforth, we refer to MJX with the Diffrax integrator and smoothed collision detection as DiffMJX.
Runtime (s)
Pareto front of gra- dient error vs forward runtime and gradient error vs back- ward runtime (4) for standard semi-implicit Euler in MJX (black) and adaptive integra- tion in DiffMJX (blue) on the cube bounce toy-example.
Figure 5:
4 CONTACTS FROM DISTANCE WITH STRAIGHHT-HROUGH ESTIMATION
While adaptive integration improves gradient accuracy, we now shift attention computing infor- mative gradients between objects that are not in contact. To illustrate why the computation of such gradients is of fundamental importance for robot learning consider the following example.
Example - Billiard shot: A billiard table is set up as shown in Fig. (top). At the first timestep, a force F is exerted the white ball such that it may hit the black ball. The optimization objective is the distance L between the black ball and the target position. If the balls collide, MJX with
5

---

yields informative gradients Published as a conference paper at ICLR 2026
adaptive integration yields non-zero gradients V.L. However, if F does not cause the balls to touch, then F is zero and therefore uninformative for optimization.
In the following, we propose contacts from distance (CFD), a method for computing contact forces for positive signed distances in the gradient computation of penalty-based simulation that yields informative gradients even if objects are not in contact. This is accomplished by (i) applying artificial contact forces between non-colliding objects, and (ti) using artificial forces only in the gradient computation to maintain simulation realism.
Creating artificial contact forces. How to generate artificial contact forces penalty-based simulator depends on the respective contact model. Below, concrete example, we propose method tailored to MuJoCo. As discussed in Section 2. the magnitude of contact forces is determined by the impedance d(r)) and position-leve' reference acceleration h(r). To enable the solver apply CFD, d(r) is augmented as shown in Fig. 6. Here, d(r) remains unaltered for r < 0 and is extended by an additional spline for O. This continuation is parametrized by solimp-CFD parameters (de, do, We: me, Pe). By default, the curve smoothly continues MuJoCo's impedance at do and tapers off de 0 to ensure smooth differentiability. The CFD-width we specifies the distance for which artificial contact forces are generated. Moreover, we soften the reference acceleration h(r) by replacing the ReLU function on the signed distance with softplus (Fig. 6). This yields modified contact forces SCrD-
Designing a surrogate gradient estimator. Naively adding CFD to a simulation produces non- physical behaviors. As shown in Fig. 7 (top), the artificial contact forces would cause quadruped to hover above the ground as a soft foam mat of thickness w, had been placed on the surface. As significantly altering simulation realism is not option, are faced with the question: Can CFD be used to obtain informative contact gradients without affecting simulation realism?
To positively answer this question, resort the straight-through-trick on the ODE level:
a(D) = sg( (F, x(1))) + (()) - (L,,
where sg is the stop-gradient operator of automatic differentiation library. Here F' denotes the original ODE obtained from MJX DiffMJX and f denotes the ODE using CFD via (CFD M + SCED). In JAX this reads as:
jax. lax Import prop_gradient
2 def Forward im: Model, d: Dats) Cata:
dajjx- forward(m, d, efd-Baise) Compute system acceleration without d (mr d, efd=true) $ compute system acceleraticn using CED
grsd_replace_jn lambda x_mix, x_dfd stop_gredient (x_mjx) x_efd stop_gradient (x_efd) return www.trre.mmppp.. replaceen,. d_jjs, e # Rerouze gradiens computation
As illustrated in Fig. 7 (bottom), the above code ensures that the forward pass uses (()))
whereas the backward pass deploys D[.00] (1, x(()). Crucially, the derivatives are evaluated at the unmodified forward trajectory r=(t). While approach requires the simulator's forward pass to be computed twice, the gradient is only evaluated once. As the gradient computation dominates the computational cost, CFD forms a practical method for improving contact gradients. Revisiting the billiard toss example (Fig. 8). using CFD with the truight-through while keeping the loss unaltered.
Note that the traight-through trick is inspired by prior work on al., 2013); see also the related work in Appendix DiffMIX and CFD introduce easy-to-use tuning knobs in MJX to control gradient quality. Appendix B explains how to tune these knobs. In Appendix C.2, we evaluate the teaighh-through "Discretize-then--ppiiiz"" and "Opiimize-theen-diseriii
estimators (Bengio
with different autodiff techniques, including

---

Published as a conference paper at ICLR 2026
MJX
DiffMjX
Experiment setup
CUD6 width in mm
no CFD CFD
Fx+++
F2+
120
real hall
real hall
100
L
80
simulated ball
simalated ball
50
gradient stop gradient stop
Figure 10: Left: If the initial ball radius p is set too small, then the gradient V,L re- sulting from a mee-step-ahead prediction zero. Right: Large initial constraint pene- tration results in large gradients Va,
Figure 9: Left: Estimation of a cube's side length in MJX via gradient descent using multi-step ahead pre- dictions. Right: Experimental setup for collecting cube toss data. Image adapted from Pfrommer al. (2021).
Before training
After training
Figure 11: Comparison between real-world cube tosses (red) and DiffMJX cube simulations (black).
EVALUATION
In the following, demonstrate the use of DiftMIX and CFD for learning physics parameters and computing robot control actions using gradient descent with simulator gradients.
5.1 PARAMETER IDENTIFICATION
Parameter identification in the presence of hard contacts remains a laborious task. If contacts are hard, learning the dynamics of a cube requires impractical amounts data for "naive" neural network regression (Parmar et al., 2021). In comparison, penalty-based -baeed simulators but the lack of correct gradients hinders efficient parameter estimation (Acosta et al., 2022). Therefore, recent work introduced intricate analytical pipelines for cube geometry estimation (Pfrommer et al., 2021; Bianchini et al., 2023) and graph-based networks for learning contact dynamics (Allen et al., 2023). In what follows, we the same real-world data used in Pfrommer al. (2021); Bianchini et al. (2023); Allen et al. (2023). We demonstrate that DiffMJXX with CDF enables simulator parameter estimation via standard gradient-based optimization.
capture hard contacts,
Dataset and training setup. We use the Contactnets dataset (Pfrommer al., 2021) which consists of 550 trajectories of a 10 cm acrylic cube that has been repeatedly tossed onto wooden table. For training, trajectories are split into segments length five such that the simulator is tasked unroll four future steps starting from the initial state. Each segment and its prediction are fed to an L2 loss whose gradient used for gradient- based optimization using Adam (Kingma Ba, 2015). For systems with stiff dynamics, we favor multi-step-alee step-ahead predictions over one-step-ahead predictions, as they capture the cumulative effects of prediction errors over time. This setup enables fairer analysis of MJX without CFD, as even for too small side length estimates as illustrated in Fig. 10, future state predictions can make contact to inform the optimization.
Training results. The training results are shown in Fig. 9. Our version of MJX with smooth collision detection and MJX with CFD both achieve good estimation results with error of around 5% relative to the ground truth. If the side length is initialized at 60 mm or 140 mm, training either stalls fully convergence is severely limited for MJX. The incorporation of CFD into MJX addresses convergence issues arising from poor initial parameters, while the integration of adaptive integration via DiffMIX significantly enhances estimation accuracy. DiffMJXX improves estimation accuracy
7

---

Published as a conference paper at ICLR 2026
In-hand manipulation Bionic tennis
Sampling
Cost
10
-2
Gradients with CFD
goal
M
10
Sampling (04/256/1024)
Gradients (no CFD)
Gradients (CFD)
10 -6
Time
1000 O
1000
Simulation Step
Figure 12: Autodiff-driven MPC with and without CFD on the in-hand manipulation and bionic tennis tasks (Fig. 13). In both tasks, only the distance between the ball and the respective goal is used as cost. Left: Simulation cost evolution of gradient-based MPC, with and without contacts from distance (CFD), vs sampling-based MPC. The number for sampling indicates the number of samples used per planning step. Sampling has difficulties solving the dexterous in-hand manipulation task; gradients without CFD cannot solve the bionic tennis task. Right: Rendering of sampling-ban based MPC (1024 samples) vs gradient- based MPC with CFD on the in-hand manipulation task. The goal swap the balls, with the cost computed as L2 distance of the ball centers to positions fixed the frame of the hand. The MyoHand model is actuated by 39 muscle-tendon units.
by dynamically adjusting the time steps during collisions, thereby mitigating time discretization errors. Further details and additional experiments in which also the contact parameters were identified are provided in E.2. Fig. 11 provides a comparison of DiffMIX's predictions after estimating contact and geometry parameters. To the best of our knowledge, we are the first to demonstrate parameter estimation of real-world cube dynamics using an automatically differentiable penalty-based simulator. While this represents a promising step forward, further experimentation is necessary to fully characterize the scope and limitations of this approach.
5.2 MODEL PREDICTIVE CONTROI.
Next, we conduct experiments on gradient- based model-predictive control. At every plan step in the MPC loop, refine a sequence of controls over 256-step horizon. In the gradient-based planner, we compute gradients by backpropagating the differentiable cost computed on the rollout of the current plan through the MJX simulator. The plan is then iteratively optimized using the Adam optimizer with a learning rate of 0.01 for 32 iterations. Finally, the resulting plan is executed for 16 steps in simulation, after which the planning procedure is repeated with the previous plan warm start. As baseline, we include a version of the predictive sampling planner from Mujoco MPC (Howell al., 2022), which every plan step samples k {64, 256, 1024) trajectories, and executes the lowest-cost plan. For enabling a fair comparison, we significantly improved the performance of this planner for muscular systems by resorting to brown noise for sampling (Pinneri al., 2020).
Models. As physical systems, we resort state-of-the-an muscle-tendon models provided by MyoSuite (Caggiano et al., 2022; Wang et al., 2022). Models include the MyoHand (Fig. 12, right) adapted from the MyoChallenge 2022, which is comprised of 29 bones, 23 joints, and 39 muscle- tendon units. We also use a bionic model (Fig. 13) modified from the MyoChallenge 2024, which is comprised of the MyoArm with 27 degrees of freedom and 63 muscle-tendon units, and the simulated modular prosthetic limb with 26 degrees of freedom and 17 motor control units.
Dexterous in-hand manipulation. First, we consider an in-hand manipulation task, where the goal is to swap two balls in the MyoHand. The cost is given by the Euclidean distance between each of the balls and the respective target location, fixed in the frame of the hand. Note that the muscle actuator implementation caused gradient errors, which we corrected using smooth functions surrogates. The results reported in Fig. 12. We find that gradient-basee based MPC can reliably based
solve this task, in contrast to the sampling- based planner. Overparameterizzaion in the muscle-tendon model with least two muscles per joint actually benefits the gradient-bact planner by helping
8

---

Published as a conference paper at ICLR 2026
start
catch
hit
goal
time
move to gool
Figure 13: Autodiff-driven MPC with CFD on the bionic tennis task. Task completion requires the racket to deflect the ball towards the MyoArm with 63 muscle-tendon actuators, which then catches the ball and moves it to the goal position. Only the distance between ball and goal is used as cost.
escape local minima, similar to its role in optimizing overparametrized neural networks. In contrast, RL and sampling-basee based planners struggle with scaling in overparametrized higher-dimensional systems (Schumacher al., 2023). Thus, first-order methods using differentiable simulation should be able to tackle more complex control problems.
Notably, this task does not require contacts from distance because hand-ball interactions are frequent due to gravity. Moreover, we identify two crucial components of the gradient- bbsed MPC loop: First is gradient clipping, which is important the scale of gradients changes massively in the presence of contacts, as illustrated in Fig. 10 (right). This technique has also been reported to be effective in previous works on differentiable simulation (Xu et al., 2022: Georgiev et al., 2024). Second, we store the rollout cost of all gradient iterations and select the one with minimal cost. This is important as the cost landscape is highly non-convex, which reflected in the non- non-monotonic cost evolution between the iterations of planning step.
Using CFD to solve complex control tasks with minimal task supervision. Next, we test more complex custom tennis task on the bionic model. The task is to move a ball that is initially moving sideways to target location below. This can be achieved by bouncing it back using a racket that is welded to the prosthetic hand, and then catching it at the target location with the musele hand. In this task, the only cost supervision is again the Euclidean distance of the ball to the target, the complicated sequential movement has to be discovered purely from this signal.
Bionic tennis:
We report our findings in Fig. 12 (left), see Fig. 13 for a rendering. By design, the task initialization such that the ball misses both hands. hence this task is not solvable by purely gradient- gradient-based MPC using vanilla MJX. On the other hand, we observe that adding the CFD mechanism allows solving this task. The sampling- based planner Initially, bouncing the ball back to the target only requires controlling the prosthetic arm, which is relatively dineennnii.n Hence, the sampling-batt based planner achieves this part easily. However, as seen in the in-hand manipulation task, it struggles with precise control of the high-dimensional MyoHand, leading to sub-optimal results in balancing the ball at the goal position.
strong baseline in this task and gets close to solving
Cube reorientation. Finally, tackle another manipulation task of reorienting cube to desired target orientation, while staying at a target position. We include both an in-hand version and an on-the-table version of this task.
We report our findings in Fig. 14. We observe that in the in-hand reorientation task, the sampling- based planner fails, while the gradient-based approach solves it. Similar to the ball-swap task, CFD does not make a big difference here, as contacts are automatically discovered from gravity. In the on-the-table reorientation task, only the gradient-based planner with CFD activated solves the task perfectly. When deactivating CFD or using sampling, the planner still manages to reorient the cube satisfactory level, but fails keep in the target location. Here, gravity alone does not result in all fingers making contact with the cuber, therefore CFD allows ligher-fidelity control in comparison to CFD disabled.
9

---

Published as a conference paper at ICLR 2026
10
Sampling
In-hand Reorientation
19'
Gradients with CFD
10
Cost
200
Simulation Step
goal
(n-the--d -o--floo Regirentation
19.
Sampling
10.
10
Sampling (64/256/1024) Gradients (no CFD)
Gradients with CFD
Gradlents (CFD)
10:
Cost
200
Simulation Step
400
0410
1000
Time
goal
Figure 14: Autodiff-driven MPC with and without CFD the in-hand and on-the-floor cube reorientation tasks. In both tasks, the loss consists of the L2 distance between the cube center and the target location, well squared loss on the difference of quaternions describing the cube and target orientation. Top: Simulation cost evolution of gradient-based MPC, with and without contacts from distance (CFD), vs sampling-basee MPC. The number for sampling indicates the number of samples used per planning step. Sampling has difficulties solving the the in-hand reorientation task; only gradients with CFD can solve the on-the-table reorientation task perfectly. Right: Rendering of sampling- sampling-based based MPC (1024 samples) vs gradient-basse on-the on-the-table reorientation task. The goal is to reorient the cubes while staying at the target location (green sphere). The MyoHand model is actuated by 39 muscle-tendon units.
MPC with CFD on the in-hand and
CONCLUSION
In this work, we tackle two standing challenges of computing gradients of penalty- simulations with hard contacts, namely gradient errors due to time discretization and zero contact gradients between non-colliding objects. A ID toy example illustrates how gradient errors are consequence of time discretization and be mitigated by reducing the integration stepsize. While reducing the stepsize improves gradient accuracy, it comes at the cost of increased simulation time and GPU memory usage. To address these limitations, we introduce DifiMJX, JAX library that incorporates the adaptive integration library Diffrax into MuJoCo XLA. The use of adaptive time integration enables the reduction of discretization errors while significantly reducing memory overhead through checkpointing. Complementing DiffMIX, we propose Contacts from Distance (CFD), an exten- sion of penalty-based contact resolution that improves gradient utility by introducing small virtual contact forces between near-interacting bodies, without affecting the forward simulation due to the We validated our contributions where DiffMJXX accurately recovered cube's geometric parameters via gradient-based optimization. Beyond this, we present the first evidence that automatic differentiation through CFD can outperform sampling-ber based predictive planning musculoskeletal manipulation tasks. We hope these results will guide future work on mitigating discretization errors in differentiable simulators and establish DiffMJX with CFD as helpful toolbox for gradient-bard based based robot learning.
real-world system identification study,
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGEMENTS
The authors thank Onur Beker, Thomas Rupf, and Nico Gürtler for providing valuable feedback.
We thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Anselm Paulus.
This work was supported by the ERC 101045454 REAL-RRL and the German Federal Ministry of Education and Research (BMBF) through the Tübingen AT Center (FKZ: DJSS180399 Georg Martius is a member of the Machine Learning Cluster of Excellence, EXC number 2064/1 Project number 390727645.
The work on this paper was supported by the Czech Science Foundation (GACR) grant no. 26-23981S.
REFERENCES
Brian Acosta, William Yang, and Michael Posa. Validating robotics simulators on real-world impacts.
IEEE Robotics and Automation Letters, 7(3):6471-6478 2022.
Kelsey R Allen, Tatiana Lopez Guevara, Yulia Rubanova, Kim Stachenfeld, Alvaro Sanchez-Gonzake Peter Battaglia, and Tobias Pfaff. Graph network simulators can learn discontinuous, rigid contact dynamics. In Karen Liu, Dana Kulic, and Jeff Ichnowski (eds.). Proceedings of The 6th Conference on Robot Learning, volume 205 of Proceedings of Machine Learning Research, pp. 1157-1167. PMLR, 14-18 Dec 2023.
OpenAl:
Marcin Andrychowicz. Bowen Baker, Maciek Chociej, Rafal Józefowicz, Bob McGrew, Jakub Pachocki, Arthur Petron, Matthias Plappert, Glenn Powell, Alex Ray, Jonas Schneider, Szymon Sidor, Josh Tobin, Peter Welinder, Lilian Weng, and Wojeiech Zaremba. Learning dexterous in-hand manipulation. The International Journal of Robotics Research, 39(1):3-20, January 2020. ISSN 0278-3649. doi: 10.1177/0278364499887774
Joachim Baumgarte. Stabilization of constraints and integrals of motion in dynamical systems. Computer methods in applied mechanics and engineering, 1(1):1-16, 1972.
Onur Beker, Nico Gürtler, Ji Shi, René Geist, Amirreza Razmjoo, Georg Martius, and Sylvain Calinon. A smooth analytical formulation of collision detection and rigid body dynamics with contact. arXiv preprint arXiv:2503. 17366, 2025.
Yoshua Bengio, Nicholas Léonard, and Aaron C. Courville. Estimating propagating gradients through stochastic neurons for conditional computation. CoRR, abs/1308.344 3432, 2013. URL
http:://rxiii org/ abs /1308 3432.
Bibit Bianchini, Mathew Halm, and Michael Posa. Simultaneous learning of contact and continuous dynamics. In Conference on Robot Learning, pp. 3966-3978. PMLR, 2023.
Frédéric Bonnans and Julien Laurent-Varin Computation of order conditions for symplectic partitioned runge-kutta schemes with application to optimal control. Numer. Math., 103(1):1-10, March 2006. ISSN 0029-599X. doi: 1010000002211000006------- URL https:///do...... 10.10071 00211 0065---66....
Karim Bouyarmane, Adrien Escande, Florent Lamiraux, and Abderrahmane Kheddar, Potential field guide for humanoid multicontacts acyclic motion planning. In 2009 IEEE International Conference on Robotics and Automation, pp. 1165 170. IEEE, 2009.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Miin and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL
http:// http://ttt./..t//// thub com hamggooglee Meroolll//o jax.
Romain Brégier. Deep regression on manifolds: 3d rotation case study. In 2021 International Conference 3D Vision (3DV), pp. 166-174, IEEE, 2021.
11

---

Published as a conference paper at ICLR 2026
Vittorio Caggiano, Huawei Wang, Guillaume Durandau, Massimo Sartori, and Vikash Kumar. Myosuite: A contact-rich simulation suite for musculoskeletal motor control. In Roya Firoozi, Negar Mehr, Esen Yel, Rika Antonova, Jeannette Bohg. Mac Schwager, and Mykel J. Kochenderfer (eds.), Learning for Dynamics and Control Conference, LADC 2022, 23-24 June 2022, Stanford University, Stanford, CA, USA, volume 168 of Proceedings of Machine Learning Research, pp. 492-507. PMLR, 2022.
M. Calvo, S. González Conzález-Pintoo and J.I. Montijano. Global error estimation based on the tolerance proportionality for some adaptive runge-kutta codes. Journal of Computational and Applied Mathematics, 218 218(():329-35))
2008. ISSN 0377-0427. doi: https://doi.org/1.... 2007.02.034. URL https: / /www. sciencedirect, com/ wwtpscceeecceeeeeeeeeeeeeeecceee aaraicceee /pii// 10377042707001 7000200.. The Proceedings of the Twelfth International Congress Compu- tational and Applied Mathematics.
Erin Catto. Box2d is a 2d physics engine for games. 2018.
Erin Catto et al. Modeling and solving constraints. In Game Developers Conference, 2009.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary
differential equations. Advances in Neural Information Processing Systems, 2018.
Erwin Coumans and Yunfei Bai. Pybullet, a python module for physics simulation for games, robotics and machine learning. http: /tppyyyybllleteeeeee org, 2016.
DeepMind, Igor Babuschkin, Kate Baumli, Alison Bell, Surya Bhupatiraju. Jake Bruce, Peter Buchlovsky, David Budden, Trevor Cai, Aidan Clark, Ivo Danihelka, Antoine Dedieu, Clau- dio Fantacci, Jonathan Godwin, Chris Jones, Ross Hemsley, Tom Hennigan, Matteo Hessel, Shaobo Hou, Steven Kapturowski, Thomas Keck, Iurii Kemaex, Michael King, Markus Kunesch, Lena Martens, Hamza Merzic, Vladimir Mikulik, Tamara Norman, George Papamakarios, John Quan, Roman Ring, Francisco Ruiz, Alvaro Sanchez, Laurent Sartran, Rosalia Schneider, Eren Sezener, Stephen Spencer, Srivatsan Srinivasan, Miloš Stanojević, Wojeiech Stokowiee, Luyu Wang, Guangyao Zhou, and Fabio Viola. The DeepMind JAX Ecosystem, 2020. URL http::/githus com/ google denple--deeinddd.
J. R. Dormand and P. J. Prince. A family of embedded Runge-Kutta formulae. J. Comp. Appl. Math, 6:19-26, 1980.
Evan Drumwright and Dylan A Shell. Modeling contact friction and joint friction in dynamic robotic simulation using the principle of maximum dissipation. In Algorithmic Foundations of Robotics IX: Selected Contributions of the Ninth International Workshop on the Algorithmic Foundations of Robotics, pp. 249-266. Springer, 2011.
John C Duchi, Peter L Bartlett, and Martin J Wainwright. Randomized smoothing for stochastic optimization. SIAM Journal on Optimization, 22(2):674-701, 2012.
Ryan Elandt, Evan Drumwright, Michael Sherman, and Andy Ruina. A pressure field model for fast, robust approximation of net contact force and moment between nominally rigid objects. In 2019 JEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. IEEE, 2019.
Roy Featherstone. Rigid body dynamics algorithms. 2014.
C. Daniel Freeman, Erik Frey, Anton Raichuk, Sertan Girgin, Igor Mordatch, and Olivier Bachem. Brax differentiable physies engine for large scale rigid body simulation, 2021. URL http: //github. wtt.coooooooo/oe
brax.
Carl Friedrich GauB, Über ein neues allgemeines Grundgesetz der Mechanik. Journal für die reine und angewandte Mathematik, 4:232-235, 1829.
Andreas René Geist, Jonas Frey, Mikel Zhobro, Anna Levina, and Georg Martius. Learning with 3D rotations, hitchhiker's guide to SO(3). In Proceedings of the 41st International Conference Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 15331-15350. PMLR, 21-27 Jul 2024.
12

---

Published as a conference paper at ICLR 2026
Ignat Georgiev, Krishnan Srinivasan, Jie Xu, Eric Heiden, and Animesh Garg. Adaptive horizon actor-critic for policy learning in contact-rich differentiable simulation. In Forty-first International Conference on Machine Learning, ICML 2024, Vienna, Austria, July 21-27, 2024. OpenReview.net, 2024.
Christian Gumbsch, Martin V Butz, and Georg Martius. Sparsely changing latent states for prediction and planning in partially observable domains. Advances in Neural Information Processing Systems, 34:17518 11-11518-17531 2021.
William Hager. Runge-kutta methods in optimal control and the transformed adjoint system. Nu- merische Mathematik, 87, 02 2001. doi: 10. 1007/00//021 10000178.
E. Hairer and G. Wanner. Solving Onlinary Differential Equations II Stiff and Differential-Algerriii Problems. Springer, Berlin, second revised edition edition, 2002.
E. Hairer, S.P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations 1 Nonstiff Problems. Springer, Berlin, second revised edition edition, 2008
Eric Heiden, David Millard, Erwin Coumans, Yizhou Sheng, and Gaurav S Sukhatme. Neural- Sim: Augmenting differentiable simulators with neural networks. In Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), 2021.
Michael Hinze, René Pinnau, Michael Ulbrich, and Stefan Ulbrich. Optimization with pde con- straints. In Mathematical Modelling, 2008. URL https: / /tpp//ap.p /api. semanticscholar.com org/
CorpusID: 1278773.
Coşku Can Horuz, Geoffrey Kasenbacher, Saya Higuchi, Sebastian Kairat, Jendrik Stoltz, Moritz Pesl, Bernhard A Moser, Christoph Linse, Thomas Martinetz, and Sebastian Otte. The resurrection of the relu. arXiv preprint arXiv:2505.22074 2025.
Taylor Howell. Mujoco warp (mjwarp). https: / / g sintthuub.. com http//google-deeemmm.e mujoco_ warp, March 2025. NVIDIA GPU Technology Conference (GTC).
Taylor Howell, Nimrod Gileadi, Saran Tunyasuvunakool. Kevin Zakka, Tom Erez, and Yuval Tassa. Predictive Sampling: Real-time Behaviour Synthesis with MuJoCo. dec 2022.
Yuanming Hu, Tzu-Mao Li, Luke Anderson, Jonathan Ragan-Kelley, and Frédo Durand. Taichi: a language for high-performance computation on spatially sparse data structures. ACM Transactions on Graphics (TOG). 38(6):201, 2019a.
Yuanming Hu, Jiancheng Liu, Andrew Spielberg, Joshua Tenenbaum, William T Freeman, Jiajun Wu, Daniela Rus, and Wojciech Matusik. Chainqueen: A real-time differentiable physical simulator for soft robotics. 2019 International conference on robotics and automation (ICRA), pp. 6265- 6271. IEEE, 2019b.
Yuanming Hu, Luke Anderson, Tzu-Mao Li, Qi Sun, Nathan Carr, Jonathan Ragan- and Frédo Durand. Difftaichi: Differentiable programming for physical simulation. In International Conference Learning Representations, 2020.
Danny M Kaufman, Shinjiro Sueda, Doug I. James, and Dinesh K Pai. Staggered projections for frictional contact multibody systems. In ACM SIGGRAPH Asia 2008 papers, pp- 1-11. 2008
Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
Patrick Kidger and Cristian Garcia. Equinox: neural networks in JAX via callable PyTrees and filtered transformations. Differentiable Programming workshop at Neural Information Processing Systems 2021, 2021.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In Yoshua Bengio and Yann LeCun (eds.), 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
Vittorio La Barbera, Steven Bohez, Leonard Hasenclever, Yuval Tassa, and John R Hutchinson. Motion tracking with muscles: Predictive control of a parametric musculoskeletal canine model. arXiv preprint arXiv:2506.23700 2025.
13

---

Published as a conference paper at ICLR 2026
Joonho Lee, Jemin Hwangbo, Lorenz Wellhausen, Vladlen Koltun, and Marco Hutter. Learning quadrupedal locomotion over challenging terrain. Science Robotics, 5(47):eabc5986, October 2020. doi: 10.1 10.1112/siirobbttcssa::3333
C. Li, M. Vlastelica, S. Blaes, J. Frey, F. Grimminger, and G. Martius. Learning agile skills via adversarial imitation of rough partial demonstrations. Proceedings of the 6th Conference on Robot Learning (CoRL), December 2022.
Jacky Liang, Viktor Makoviychuk, Ankur Handa, Nuttapong Chentanez, Miles Macklin, and Dieter Fox. Gpu-accelerated robotic simulation for distributed reinforcement learning, 2018.
Min Liu, Gang Yang, Siyuan Luo, and Lin Shao. Softmac: Differentiable soft body simulation with forecast-based contact model and two-way coupling with articulated rigid bodies and clothes. In 2024 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 12008- 12015. IEEE, 2024.
Miles Macklin. Warp: A high-performance python framework for gpu simulation and graphics. https: /tapps://gtttiie thub. www/moimmiiciieeeeeeoo /warp, March 2022. NVIDIA GPU Technology Confer- ence (GTC).
Igor Mordatch, Zoran Popović, and Emanuel Todorov. Contact- optimization for hand manipulation. In Proceedings of the ACM SIGGRAPH/Eunogaphiis symposium on computer animation, pp. 137-144, 2012a.
Igor Mordatch, Emanuel Todorov, and Zoran Popović. Discovery of complex behaviors through contact-invariant optimization. ACM Transactions on Graphics (ToG), 31(4):1-8, 2012b.
Rhys Newbury, Jack Collins, Kerry He, Jiahe Pan, Ingmar Posner, David Howard, and Akansel Cosgun. A Review of Differentiable Simulators. IEEE Access, 12:97581-97604, 2024. ISSN 2169-3536.
Tao Pang, HJ Terry Suh, Lujie Yang, and Russ Tedrake. Global planning for contact-rich manipulation via local smoothing of quasi-dynamic contact models. IEEE Transactions on robotics, 39(6): 4691 4711, 2023.
Mihir Parmar, Mathew Halm, and Michael Posa. Fundamental challenges in deep learning for stiff contact dynamics. In 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5181-5188. IEEE, 2021.
Samuel Pfrommer, Mathew Halm, and Michael Posa. Contactnets: Learning discontinuous contact dynamics with smooth, implicit representations. In Conference on Robot Learning, pp. 2279-2291.
PMLR, 2021.
Cristina Pinneri, Shambhuraj Sawant, Sebastian Blaes, Jan Achterhold, Joerg Stueckler, Michal Rolinek, and Georg Martius. Sample-efficient cross-entropy method for real-time planning. In Conference on Robot Learning (CoRL) 2020, volume 155, pp. 1049-1065. PMLR, 2020. URL https::/com corl:://ccoccoooooo/ ..f.fiiiie github. i paper_2111
Ilija Radosavovic, Tete Xiao, Bike Zhang, Trevor Darrell, Jitendra Malik, and Koushil Sreenath. Real-world humanoid locomotion with reinforcement learning. Science Robotics, (89))aa April 2024. doi: 10.1 12//scinohooicssadi959999
Russ Tedrake, and the Drake Development Team. Drake: Model-based design and verification for robotics, 2019. URL
wwt.mit.......m................. edu.
Subham Sekhar Sahoo, Anselm Paulus, Marin Vlastelica, Vít Musil, Volodymyr Kuleshov, and Georg Martius. Backpropagation through combinatorial algorithms: Identity with projection works. In The Eleventh International Conference Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023, OpenReview.net, 2023, URL https: : bttps::/openzee//// /ttp///pwwweeeeeeeeeeeeee net for ummm? dd=JZMR7 7 70299.
14

---

Published as a conference paper at ICLR 2026
Adrian Sandu. On the properties of runge-kutta discrete adjoints. In Proceedings of the 6th International Conference on Computational Science Volume Part IV. ICCS 06, pp. 550-557, Berlin, Heidelberg, 2006. Springer- Verlag. ISBN 3540343857. doi: 10.1007/111585499 76. URL
https: /doi. .rrg 0 1007 11154717759999
Pierre Schumacher, Daniel F.B. Haeufle, Dieter Büchler, Syn Schmitt, and Georg Martius. Dep-rl: Embodied exploration for reinforcement learning in overactuated and musculoskeletal systems. In The Eleventh International Conference on Learning Representations (ICLR), May 2023.
Clemens Schwarke, Victor Klemm, Jesus T Jean-Pierre Sleiman, and Marco Hutter. Learn- ing quadrupedal locomotion via differentiable simulation. arXiv preprint arXiv:2404.02887 2024.
Gustaf Söderlind. Automatic control and adaptive time-stepping. Numerical Algorithms, 31:281-310, 2002.
Gustaf Söderlind. Digital Filters in Adaptive Time-Stepping. ACM Transactions Mathematical Software, 20(1):1-26, 2003.
Yunlong Song, Sang bae Kim, and Davide Scaramuzza. Learning quadruped locomotion using differentiable simulation. In 8th Annual Conference on Robot Learning, 2024.
Hyung Ju Terry Suh, Tao Pang, and Russ Tedrake. Bundled gradients through contact via randomized smoothing. IEEE Robotics and Automation Letters, 7(2):4000-4007, 2022.
Jie Tan, Tingnan Zhang, Erwin Coumans, Atil Iscen. Yunfei Bai, Danijar Hafner, Steven Bohez, and Vincent Vanhoucke. Sim-to-real: Learning agile locomotion for quadruped robots. In Proceedings of Robotics: Science and Systems, Pittsburgh, Pennsylvania, June 2018. doi: 10. 10.15607/RSS.2016 XIV.010.
Yuval Tassa and Emo Todorov. Stochastic complementarity for local control of discontinuous dynamics. In Robotics: Science and Systems, volume 6, 2010.
A. Howell Taylor, Simon Le Cleac'h, Zico Kolter, Mac Schwager, and Zachary Manchester. Dojo: differentiable simulator for robotics. arXiv preprint arXiv:2203.0000 000006 2022.
Emanuel Todorov. A convex, smooth and invertible contact model for trajectory optimization. In 2011 IEEE International Conference on Robotics and Automation, pp. 1071-1076, May 2011.
ISSN: 1050-4729.
Emanuel Todorov. Convex and analyticall--ineriible dynamics with contacts and constraints: Theory and implementation mujoco. 2014 IEEE International Conference on Robotics and Automation (ICRA), pp. 6054-6061. IEEE, 2014.
Emanuel Todorov, Tom Erez, and Yuval Tassa. MuJoCo: A physics engine for model-based control. In 2012 JEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 5026-5033, October 2012. ISSN: 2153-0866.
Tsitouras Runge-kutta pairs of order 5 (4) satisfying only the first column simplifying assumption. Computers & Mathematics with Applications, 62(2):770-775, 2011.
Andrea Walther. Automatic differentiation of explicit runge-kutta methods for optimal control. Com- putattioall Optimization and Applications, 36:83-108, 11 2007. doi: 10.1007/111599-0000003..7
Huawei Wang, Vittorio Caggiano, Guillaume Durandau, Massimo Sartori, and Vikash Kumar. Myosim: Fast and physiologically realistic mujoco models for musculoskeletal and exoskeletal studies. In 2022 International Conference Robotics and Automation, ICRA 2022, Philadelphia. PA. USA, May 23-27, 2022, pp. 8104-8111. JEEE, 2022.
Keenon Werling, Dalton Omens, Jeongseok Lee, Ioannis Exarchos, and C Karen Liu. Fast and ceature-complete differentiable physics engine for articulated rigid bodies with contact constraints. In Roboties: Science and Systems, 2021.
15

---

Published as a conference paper at ICLR 2026
Jie Xu, Viktor Makoviychuk, Yashraj S. Narang, Fabio Ramos, Wojciech Matusik, Animesh Garg, and Miles Macklin. Accelerated policy learning with parallel differentiable simulation. In The Tenth International Conference Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net, 2022.
Jijie Xu, Tak-Kuen John Koo, and Zexiang Li. Sampling-based finger gaits planning for multifingered robotic hand. Autonomous Robots, 28(4):385 4022, 2010.
Gang Yang. Siyuan Luo, Yunhai Feng, Zhixin Sun, Chenrui Tie, and Lin Shao. Jade: differentiable physics engine for articulated rigid bodies with intersection- frictional contact. In 2024 IEEE International Conference on Robotics and Automation (ICRA), pp. 16915-16922. IEEE, 2024.
Yunbo Zhang, Deepak Gopinath, Yuting Ye, Jessica Hodgins, Greg Turk, and Jungdam Won. Simula- tion and retargeting of complex multi-character interactions. In ACM SIGGRAPH 2023 Conference Proceedings, pp. 1-11, 2023.
16

---

Published as a conference paper at ICLR 2026
A RELATED WORK
Differentiable simulators are an active area of research spanning multiple fields of physics, including elastic object and fluids modeling Hu al. (2019b); Macklin (2022); Hu et al. (2020); Liu al. (2024) and modeling rigid body collisions Todorov et al. (2012); Hu et al. (2020); Taylor et al. (2022); Howell (2025) (see Newbury et al. (2024) for a recent overview). As outlined Fig. 15, we focus on robotics simulators involving hard contact collisions, with the hierarchical objectives of (i) accurately simulating dynamics, (ii) computing correct simulator gradients, and (iii) obtaining informative gradients between non-colliding objects.
Goal 1: Realistic simulation
Simulating contact dynamics. Commonly deployed robot simulators can be categorized according to their contact model into impulse-based Catto (2018); Freeman et al. (2021); Hu et al. (2020), camplementary-basee Taylor et al. (2022); Wer- ling et al. (2021); and penalty-based Todorov et al. (2012) approaches. Due to their simplicity and speed, impulse-bag based simulators (Catto et al., 2009), such as those using DiffTaichi (Hu et al., 2019a; 2020), are widely used in game development but remain currently uncommon in robot controller synthe- sis. Complementary-busee contact models compute constraint forces as a solution to a constrained optimization problem ei- ther in the form of nonlinear-com complmentrr problem (NCP)
Salution 1: Exact contact solving
Soluution 2: Compliant contacts with stiff settings
Goal 2: Accurate gradients
Solution 1: Randomized smoothing
Salution differances
Solution 2: Analytic smoothing
Solution 4 (ours): Adaptive integration
Goal 3: Contact invariance
Solution 1: Inter-body distances
(Taylor et al., 2022) or limear-complementary problem (LCP) (Russ Tedrake, 2019; Coumans & Bai, 2016; Heiden et al., 2021; Yang et al., 2024). Exactly solving NCP is an NP-hard problem Kaufman et al. (2008) that is considerably difficult to solve even approximately. The NCP formulation in Dojo is physically more accurate than LCP formulations Taylor et al. (2022), but lacks support for parallel computation limiting its utility for controller synthesis. That said, it is currently not clear which complementary formulation is best suited for robotics.
Contact-invariant optimization
Solution foursh: Contacts from distance
Straight-trough
Figure 15: Overview hierarchi- cal goals that need to be accom- plished to use simulator gradients for robot controller synthesis.
To improve computational efficiency, MuJoCo (Todorov et al., 2012: Todorov, 2014) reformulates the complementary constraint problem as a convex optimization problem. This reformulation builds on prior work on complementarity frre approaches (Todorov, 2011: Drumwright & Shell, 201 1). Briefly, contact forces are computed by solving global optimization problem, in which the ability of a constraint to generate force is modulated by the geometric penetration depth at the contact interface. As an alternative, Drake Russ Tedrake (2019) also adopts a soft patch contact model Elandt et al. (2019), Pang et al. (2023). In turn, these approaches to contact resolution soften the dynamics.
Contact-invariane optimization. To guide an optimizer toward using body collisions for task facilitation, recent works in controller synthesis and path planning incorporate inter-body distances into the optimization objective. In RL, Zhang et al. (2023) added inter-point distances to rewards using softmax function. Alternatively, the algorithms in MuJoCo MPC (Howell et al., 2022) rely on inter-point distances combined with trajectory samples establishing contact to inform the optimization. However, depending on the task hand, setting up distance- loss terms to inform an optimizer about the necessity for specific object collisions can quickly become cumbersome. As an alternative approach, Mordatch et al. (2012a) proposed the framework of "comtact-invarianl optimization" " (CIO) in which a simulation is adjusted to also apply contact forces between non-colliding bodies. While the simulation can apply non-physical contact forces, the optimization is encouraged through the addition of several loss terms to minimize the usage of these forces. In Mordatch et al. (2012b), CIO was extended to in-hand manipulation of objects. Our proposed extension of contacts from distance (CFD) is inspired by CIO. Yet, while CIO applies artificial contact forces in the forward simulation, CFD resorts to the straight-through-tiie is solely used for gradient computation. In turn, CFD does not require the addition of regularization terms to a loss when used for planning. Recently, Beker et al. (2025) proposed soft signed distance fields an alternative geometry representation, while also adjusting collision detection and contact force computation to be inherently soft. This approach which routes at the geometry level of robot simulation can be seen as alternative avenue towards CIO.
to only exert contact forces in a simulation that
17

---

Published as a conference paper at ICLR 2026
CFD's mechanism to create "forces from distance" (FFD) is also conceptually similar Pang et al. (2023). Drawing inspiration from works on creating FFD. this work smoothens the contact model of quasi-static penalty-basee based simulator using log-barrier functions and analyzes FFD in the context of quasi-static dynamics and implicit differentiation. In comparison, CFD uses MuloCo's smooth contact model based polynomial splines to obtain forces from distance. We note that both quasi-static dynamics and implicit differentation would be useful extensions to DiftMMXX and CFD.
Straight-through estimation. The use of the straight-through--riik in CFD is inspired by plethora of works in robotics and beyond (Bengio al., 2013; Sahoo al., 2023). The code underlying the traaighttthrouge trick is reported in JAX (DeepMind et al., 2020) documentation under the collective heading of traight-through estimation. Straight-through estimation has been used Gumbsch al. (2021), where the backward pass of Heaviside function is replaced with identity function to obtain sparse RNNs for model-based RL. Horuz et al. (2025) demonstrate that using a ReLu activation with custom backward function rivals smooth ReLu surrogates. Recently, Song al. (2024) replaced simulator dynamics obtained from IsaacGym (Liang et al., 2018) with single-rigid body dynamics for the gradient computation in the backward pass. This approach is closely related to our work, but unlike CFD. Liang et al. (2018) does not utilize the forward pass simulation in the backward pass, limiting its utility for complex control synthesis tasks and system identification.
def step_ideal(q, v, h): # Ideal elastic collision simulation
a 0.0 & Contact-free acceleration
q, v jax. cond(q 0.0, lambda: (q. -V), lambda: (g, v)) Semi-inplicit Euler
v veve q a-q i-qiv q,
point mass
4o
target
def step_penaltyye h): Penalty collision simulation
Vo
0.0 Contact-free acceleration
QN
a_ref www.lax.condde. 0.0, lanbda: -1.8*v -50.0*q, lambda: 0.8; == a ref)
implicit suler
return 9.
Øjax wwlue_alddddaeaeeeeeee and_grad def unrollfq0, vê, n_steps, step_fn, h):
b=d
body(carry, _): carry 9. step_fnlg. V, h) return (9, v),
surface
(qf, vf), lax. sca. (g0, vB), lengthen_stops) e_target 1.0 loss [np q_target)
velocity change due contact
return
Figure 16: Jax code for minimal collision simulation. The code simulates a minimal version of penalty-based ideal-elastic simulator and is used in Fig. 4 to illustrate the TOI discretization emros as shown in Fig. 4.
18

---

Published as a conference paper at ICLR 2026
10
10a
10
10a
Loss RMSE
Grad RMSE
10-1
10-1
10
10-2_
10
10
10a
10
10
10a
Runtime Forward
Runtime Backward
10a
101
10
100
Loss RMSE
Grad RMSE
10
10
10-4
10
10
4 10
1x
10a
10'
6x10!
10²
JIT Time Forward (s)
JIT Time Backward (s)
MIX (EuleiSemilmplict
DirMMJX (Euler) DirmiX (Euler-Semllpplcct
[MJX* (Heun2)
DirMjJX (Dopri5)
(RK4) MJX (Impli(fFast)
mAjx (Bosh3)
DitMjx- (recommended)
Figure 17: Error of loss and gradient vs runtime and compilation time for different integrators. The loss and gradient are computed for the cube toss at an initial velocity of -2.0 with contact settings solceš=10_00 2.0] and selimp=[O.e 0.95 0.001 0.5 2] as depicted in Fig. 21 and the top of this figure. The ground-truth gradients are computed from finiteddifferenees using an adaptive solver with very low tolerance (10-12) to simulate the rollout. Standard fixed-stepsize MJX integrators have red colors, fixed-stepsize integrators in Diffrax have green color, and adaptive-stepsize integrators in Diffrax have blue color (also marked by *). For adaptive stepsize control we use the Difffrax PID controller with p 0.2, 0.4, D 0.0, as recommended for stiff ODEs. We observe that the pareto-front of adaptive solvers is shifted, allowing lower loss and gradient errors with less runtime. The best- performing solver is the Tsit5 solver (Tsitouras, 2011) 1), which is 5th order explicit Runge-Kutte Kutta method with an embedded 4th order method for adaptive step sizing. This is also the default solver used in other experiments.
19

---

Published as a conference paper at ICLR 2026
B TUNING DIFFMJX
Through Diffrax, DifiMIX provides a whole suite of additional tuning knobs that can be used to
improve gradient computation shown in Fig. 18.
DiffMJX
CFD
Soft collision detection
Artificial contacts gradient-basee optimization
MJX
Diffrax
Set contact stiffness
Set gradient accuracy
CFD- so mp d_ W. M power
Timestep
Error tol
solref
+
Checkpoints
solimp
Adaptive
CB Adjoint
integration
unroll backssotve
Figure 18: DiffMJX and CFD add tuning knobs to MJX's gradient computation. DiffMJX adds support for adaptive integration atop MJX. The integrator's 's tolerances trades computational speed for improved gradient accuracy. The number of checkpointing steps trades computation speed with GPU memory consumption. CFD adds parameters that the distance and magnitude which artificial contact forces between non-colliding objects are applied in the gradient computation.
Error tolerances set gradient accuracy. In this work, showed that the errors of penalty-based simulator gradients depends on
the simulation stiffness, aka the hardness of contacts and joint limits,
2. the contact velocities in normal direction, aka how fast objects collide into each other, 3. and the integration stepsize.
While simulation stiffness and contact velocities are system dependent, the integration stepsize is tuming parameter. To obtain accurate gradients, one sufficiently small stepsize which inevitably results in slow integration.
constant stepsize integrator with
To speed up computation, adaptive integrator solely reduces the stepsize if necessary function of the integrator's s error tolerances. In turn, as the integrator's ' error tolerances directly determines the gradient errors, it forms the most important hyperparameter.
To speed up computation, the largest error tolerance that still produces suffi- ciently accurate gradients.
Figure 17 shows the error vs runtime trade-off for different integrators in MJX and DiftMUX. We observe that in the case of the cube toss, the adaptive integrators reduce the runtime required to achieve gradients with low error.
Number of checkpoint trades GPU memory for compute time. Thanks Diffrax's advanced techniques for checkpointing, DiffMJX memory consumption is significantly reduced compared to MJX. As noted Diffrax's 's documentation, the memory used approximately equals the number of checkpoints multiplied by the size of the ODE's state.
Increasing the number of checkpoints speeds computation. Therefore, use as many checkpoints as possible given the available GPU memory.
By reducing the number of checkpoints, we were able run DiffMJX simulations that MJX could not, because MJX exceeded the available GPU memory.
Maximum number of ODE steps affects JIT compilation time. As shown in Fig. 17, there is a notable increase in JIT compilation time when using Diffrax compared to MuJoCo's semi-implicit Euler integrator. This is due to the program being compiled for the maximum number ODE steps that is set by the user. Therefore, this parameter should be chosen small as possible.
20

---

Published as a conference paper at ICLR 2026
C METHOD DETAILS
C.I REMOVING DISCONTINUITIS FROM COLLISION DETECTION AND MUSCLE ACTUATORS
The collision checking in MJX implemented for pairs of several different geometries. Many of the collision check implenentations rely on several case distinctions, e.g. the collision of a cylinder and a plane is separated into the case of the cylinder being parallel to the plane vs not. These hard case distinctions can lead to errors in the gradients, they introduce discontinuities in the dynamics and its gradient. We adopt the typical approach for smoothing such non-differentiabiliies by replacing any hard case distinctions by smoothly interpolating between the cases with weights computed from a sigmoid. Specifically, we end up modifying the collision detection between plane-cylinder pphere-capsule, capsule-capsule, plane-capsule and plane-box pairs. More complicated collisions like mesh-mesh collisions are currently not modified as they are not relevant in our experiments, but we plan on supporting these in a future version.
As an example, consider our soft version of a function for plane- box collision detection. More details on other smoothed functions will be available in the released codebase.
_pw.nc_boxttccccee ipiine GeonInta, **"Gaccllate contacte between pLane convex object.
ConvexIntor collision:
vert convex.vers
points in the planc_pos ccavex.mat. (planc.pcs convex.pos!
trame
convex.nst. 3 plane.nat 2] support (p) ane pos vert) 8
for manifold points within cpsilon skin
threshold tp.naxtmumr,, support.naxl 1e-3) sott poly nack math. stg.sss lowetneettold - 1e-3, highethresholdl support_necked math. wtt.:ooof_ccccoooott__ooo
11
12
support, 0)
13
14 dist_noe, idx 1ax. (ewp.ooopollppoott__sseee 4:
15 dist -dist_neg
16 poz - vert
17
18
convert
frame
19
pos convex.pos + por 8 convex.nat.I
20
plana.nat(t,
21
22
frame jp.stack (wwch..kaaaaa.o..oo... frane (n)] 4, axia=c)
23
pes # pes - 0.5 + dist Ranc]
24
return dist, pac, frame
shtt_ tcondition, VIa seturn concition
27
11 conditicni
We also encountered differentiation errors in MJX environments with muscle-tendon units. The source of these issues that the functions for computing wrapping of muscles around objects rely on nonlinearities such as arcsin, which has an undefined gradient at some points. A practical remedy is the "double-where" trick, which prevents gradient flow at those singularities.
Concurrently with our work, La Barbera et al. (2025) employed MJPC Howell et al. (2022) -a library for CPU-based MPC MuloCo via finite differences - to compute ideal controls of a musculoskeletal dog model. This seminal work also reports that modifying the muscle dynamics to be differentiable was a key contribution for improving MPC performance.
C.2 ODE DIFFERENTIATION IN DIFFRAX
Diffrra: Numerical integration with JAX The ODE solved by the Mujoco simulator can be
written
x0))= To
i(t) Fo(t, x(1),,
(5)
where 9 incorporates all model parameters and is the full state of the system.
In this work, we resort to the Diffrax library (Kidger, 2021; Kidger & Garcia, 2021) for numerical integration in Jax to solve the above ODE. This powerful library provides efficient implementations of many fixed and adaptive timestep solvers (Tsitouras, 201 Dormand & Prince, 1980). It also allows easily switching between different modes for backpropagation: (also referred to as unrolling) and optimize- thent-t-ttee-dieettiiiie (also referred to as backsolving) (Chen et al., 2018 8; Kidger, 2021). Overall, the choice of solver, stepsize controller and differentiation
21

---

Published as a conference paper at ICLR 2026
technique is typically application-depeee dependent. Hence, implement the general Diffrax integrator as an easy-to-use alternative for the existing fixed-stepsize integrators in MJX, offering the full flexibility of the Diffrax library to the user. Notably, we adjust the Diffrax solvers to accommodate for exact integration of quaternions and stateful actuators, similar to the Runge-Kutte implementation of MJX, which reduces the number of required integration steps.
The first notable differentiation mode is discretize- is the result of discretizing the forward ODE with a numerical integrator and then computing the gradients of the discretized forward ODE by unrolling the computation graph. This approach aligns with the paradigm of differentiable programming used in autodifferentiation libraries such as JAX, hence this is also the approach taken in standard MJX. The problem here is that the memory requirements scale linearly with the number of solver steps. For many applications this is a serious bottleneck in standard MJX, and the problem is exacerbated when relying on adaptive integration, as it typically involves storing even more substeps. Fortunately, Diffrax implements an optimal gradient checkpointing scheme that allows reducing memory requirements from ((n) to O(1) in exchange for increasing runtime from (n) to On log n.), where is the maximum number of solver steps.
Discretiz-tteenoptiiize
It
Another potential issue is that discretize- wwwreeiiee-teen-opiie then-optimize inherently computes gradients of the dis- cretized dynamics, rather than the true continuous dynamics. This can lead to interesting failure cases in which even in the limit of the stepsize going to zero, the computed gradient different from the true gradient of the continuous dynamics. This is the result of the basic fact that the convergence of a parameterized function (the discretized ODE) to a limit (the continuous ODE) does not imply the convergence of its gradient. Note that this is not mere mathematical corner case, instead it is the underlying reason for time-of-impact oscillations as described in Fig. 4 and Hu et al. (2020); Schwarke et al. (2024).
Fixing this requires (approximately) computing the gradients of the true continuous forward ODE. This is achieved by analytically computing the gradient as the solution to the continuous adjoint equations. Assuming loss I. [xx(T)) only on the final state w.l.o.g., the continuous adjoint equations (Kidger, 2021, Theorem 5.2) can be written as:
Optimize-then-disreee
dL
åz = (t) Oo
ax (T( dx(
(t, x(t))
(6)
Ox
oog
ag (T)
(t( = az (t)
(t, xt))
(7)
d0
Solving these using any numerical solver on the backward pass, as offered readily by the Diffrax dL
library, yields the desired gradients Im(e] (t) and di de (0).
This approach is called "optimize- cewemiii--tte--diicetttt the memory requirements can again be reduced to constant in the number of solver steps, by solving the forward ODE "backwards in time" together with the adjoint ODE. petimize-then-discreee contrast to unrolling, allows to directly specify tolerances on the error in the gradients.
BacksolveAdjoin in Diffrax. Note also that
in
One of the intuitive potential benefits of ceeiiiiz--tenndiseeeee for solving the adjoint ODE now has the ability to adapt the stepsize when computing the adjoint derivatives, rather than being constrained to unrolling the steps that were taken during the forward solve. However, in practice it is typically the case that the forward ODE and the adjoint ODE are stiff "in the same regions", i.e. the state derivative changes quickly when the state itself changes quickly. Hence, when unrolling the forward solve, the steps are already small in the regions where the adjoint ODE is also stiff. Therefore it typically suffices to use the less complicated discretize-then- -pttioie which has a less complex computational graph and therefore comes with lower compilation times. The above explanation of course only holds when the forward ODE and the adjoint ODE "match": when modifying the adjoint ODE (as will be discussed in the following section) it may be beneficial to use optimize- teeeeetthn-ddee discrette.
is that the adaptive stepsize controller
C.3 THEORETICAL RESULTS ON ADAPTIVE ODE INTEGRATION
There exists a rich history of results on analyzing discrete adjoint sensitivities, supporting the claims on gradient accuracy improvements from aon-adaetive solvers. For instance, Sandu (2006); Walther
22

---

Published as a conference paper at ICLR 2026
0,050
0.05
-6.05
0.04
0.64
Oraliirt
0.04
Gradlant
0.03
0.03
0.00
0.03
0.00
-0.05
0.05
0,050
0.01
Veliciity
in the presence of CFD. Results for reduced billiard example in Fig. without table-ball contacts and gravity. Left: DifIMJX adaptive integration results in well-behaved gradients, but gradient is zero if objects do not collide. Middle: Using CFD with discretize seneeiize-hee--oee oettinii results in small gradient oscillations. The adaptive integrator selects large stepsizes as no contacts are happening (forward ODE is non-stiff), whereas in the gradient computation the CFD are added (adjoint ODE is stiff). In turn, unrolling the large stepsizes cause discretization errors. Note that this phenomenon persists no matter how low the tolerance of the forward ODE solver set, as the stepsize controller is unaware of the CFD mechanism that is only used in the gradient computation. Right: Using optimize-then- thssericceteee allows the solver for the adjoint ODE to select small stepsizes when the adjoint ODE stiff due to CFD. This fixes the gradient oscillations.
Figure 19: Optinizz-then-disereiize vs
(2007) establish that for one-step Runge-Kutte Kutta methods of order-p. the reverse-mode sensitivities obtained by differentiating the discrete solve are themselves order-p consistent and stable. This directly links reductions in forward global error to reductions in gradient error the same order. Calvo et al. (2008) show that (under standard Lipschitz/stability conditions) the global state error scales like the (relative/absolute) tolerance enforced by the controller; the adjoint/sensitivity error inherits this scaling with constants depending on the assumed Lipschitz constant. As our setting is an application of adaptive integration to the ODE given by the simulator, these results apply out of the box and establish a theoretical basis for the observed improvements in gradient accuracy from the use of adaptive solvers. We hope that these changes will strengthen the presentation.
The advantages and disadvantages of discretize- theeoiizz-the---oiimizeee (DTO) vs optimize ceemiize-then--iine then-discretize (OTD) have also been researched for many years. This includes Hinze et al. (2008) comparing DTO and OTD, giving conditions for equivalence, and providing error bounds in the discretized optimality system. Hager (2001) proves accuracy of discrete optimality systems under Runge-Kutta time discretization and establishes conditions ensuring correct gradients. Moreover, Bonnans & Laurent- Varin (2006) provide order conditions for (partitioned) Runge-Kutt Kutta schemes in optimal control and implications for obtaining matching accuracy in the adjoint/gradient after discretization.
23

---

Published as a conference paper at ICLR 2026
C.4 CONTACTS FROM DISTANCE WITH THE STRAIGHT-THROURE TRICKK
In discretize-then-opiiiie automatic differentiation automaticaly uses the vector-jacobian-prodde derived from Fo instead of Fo. In cennnii--thenddicceetiie with the trrighh-throogh--rick applied to F, replace the jaa-jaaannannuuuuettt in the adjoint equations Eq. (6) as follows
dL
(t) (t) aoP 8x
(T) (T)
(t. x(t))
(8)
a (1, x()).
0g (T) =
(s( - -z (1)
do
Solving these altered adjoint equations gives replacements for the respective gradients, which in our case incorporate gradient information contacts from distance. Again, this approach does not require any additional implementation effort, as the neaight-through-trice allows to just use the existing Diffrax implementation of optimize-then--eee then-discretize.
-ptimize-then-diierettieee vs fiscretize-then-opiiiiz in the presence of CFD. As discussed in the previous section, adding contacts from distance may change the preferred choice of gradient computation, as the forward and adjoint ODE may have different stiffness values in different regions.
To illustrate this, we again consider version of the billiard-toy example from Fig. 8. This time, turn gravity off and disable the table-ball contacts to isolate the ball-ball contact, which means that the ODE linear before and after the collision. We again compute the loss and gradient over a range of initial parameters, the results are reported in Fig. 19, As before, we observe that with the DiffMJX adaptive integrator, we get well-behaved gradients, but the gradient becomes zero when the two balls do not collide.
Now, we activate the CFD mechanism and repeat the experiment. In this case, the gradient show oscillations due to the different stiffness settings of the forward and adjoint ODE: The solver for the forward ODE takes very large steps when no ball-ball collision happens, but the adjoint ODE incorporates the collision signal and hence is stiff. Simply unrolling the forward integration therefore makes a discretization error, resulting in the oscillations. Note, that we did not observe this phenomenon in Fig. 8, as this experiment includes ball-table collisions which cause the adaptive solver of the forward ODE to take small steps even when no ball- ball-ball collision is happening. The solution is to optimize- cirizzee-hen-deee then-discre Here, the adaptive solver for the adjoint ODE can select stepsizes according to the stiffness of the adjoint ODE. This is also confirmed by the results in Fig. 19 (right).
This experiments highlights a specific corner case, and we did not observe issues in our other experiments when using peccretiz--hennoptiiizz with CFD, we typically do have such extreme variations of the stiffness in the forward ODE. However, we believe it is beneficial to be aware of this potential caveat and how one can resolve it with the tools available in DiffMJX.
D DIFFMJX & CFD: THE SHARP BITS
CFD adds additional contacts to the solver which can slow down computation. The documen- tation of MuJoCo XLA advises to use meshes with 200 vertices or fewer due to computational limitations on the current implementation of collision detection. If MJX is used with CFD, then every contact is added to the contact solver for which r < We- In turn, for dense meshes and large We computation time notably increases. The parameter We can be made arbitrarily small such that, in its limit for w. 0, we retain vanilla MJX with our improvements the collision detector. That said, for in-hand manipulation tasks and robot locomotion tasks, the number of objects colliding at one instance is small enough to render CFD useful tool for easing optimization with penalty-based simulators.
the differentiability of
Adaptive integration slows down training if tolerances are set too small. For the very large initial cube side length initial value of 140 mm, DiffMUX saw a significant drop in computation speed forcing to abort these runs. This is not surprising penetrations amounting to 40% of the
"https://mu readthedocs. e/ee//eeeee/ee staale/mjjx html #mix- thhessharp----
24

---

Published as a conference paper at ICLR 2026
total cube's width cause huge contact forces that significantly stiffen the ODE. In Appendix E.2, we reduce the integration error tolerance which allowed optimization to converge even for large initial parameters, Alternatively, one could speed up training for these too large geometry parameter initializations by first using soft impedance settings that during training hecomes annealed to become increasingly stiff. The same principle applies to scenes with multiple objects: overly tight integration error tolerances can substantially slow down the simulation. As integration error scales with relative contact velocity, scenarios involving large state changes along contact normals trigger smaller step sizes. Notably, while the ball in the billiard shot example (Fig. 8) was rolling over the table, the integration stepsize did only marginally decrease. Whereas during the collision between balls, the integrator adjusted the stepsize notably.
Finite differences can be viable alternative for computing low-dimensional gradients. In general, gradient computation with MJX via automatic differentiation is notably slower than running only the forward simulation. This makes zeroth-order methods such as predictive sampling favorable whenever the task can be solved with them. This is usually the case for low-dimensional systems or computing gradients of only a few parameters, e.g. compting gradients via finite differences for less then ten parameters seems to work well in practice. As sampling suffers from the curse of dimension- ality, automatic differentiation for computing gradients scales to extremely high-dimensional control problems which enabled us to do gradient-based MPC on the musculoskeletal systems.
CFD should be combined with sampling to solve nonconvex optimizations. For MPC, our CFD mechanism tends to be effective when the object needs to be pushed in a specific direction by the hand, but it performs poorly in grasping tasks. This limitation arises because, at distance, the cumulative force arising from the hand's CFDs conglomerate into a force that pushes objects away. However, the gradient does not encode the possibility of grasping, as this signal only emerges when the object is inside the hand. Fundamentally, this issue reflects the non-convex nature of the optimization landscape. However, we do not believe that this is a problem that should be solved using CFD mechanism; rather, in such tasks, the best approach would be to combine CFD with sampling to overcome the non-convexity while maintaining the favorable scaling of gradients with dimensionality.
Parallelizing DiffMJX simulations. Locomotion policies trained with reinforcement learning typically rely on domain randomization with many environments executed in parallel. DiffMJX supports such batched execution via JAX's vmap, analogous MJX. However, a single policy induces large contact forces, then this may trigger very small integration steps, which results in the other simulations to wait for a slow simulation to finish. This slowdown is mitigated by increasing the integration error tolerances or by carly-terminating the respective simulation once predefined solver-step budget is exceeded. However, it important to note that such a slowdown is only observed when many contacts happening in each parallelized environment, there is no additional slowdown due to vmapping. This is because the selected adaptive stepsize in Diffrax is not global across the vmapped environments, it is fully vmapped and can take different values across environments. WE experimentally verify this in Fig. 20, see the caption for details.
happening in single environment. On the other hand, when few contacts
E EVALUATION DETAILS AND FURTHER EXPERIMENTS
E. GRADIENT ANALYSIS OF MJX
In Fig. 4, the time discretization errors arising in panelty- simulators are first illustrated minimal toy example. The code to generate this toy example is shown in Fig. 16. After illustrating on this toy example that time discretization errors can be mitigated by reducing the stepsize, it is shown in Fig. that these errors also occur in MuloCo. Fig. 21 is an extended version of Fig. 3 containing all of MuloCo's basic shape primitives.
25

---

Published as a conference paper at ICLR 2026
DiffMJX simulations of spheres colliding with ground increasing initial heights run simultaneously on GPU
100
Integrator stepsize (5)
Each simulation has Individual stepsize
10-4
Sim ID Solver steps
10-2
2/40
12:44
4/40
14/40
15/50
10a3
7/33
17/42
18/41
9/36
19/41
10/43
20/48
10a4
0.0
0.2
0.4
0.6
0.8
1.0
Simulation
Figure 20: We vmap the DifIMJX simulation of ball drop across 20 different environments. In each environment, a single ball is dropped (with an initial downwards velocity and no gravity) onto plane, from which it bounces back up. The 20 environments start with the ball different heights, such that the contact happens at a different simulation time for each environment. We then plot the selected adaptive stepsize of the solver over the simulation time for each environment, and observe that the stepsize indeed independent for each run, with a reduction in stepsize only occurring at the single contact in specific environment. Therefore, when vmapping across environments, contacts in one simulation do not slow down other simulations.
box
cylinder
capsule
ellipsold
sphere
Central differences
MIX Gradient
DiMMjx Gradient
Ix10*
-0.2
Gradient
-0.0
Loss/ TESS distance
-50
-0.75
Gradient
Gradienc scaled by 100
Gracient scaled 10c
Cradient scaled
Cradient 10000
DirMyx
dientttt.. 2
-0.4
-0.6
Groddont scaled 15
srodiantscaleddyy 3
Gradiest soard 700
Initial velocity V,
Initial velocity v.
Initial velacity
Initial velocity v
Initial velacity v.
Figure 21: Simulation of geometric primitives thrown onto surface using MJX or DiffMIX. Contacts at default stepsize 0.002 soires-ro saree-co.003 co.035 1.0), and colimp- (0.0 0.001 cause MJX's gradients of the toss distance to deviate significantly from central difference gradients, while DifIMJX maintains close agreement.
E.2 SYSTEM IDENTIFICATION EXPERIMENTS
Contactnet dataset The contactnet dataset² consists of 550 trajectories of an acrylic cube that has been repeatedly tossed onto a wooden table. As reported in Pfrommer et al. (2021); Acosta et al.
26

---

Published as a conference paper at ICLR 2026
Thaj ectory emor
Computation time
Sidelength
Mass
25.39
120
0.38
0.38
0.35
Trajecturyy
0.34
0
200
gradient step
gradient step
gradient stop
gradient stop
salref[0]
solref[
salimplo) (do)
solimp[1] (dw)
Trajectory ertrr
0.35
1,00
200
gradient
gradient
gradient
gradient
solimp[2] (width)
solimp133 (mid)
Triction(0]
Triction[1]
friction[2]
0.05
0.02
0.08
Tapestory aro
0.56
0.40
0.35
0.50
0.00
step
Trajectory error
Computation
Sidelength
Mass
0.45
31.91
180.7
0.35
0.26
0.10
0.0
0.28
gradient step
step
gradient step
grad step
Solre([0]
Solref[J]
solimp[O] (d0)
solimp[L] (dw)
0.09
0.99
0.40
0.0
0.1
100 gradient step
100 gradinnt stap
100 gradient stop
100
gradient step
solimp[2] (wisth)
solimp[3] (mid)
friction|0]
friction[1]
friction(2)
0.45
Feeeorr enor
0,30
0.00 100
0.0
0.00
0.00
100
200
100
100 gradiant stop
gradient stop
gradient stop
gradient stop
gradient stop
Figure 22: Top: Identification of MuJoCo parameters with DiffMJX on Contactnets cube toss dataset starting from pre-tuned inital conditions. Bottom: Identification starting from random initial parameters. The cube's geometry has significant effect on the time of impact such that the optimization converges to a sidelength close to 100 mm to reduce the train loss. The effect of other model parameters on the system dynamics is more inter-twined. For example, a too large mass resulting in larger penetrations during impact can be compensated by increasing the constraint stiffness via cairer(0) or colimpto).
(2022), data has been collected at 1480 Hz and the cube's physical parameters amount side length of 10 em, mass of 0.37 kg, inertia of 0.0081 kg m?, friction coefficient 0.18, and restitution of 0. 0..25.
27

---

Published as a conference paper at ICLR 2026
Training setup The trajectories split into segments of length five. In the supervised train loop, given a trajectory segment's 's initial state, each simulator is unrolled for four steps. Subsequently, an mean square error (MSE) loss between the segment's s states and predicted states is computed. Before being fed to the loss function, each state's quaternion is converted to a rotation matrix to avoid representation singularities negatively affecting training Geist et al. (2024); Brégier (2021). For adaptive integration, DifIMIX is set to use "RecusiveCheecpointii Runge-Kutta method ("Tsit5" in Diffrax) with PID error tolerances of le-5, We use the Adam Kingma & Ba (2015) implementation from the Optax library DeepMind et al. (2020) for gradient- based optimization with parameters [bi 0.5, 82 D.S, egs 1w-61. MuJoCo's parameter are set to timestep of 0.006767s (the data collection sampling time), steratione-t (contact solver iterations), sterstione-:0 (solver steps), to tollerance-ee-5 imprationl,0, rol 02 11, solimp-1C.00 0,95 0.0; a.s 21, and friction imp--ii-- 10.4 0-01 0.11. The cube's mass is set to the reported groundtruth value and its inertia is computed by MuJoCo's default equal density approximation. DiffMJX uses the same MuJoCo settings with the exception of iscrationn-2 (adaptive integration does not require as many solver iterations). The CFD impedance function set to.c, c.01, 0.01, 1.0, 4.01 that is de = 0.01, dg 0 with u 1 m.
and a 5th order explicit
E.2.1 IDENTIFICATION OF ADDITIONAL MuJoCo PARAMETERS
In this section, we extend the experiment from Section 5.1 and use DiffMJX to estimate MuJoCo's simulation parameters alongside the cube's side length. In these experiments, the same setup is used for DiffMJX as detailed in the previous section. In comparison to the experiment in Section 5.1, the integration error tolerance is increased to 1e-s to reduce computation time. Moreover, the parameters constrained using either a softplus function parameter has been carefully tuned. While the models are trained on an MSE between trajectory segments, we evaluate the model the "trajectory error" being the mean absolute error between the first forty trajectories in the dataset and DiffMJX's trajectory prediction. In the these experiments, DifIMJX estimates the following parameters simultaneously:
softclip function where the softening hyper-
sidelength also estimated in Section 5.1,
mass from which MuJoCo automatically computes the cube's inertia,
solref parameters (time constant tc, damping ratio Φd) determining the constraint stiffness,
solimp parameters determining the constraint's 's ability to generate force,
friction parameters determining the extend of the contact friction cone.
In total, we conduct two additional experiments with the Contactnets dataset:
Starting from pre-tuned initial conditions: The optimization is started from the same pre-tuned initial conditions as used in Section 5. The prediction horizon set to N = 10. The training results shown in Fig. 22 (Center). Each model requires around 10 seconds of computation time per gradient step.
Starting from random initial conditions: The optimization is started with random initial parame- ters. The prediction horizon is set to N = 4 such that the majority of runs require around seconds for a gradient step. The training results are shown in Fig. 22 (Bottom).
In both experiments, the trajectory error remains above 0.25. As noted by Parmar et al. (2021), the dynamics of hard cube exhibit degree of chaotic behavior. Consequently, small variations in initial conditions potentially arising from sensor noise or external disturbances such as wind can substantially alter the cube's 's trajectory after one of its corners hits the table. Nevertheless, the reduction in trajectory error from approximately 0.32 to 0.26 results in the simulated cube notably more closely matching the cal-worldd data, as shown in Fig. 22.
E.3 MPC EXPERIMENT
Experimental setup. The MyoHand and MyoArm collision parameters are set solsez=ço.02 1.01, and 2ollap=to.0 0.95 0.001 o.s 21, the racket collision parameters are set to se1xet=t-10c00 or for clastic collision. When using CFD, we set sot 100 0.95
0.001 0.5 and use eal imp-or 0.1 1.0 1.0 4).
28

---

Published as a conference paper at ICLR 2026
The timestep of the simulation is set to 0.0025, the constraint solver is the Newton solver with 4 iterations and 16 linesearch iterations. All models use our refinements for improving differentiability of MuJoCo's 's collision detector.
All MPC experiments are performed on NVIDIA GeForce RTX 3060. In the in-hand manipulation task the runtimes for the gradient-based MPC are 2.9h (+2.2h JIT) without CFD and 6.4h (+2.11h JIT) with CFD. The sampling with 1024 samples takes 35min (+-3iin JIT). In the bionic tennis task, the runtimes for the gradient-based MPC 9.7h (+1.2h JIT) with CFD. The sampling based MPC runs for 1.4h with 2048 samples.
Note, that did optimize any of the methods for speed in this experiment, i.e. in baoding with gradients the task was solved after less than half the total timesteps. The main point here is that we can solve the challenging task using a purely gradient-bases method that scales to high-dimensional systems. The large runtime of the gradient-based approach is mainly due to inefficient automatic differentiation of the MJX simulation, which remains the strongest limitation on gradient-based approaches at this point. One potential improvement could be to implement implicit differentiation of the constraint solver MJX, which could allow for much more efficient gradient computation.
29


Reviews:
[{"review_id": "lbgfMUSYzM", "reviewer": "Reviewer_SE5u", "summary": "The paper addresses issues that exist in using simulator gradients for robot learning. More specifically, it tackles two major issues: (a) the time quantization and discontinuities that happen at contact cause erroneous gradients; and (b) objects that are not in contact have no gradients (even when very close), thus making learning more difficult. For the first problem, this work presents an adaptive-timestep integration method that uses the disagreement between two numeric integrators of different orders to determine whether the step size needs to be reduced. The second issue of no gradients outside of contact is addressed by adding contacts from a distance, but this is avoided in the simulation, as it would yield unrealistic behaviour. Finally, they validate on some toy examples with a bouncing cube and some more realistic applications involving a tendon-driven arm and a ball (with MPC).", "questions": "- In Appendix D, the authors mention that one could reduce the tolerance, but then does this not simply result in the current setup of choosing a fixed step size that is a compromise between accuracy and time? Besides MPC and small-scale applications, where do you foresee the most uses in practice of the variable integration step size?\n\n- How hard is it to tune the gradient computation method for contacts from a distance? And how is this comparable to different loss formulations that exist, where these far-away contacts are added to the MPC/RL problem itself?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- The work addresses very relevant problems in simulators for robotics, and especially for using differentiable simulation to improve training.\n\n- There is a lot of engineering effort under the hood for both the simulator integration and the robot experiments.\n\n- The technical details in the paper are well presented and easy to follow.", "weaknesses": "- The variable step size on demand (i.e., during collisions) is a very poorly scaling option. It works well in the very specific cases of dropped objects or a single robot in contact with a ball that the authors present. However, in most applications of RL with multiple (1000s) parallel environments, the number of contacts would cause the simulator to continuously operate at a minimal step size, thereby bogging it down indefinitely.\n\n- The contacts from a distance approach introduced by the authors has been done under various forms as part of loss functions in RL (e.g., for dextrous manipulation [a]), where rewards are given based on distance and not contacts to avoid the discretization problem. The paper itself notes that the gradients are not used in the forward step but only in the backward pass to maintain realism. How these yet-to-be-done contacts are modelled on the gradient side is more of a task-level issue rather than the simulator providing gradients for a non-real physical effect.\n\n[a] Wang, Ruicheng, et al. \"Dexgraspnet: A large-scale robotic dexterous grasp dataset for general objects based on simulation.\" arXiv preprint arXiv:2210.02697 (2022).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PmTKGipaML", "reviewer": "Reviewer_S6dm", "summary": "This paper addresses one of the core challenges in differentiable physics: producing stable, accurate gradients through intermittent hard contacts such as impacts, frictional sticking, and restitution.\nThe authors first analyze why automatic differentiation in penalty-based engines (e.g., MuJoCo, Brax) produces pathological gradients near contact discontinuities. Together these techniques yield physically faithful forward dynamics and smooth, informative gradients for policy optimization and parameter identification. Experiments include gradient diagnostics on rigid-body impacts, inverse parameter identification, and manipulation policy learning, showing consistent improvements in convergence speed and stability.", "questions": "Please see the weaknesses sections. Thanks!", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper addresses a real and unsolved bottleneck—gradient instability under hard contacts.\n\n- The integration with MuJoCo-XLA makes the method directly applicable to a widely used robotics stack.\n\n- The numerical experiments (bouncing-ball, simple manipulator) clearly illustrate gradient discontinuities and the partial improvement achieved by DiffMJX.", "weaknesses": "**1.Relation to prior “force-from-distance” work.**\nThe proposed “Contacts From Distance (CFD)” estimator appears conceptually related to earlier “force-from-distance” smoothing used in quasi-dynamic contact models, such as Pang et al. (Global Planning for Contact-Rich Manipulation via Local Smoothing of Quasi-Dynamic Contact Models, arXiv:2206.10787). That paper also introduces continuous, distance-based contact forces that activate before penetration to improve numerical stability.\n\nHowever, it should be clarified whether the current work’s CFD mechanism differs in how gradients are propagated. The cited prior work modifies the forward dynamics to include soft forces, whereas this paper claims a backward-pass estimator that injects pre-contact gradients without changing the forward model. If that is indeed the novelty, it needs to be explicitly contrasted both conceptually and empirically.\n\n**2. CFD stability under multi-contact scenarios unclear.**\nCould the authors elaborate on whether CFD can handle overlapping distance fields or multiple contact pairs? What do you see as a new contribution here by CFD?\n\n**3. More realistic contact-rich scenarios will add value.** It will be great if the authors could showcase multi-fingered robotic in-hand reorientation, or a toy version of this task (see the 2D two-finger task in [1]). Could the authors discuss the advantages and limitations in applying CFD for trajectory optimizations of these tasks?\n\n[1] Pang et al., Global planning for contact-rich manipulation via local smoothing of quasi-dynamic contact models, 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "R0luHvZ93v", "reviewer": "Reviewer_7rxB", "summary": "This paper attempts to address efficient differentiable simulation under hard contacts, which is beneficial for accelerated learning and control. Traditional MuJoCo’s penalty-based contact resolution methods either (1) provide inaccurate gradients under hard contacts or (2) require excessively small $\\Delta t$ for more accurate gradient estimates, leading to unacceptably slow simulation runtime. Another issue with penalty-based contact resolution methods is uninformative gradients ($\\nabla = 0$), which does not provide any feedback in unrealized contacts (e.g. to steer away from / toward collision).\n\nThe paper proposes two core contributions: (1) combine adaptive time stepping and differentiable collision detection to provide more accurate gradients under hard contacts while keeping simulation efficiency; (2) CFD: simple modification of $d(r)$ and $h(r)$ for small non-zero artificial contact forces in unrealized contact scenarios. Implementation-wise, CFD is applied backward pass only via straight through tricks to maintain forward physics realism. The experiments are comprehensive and validate the main claims in the paper.", "questions": "1. How sensitive are results to CFD parameters ($w_c$, $d_c$, etc.)? Can you provide ablation studies or guidelines for setting these?\n2. Can you characterize more precisely when CFD helps vs hurts? The paper shows CFD works for directional pushing (billiard, tennis racket deflection) but acknowledges it fails for grasping. What about intermediate manipulation tasks such as pivoting or reorienting objects on a surface?\n3. You mentioned that you smooth the collision detection to make it differentiable in discrete case distinctions. Although some details were provided in the appendix, including some brief formal formulation or implementation details in the main text would be helpful for the readers.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The problem is interesting with potentially high, practical impacts on accelerated differentiable learning and control, as hard contacts are involved in almost all control problems.\n2. The proposal of adaptive time stepping is well-motivated and novel in the context of accelerated learning and control. Experiments sufficiently show that it helps improve gradient estimate accuracy under hard contacts while maintaining the simulator's computational efficiency. Although CFD is simple, it is well-motivated, has strong empirical performance, and is well-executed with straight through tricks. \n3. The experiments are comprehensive with diverse control tasks (controlled toy setups, on real cube-toss parameter identification, and gradient-based MPC on complex musculoskeletal systems). Comprehensive and through ablation studies on failure modes were included in the appendix.", "weaknesses": "1. Limited theoretical analysis of several claims. For example, while empirical results show that adaptive integration improves gradient accuracy, it lacks rigorous theoretical analysis of when and why this approach guarantees correct gradients. The authors can provide theoretical bounds on gradient error as a function of integration tolerance, or at minimum, a more rigorous analysis of when adaptive integration helps.\n2. Many explanations remain at the high level of intuition instead of formal theoretical analysis (e.g. the discussion of Optimize-then-Discretize benefits in gradient convergence). Thorough theoretical analysis of such cases would make the paper stronger.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)