cs.ROOct 26, 2026

Differentiable Contact Dynamics for Legged Locomotion: Policy Learning and Contact-Implicit Trajectory Optimization

Authors: Gijeong Kim, Seungbin Oh, Dongyun Kang, Hae-Won Park

Organizations: The authors are with the Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Daejeon 34141, South Korea.

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

This paper proposes a smooth analytic gradient for frictional contact and applies it to policy learning and contact-implicit trajectory optimization for legged locomotion. The exact contact sensitivity is available in closed form but describes changes within the current contact mode, providing no first-order information toward neighboring contact modes. The proposed gradient replaces this sensitivity only in the backward pass, leaving the hard-contact forward solve unchanged and the realized rollout feasible. A central-path relaxation of the de Saxcé-corrected second-order-cone complementarity problem for Coulomb friction smooths the conditions defining the separating, clamping, and sliding modes. The relaxed sensitivity is evaluated at a point on the central path obtained by a Newton iteration initialized from the hard-contact solution. The sensitivity remains continuous through contact-mode transitions and converges to the exact sensitivity as the relaxation vanishes. Using 100,000 samples and less than 2 min of training on a desktop PC, policy learning with the proposed gradient produces a policy that transfers to hardware without fine-tuning. On hardware, the robot under a trotting policy withstands pushes and traverses debris and slippery ground. In trajectory optimization, the proposed gradient reaches contact-rich solutions missed by exact and Signorini-relaxed gradients. Project page: https://diffcontact.github.io

Figures & tables

Fig. 1 : Converged iLQR solutions for three CITO tasks over a 30-step horizon with hard-contact forward dynamics. The solutions differ only in the backward gradient: mode-dependent exact, Signorini-relaxed, or proposed full-cone relaxation. (a) Lift-off: the mass starts clamped on the ground with the target above its initial position. (b) Traverse: the mass starts clamped with a lateral target. (c) Friction braking: without gravity, the mass hovers above the ground at the target position with an initial lateral velocity and returns after moving away. All rollouts start from zero control. An MIQP over the contact-mode sequence, solved to zero optimality gap, provides the reference trajectory (dashed) and cost normalization. Right: final cost normalized by the MIQP cost and contact-mode sequence over the horizon; the green marker denotes the MIQP solution.
Fig. 2 : Central-path relaxation of the exact contact solutions in the clamping, separating, and sliding regimes for a two-dimensional contact with μ=0.8\mu=0.8 . The top row shows the contact impulse λ∈Kμ\bm{\lambda}\in\mathcal{K}_{\mu} , and the bottom row shows the de Saxcé-corrected contact velocity σ+Γμ(σ)∈Kμ∗\bm{\sigma}+\Gamma_{\mu}(\bm{\sigma})\in\mathcal{K}^{*}_{\mu} . At ρ=0\rho=0 , the clamping regime has σh=0\bm{\sigma}_{h}=\bm{0} , the separating regime has λh=0\bm{\lambda}_{h}=\bm{0} , and, in the sliding regime, both λh\bm{\lambda}_{h} and σh+Γμ(σh)\bm{\sigma}_{h}+\Gamma_{\mu}(\bm{\sigma}_{h}) lie on their respective cone boundaries. For ρ>0\rho>0 , the corresponding relaxed variables lie in int⁡Kμ\operatorname{int}\mathcal{K}_{\mu} and int⁡Kμ∗\operatorname{int}\mathcal{K}^{*}_{\mu} . Black markers denote the exact solution. Blue and orange markers denote the relaxed impulse λ⋆\bm{\lambda}^{\star} and corrected contact velocity σ⋆+Γμ(σ⋆)\bm{\sigma}^{\star}+\Gamma_{\mu}(\bm{\sigma}^{\star}) , respectively, for ρ∈(0,0.5]\rho\in(0,0.5] . Darker markers correspond to larger ρ\rho , and dotted curves trace the relaxed-solution loci as ρ\rho increases.
Fig. 3 : Sensitivities of a single contact under the exact, randomized-smoothing, and proposed gradients: the impulse and contact-velocity components with respect to the free normal velocity bNb_{N} . The point-mass model has diagonal AA and nonzero free tangential velocity bTb_{T} . As bNb_{N} increases, the hard-contact solution changes mode from clamping to sliding and then to separating. The exact sensitivity (black, ρ=0\rho=0 ) jumps at the mode boundaries and contains zero components resulting from differentiation of the contact constraints with the mode held fixed. Randomized smoothing with perturbation standard deviation 0.60.6 yields zeroth- and first-order estimates [ 16 , 30 ] using N=128N=128 perturbed solves per point. The proposed sensitivity ( ρ=0.3\rho=0.3 ) is evaluated deterministically at a single Newton-solved central-path point for each bNb_{N} (Section V-B ) and remains continuous across these mode transitions.
Fig. 4 : Effect of the normal free-velocity offset δ\delta in the proposed relaxation. (a) Increasing δ\delta translates the dynamics relation (dashed). Its intersection with the relaxed contact law λNσN=ρ\lambda_{N}\sigma_{N}=\rho moves along the curve toward lower λN\lambda_{N} and higher σN\sigma_{N} . (b) Normal sensitivity ∂σN/∂bN\partial\sigma_{N}/\partial b_{N} as a function of the unshifted bNb_{N} . Increasing δ\delta shifts the sensitivity curve to more negative bNb_{N} .
Fig. 5 : Overview of the proposed differentiable contact framework. Given (xt,ut)(x_{t},u_{t}) , the hard-contact forward solve determines the contact impulse λ\lambda , contact velocity σ\sigma , and next state xt+1x_{t+1} under the realized sliding, clamping, or separating mode. For the backward pass, the Signorini condition, Coulomb friction cone, and maximum-dissipation principle are expressed as a single conic complementarity condition using the de Saxcé correction. The hard-contact solution at ρ=0\rho=0 is relaxed along the central path to ρ>0\rho>0 , and analytic differentiation yields the relaxed contact sensitivity and transition Jacobians (fx,fu)(f_{x},f_{u}) . These Jacobians are used for BPTT policy learning and DDP-based trajectory optimization.
TermWeightTermWeight
Base Joints 
Position, x,yx,y00Position5×10−15\times 10^{-1}
Height4×1014\times 10^{1}Hip-roll position11
Orientation, roll–pitch2×1012\times 10^{1}Velocity5×10−35\times 10^{-3}
Orientation, yaw00Feet 
Linear velocity, x,yx,y44Position, xx33
TABLE I : Running-cost weights for policy learning
Fig. 6 : Training with three backward contact gradients (no relaxation, Signorini-relaxed, and proposed) on four locomotion tasks (10 seeds). Training uses 16 environments, with friction coefficient μ∈[0.2,0.4]\mu\in[0.2,0.4] for low-friction gait-free locomotion and μ∈[0.6,1.0]\mu\in[0.6,1.0] otherwise. (a) Training reward. Curves show the mean and shaded bands denote ±1\pm 1 standard deviation across training curves smoothed separately for each seed. (b) Violin plots of evaluation cost. Each violin shows evaluation costs from 2,250 in-distribution rollouts across commands, friction coefficients, and seeds. Each cost is averaged over pre-fall steps. White markers and thick bars denote medians and interquartile ranges, respectively. Median costs are listed below the violins, with the lowest per task in bold. (c) Evaluation fall rate over an extended friction range; gray shading denotes out-of-distribution friction.
Fig. 7 : Hardware deployment of BPTT-trained locomotion policies. Snapshots show the hardware experiments, and the plots show estimated foot-contact patterns. (a) Trot, (b) pace, and (c) gait-free locomotion. Training uses 100,000 control steps in 108 s for the trot and pace policies, and 160,000 steps in 172 s for the gait-free policy.
Fig. 8 : Hardware robustness tests of trotting policies trained with 100,000 samples in under 2 min. Snapshots show representative trials for each experiment. (a) Traversal over scattered debris. (b) Recovery from external pushes. The plot shows lateral-velocity tracking error over 14 trials; gray curves show individual trials, the green curve shows the median, and the shaded interval marks the push. (c) Locomotion on a low-friction surface. Base roll and pitch and estimated foot-contact modes are shown during pronounced foot slippage; the highlighted interval marks the slippage event.
Fig. 9 : Training-reward comparison between BPTT and PPO on the trot command-tracking task. (a) BPTT and (b) PPO training reward over policy iterations for different numbers of environments; curves show the mean and shaded bands denote one standard deviation. (c) Training reward versus cumulative environment steps for BPTT with one environment and PPO with 64 and 128 environments. Among the tested configurations, BPTT learns the trot with one environment, whereas PPO first learns it with 128 environments; PPO with 64 environments does not learn the trot.
Fig. 10 : Comparison of BPTT and PPO in trained-policy evaluation and training wall-clock time versus the number of environments (mean across seeds; shading denotes one standard deviation). Policy evaluation uses 1010 -s trot rollouts over sampled velocity commands and in-distribution friction values. BPTT is shown through 64 environments and PPO through 1024. (a) Linear-velocity tracking RMSE. (b) Angular-velocity tracking RMSE. (c) Evaluation cost. (d) Training wall-clock time on a desktop (log scale); dashed lines mark 2 and 20 minutes.
Fig. 11 : Newton convergence and root selection in the backward contact evaluation during BPTT trotting training. (a) Cumulative percentage of all contact evaluations reaching the residual tolerance ∥Fρ∥<ε\|F_{\rho}\|<\varepsilon within kk Newton iterations of Algorithm 1 . (b) Percentage of all contact evaluations that satisfy the residual tolerance but fail the cone-feasibility check.
Fig. 12 : Comparison of the three backward contact gradients for CITO base-pose reaching. (a) Translation : the reference distance ranges from 0.20.2 to 1.61.6 m. (b) Yaw rotation : the reference yaw ranges from 0.170.17 to 3.143.14 rad. (c) Pitch rotation : the reference pitch ranges from 0.40.4 to 1.21.2 rad. For each target, the finite-horizon objective is averaged over 400 successive solves along a 10-s trajectory. Lines and shaded regions show the mean and one standard deviation over eight translation directions in (a), ten planar offsets within 0.1 m in (b), and ten planar offsets within 0.3 m in (c). The top row shows the initial standing and reference poses as solid and transparent robots, respectively.
Fig. 13 : Comparison of a single CITO problem over DDP iterations under the three contact gradients. All optimizations use the same current state, base-pose reference, and initial state and control trajectories for the 0.750.75 -m rightward Translation case in Fig. 12 (a); the initialization is taken at t=1.30t=1.30 s of the exact-gradient rollout. (a) Cost versus DDP iteration. (b)–(d) Initial guess and optimized trajectories after 25 iterations over the 0.50.5 -s horizon: (b) base roll angle, (c) RR-foot contact mode, and (d) RR-foot lateral displacement.
Fig. 14 : Comparison of the three contact gradients in CI-MPC under unmodeled external pulls. (a) Maximum pull magnitude over eight nominal directions. Each value is the largest pull magnitude for which at least 90%90\% of rollouts with perturbed pull directions complete without failure. (b)–(d) Responses to a 250250 -N forward pull: (b) base tilt, (c) planned front-right (FR) foot contact mode, and (d) base and FR-foot xx -positions. Crosses in (d) mark failures due to non-foot body contact with the ground.
Appendix figures & tables1 asset

Supplementary material from the paper’s appendix.

Appendix

optimizer / learning rateAdamP / 5×10−45\times 10^{-4}
discount / GAE λ\lambda0.99 / 0.95
clip ϵ\epsilon0.2
epochs / minibatches per update10 / 2
value / entropy coefficient0.5 / 0.01
gradient-norm clip0.5
action std (init / bounds)1.0 / [0.2, 1.5][0.2,\,1.5]
Appendix
TABLE II : PPO settings.
More in this topic

Differentiable Physics

Differentiable physics makes physical simulators usable for gradient-based learning, parameter estimation, design, and control.

See all 66 papers

Explore similar work

CardsList
Sep 25, 2026

Bundled Contact Gradients: Stabilizing Differentiable Simulation for Deployable Dynamic Tasks

Differentiable simulation provides analytic gradients of robot dynamics, enabling fast and sample-efficient first-order policy optimization. However, obtaining smooth and informative gradients through rigid-body contact typically requires softened contact models, often at the expense of physical fidelity and thereby limiting learned policies largely to simulation. This trade-off becomes particularly consequential for dynamic humanoid motions, where accurate contact dynamics are critical for transferring policies to the real world. Increasing contact stiffness in rigid-body simulation improves the fidelity of interactions, but also makes the dynamics increasingly sensitive to small state perturbations, producing high-variance gradients that can destabilize first-order policy learning. To address this, we propose Bundled Contact Gradients (BCG), a contact-local randomized smoothing framework for differentiable policy learning. When stiff contact is detected, our method evaluates a local bundle of randomized perturbation rollouts around the stiff contact configuration and aggregates their gradient signal thereby reducing gradient variance. We demonstrate the effectiveness of our method by successfully training and transferring dynamic motions zero-shot onto a real-world Unitree G1 humanoid platform. Videos and supplementary information can be found at https://bundledcontactgradients.github.io/

Dyuman Aditya, Jin Cheng, Clemens Schwarke +4

University of Southern California, LA, USA. · ETH Zurich, Switzerland. · NVIDIA. +1

Robotics SimulationDifferentiable PhysicsRandomized SmoothingPolicy Gradient MethodsRobot Locomotion
AskPDF
Oct 14, 2026

LEAP: Gradient Flow Augmented Lagrangian Trajectory Optimization for Contact Constrained Robots

Generating dynamically feasible motions for legged and manipulating robots requires solving a large, nonlinear, non-convex optimal control problem (OCP) in which contact wrenches must satisfy friction-cone constraints throughout the motion, not only at isolated instants. Existing solvers are sensitive to the initial guess and scale poorly with the number of joints. This paper introduces LEAP, which transcribes the OCP with pseudospectral collocation and, instead of a Newton or SQP solve, integrates a primal--dual gradient flow of a Powell--Hestenes--Rockafellar augmented Lagrangian to compute a solution. Because the flow is stiff, this paper designs a linearly implicit integrator with a Gauss--Newton metric, momentum, and a merit-based step controller that keeps the augmented Lagrangian monotone. Across locomotion benchmarks on five robot platforms, from planar hoppers to 32-joint humanoids, LEAP is the only method to solve every tested problem, achieves substantial speedups over Crocoddyl and Aligator on several benchmarks, and reaches trackable references at smaller compute budgets. The resulting trajectories are validated by tracking them in a high-fidelity simulator and on hardware. We will make the code public upon final acceptance of this paper.

César E. Ramos Chuquiure, Ryan Tsai, Yutong Zhou +8

Department of Robotics, University of Michigan, Ann Arbor, MI, USA. · Computer Science and Engineering, University of Michigan, Ann Arbor, MI, USA. · Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI, USA.

Hypergradient OptimizationTrajectory OptimizationRobotic OptimizationAugmented Lagrangian MethodRobot Locomotion
AskPDF
Sep 8, 2026

Ostrich: Taking Large Strides Through Stiff Contact in Differentiable Dynamics

Three properties determine whether a differentiable simulator can drive gradient-based optimization through contact: simulation accuracy, gradient reliability, and per-iteration cost. Tape-based engines such as MJX and Newton Semi-Implicit require timesteps small enough to keep contacts numerically tractable, and their backpropagation memory grows linearly with the number of timesteps T. Surrogate models bound memory by approximating contact away, but the resulting gradients lose the geometry the optimization depends on. We present Ostrich, a GPU-accelerated rigid-body simulator that resolves hard contacts and friction with non-smooth Newton iteration at large timesteps (h ~ 0.1 s), and differentiates the converged residual via the implicit function theorem, reusing the forward Schur complement to compute the adjoint at O(1) memory per timestep. On real-robot trajectories over a pallet obstacle, Ostrich holds MuJoCo's sim-to-real accuracy up to a 50x larger timestep. Its gradients converge from random initializations where MJX descends slowly and Newton Semi-Implicit stalls; a warm iteration runs 211x faster than MJX's and 4.7x faster than Semi-Implicit's. On the same scene Ostrich differentiates 8,192 parallel worlds on a single 24 GB GPU, sustaining 29x checkpointed MJX's optimization throughput; without checkpointing both baselines exhaust memory at far fewer worlds. We close with a gradient-based trajectory optimization demonstration over triangle-mesh terrain across a 10 s horizon, a setting where prior engines either restrict to primitive geometry or face the convergence and memory limits shown above.

Aleš Kučera, Karel Zimmermann

Department of Cybernetics, Faculty of Electrical Engineering, Czech Technical University in Prague, Prague, Czech Republic

Trajectory OptimizationContact-Rich ManipulationMultibody DynamicsDifferentiable ProgrammingDifferentiable PhysicsPhysics-Based Simulation
AskPDF