CENIC: Convex Error-controlled Numerical Integration for Contact

Vince Kurtz; Alejandro Castro

cs.RONov 11, 2025

CENIC: Convex Error-controlled Numerical Integration for Contact

Authors: Vince Kurtz, Alejandro Castro

Abstract

State-of-the-art robotics simulators operate in discrete time. This requires users to choose a time step, which is both critical and challenging: large steps can produce non-physical artifacts, while small steps force the simulation to run slowly. Continuous-time error-controlled integration avoids such issues by automatically adjusting the time step to achieve a desired accuracy. But existing error-controlled integrators struggle with the stiff dynamics of contact, and cannot meet the speed and scalability requirements of modern robotics workflows. We introduce CENIC, a new variable-step integrator that brings together recent advances in convex time-stepping and error-controlled integration, inheriting benefits from both continuous integration and discrete time-stepping. CENIC runs at fast real-time rates comparable to discrete-time robotics simulators like MuJoCo, Drake and Isaac Sim, while also providing guarantees on accuracy and convergence.

Figures & tables

Panel of Figure 1 — Fig. 1 : A challenging bimanual manipulation scenario in simulation (top) and on hardware (bottom). A teleoperated robot drops several objects (spatula, spoon, mug) into a thin-walled bin. It then dumps the contents of the bin on a dish rack composed of thin wires. This scenario is particularly challenging: complex geometry, thin objects, heavy arms, and stiff joint controllers create a wide range of time scales that pose significant numerical difficulties. Discrete time-stepping methods struggle with passthrough, rattling, and unstable contact forces; CENIC eliminates these artifacts while completing the simulation at a real-time rate over 300%.

Fig. 2 : CENIC’s regularized friction model captures both static (

\mu_{s}

) and dynamic (

\mu_{d}

) friction regimes. Our default stiction tolerance is

v_{s}=0.1\,\text{mm/s}

, with a transition width of

\Delta=10

Fig. 3 : In practice, robotics simulations involve multibody dynamics connected to controllers, planners, and other components. We abstract these elements as a single external dynamical system. Unlike discrete formulations, CENIC integrates external systems implicitly, an effective strategy even for stiff feedback dynamics.

Fig. 4 : A simple nonlinear controller ( 30 ) induces a series of stable and unstable equilibria in a torque-controlled pendulum. CENIC’s linearization strategy integrates the controller implicitly in the stable (shaded, negative slope, negative feedback) regions. To maintain convexity, we fall back to an explicit scheme in the unstable regions. Enforcing convexity everywhere while performing implicit integration whenever possible enables faster and more reliable simulation of complex nonlinear controllers, see Sec. VII-D1 .

Fig. 5 : Vertical position and velocity of a ball dropped on a table, showing differences between solutions computed with a 1 ms (solid line) and a 5 ms (dashed line) time step. Velocity error spikes, but position error does not.

	Discrete Solvers [ 1 , 2 , 3 , 4 , 5 ]	Error-Controlled Implicit Integrators [ 15 ]	CENIC
Simulate Stiff Frictional Contact	Fast	Very Slow	Fast
Guaranteed Accuracy	No	Yes	Yes
Convex	Sometimes	No	Yes
Geometry Queries	One per step	One per iteration (many per step)	Two per step
Reduce Timestep	Never	For accuracy and convergence failure	For accuracy only
Implicit External Systems	No	Yes	Yes

TABLE I : Comparison of key characteristics of different simulation strategies. CENIC brings together the best of classical error-controlled integration and discrete solvers often used for robotics simulation.

Fig. 6 : Solver convergence during hard clutter simulation, illustrating the impact of Hessian reuse. Red dashed lines indicate a recomputed Hessian, while light grey lines separate time steps. ( 36 ) allows us to reuse the same Hessian when possible (first ten steps). When necessary, a few fresh factorizations enable rapid Newton-type convergence (final step).

Name	Degrees of Freedom	Max. Num. Contacts
cylinder	6	1
gripper	7	3
franka	16	2
soft clutter	120	31
hard clutter	120	77

TABLE II : Scale summary (degrees of freedom and maximum number of simultaneous contacts) for the test cases shown in Fig. 7 .

Panel of Figure 9 — Fig. 7 : Simulation test cases, ranging from least to most complex. These test cases highlight different challenges in robotics simulation: sharp corners ( 7(a) ), stable force closure ( 7(b) ), stiff nonlinear controllers ( 7(c) ), many degrees-of-freedom ( 7(d) ), and complex contact interactions ( 7(e) ).

Fig. 8 : Contact force magnitude for the static gripper example at accuracy

\varepsilon_{\text{acc}}=10^{-3}

. CENIC (top) produces smooth, stable contact forces, while RK3 (bottom) produces noisy, unstable results. Deceivingly, positions might look visually accurate with RK3, until later in the simulation when RK3 mistakenly drops the object.

Fig. 9 : Energy conservation error (percent of energy lost after 10 seconds) for a 0.1 kg bouncing ball with zero dissipation. The ball is relatively soft (contact stiffness

10^{3}

N/m) and bounces 11 times in the 10 second simulation. Potential energy is defined such that total energy is zero when the ball is at rest on the ground. The trapezoid method (Sec. V-B ) shows second-order convergence, while step-doubling (Sec. V-A ) is first order.

Fig. 10 : Work-precision plots for hard clutter (lower is better). Despite requiring only two convex solves per step versus three for step-doubling, the trapezoid method’s lack of stability hinders performance.

Fig. 11 : Work-precision plots for the test cases in Fig. 7 , using error control on the full state

\mathbf{x}

(top) and on positions

\mathbf{q}

only (bottom). Lower is better. Wall times are normalized to one second of simulation time. Missing data points indicate solver failure or timeout after 100 seconds (real-time rate

<1\%

). CENIC performs particularly well with many contacts and at loose accuracies, where convergence failures force other implicit schemes to take small steps.

Fig. 12 : Simulator speed comparison: wall times to complete a 10-second (top) and 3-second (bottom) simulation. Three mesh bunnies (top) and an unactuated quadruped robot (bottom) are dropped on the ground. CENIC is competitive with state-of-the-art robotics simulators in both fixed-step and error-controlled mode. Note that these benchmarks measure speed only, not simulation quality. Also note that the GPU simulators are optimized for simulating many parallel scenes, but these tests use only a single scene.

Fig. 13 : Number of time steps required to simulate the franka task for 10 seconds with different controller stiffnesses (

=K_{p}=K_{d}

) and accuracy

\varepsilon_{\text{acc}}=0.01

. With a standard explicit treatment, step size must decrease to keep stability under control. In contrast, our proposed implicit scheme (Sec. V-C ) is stable even with extremely stiff controllers.

Panel of Figure 16 — Fig. 14 : Suction gripper after 0.5s of simulation. For discrete time-stepping, we use $\delta t=0.01$ s. For CENIC, we request accuracy $\varepsilon_{\text{acc}}=10^{-3}$ . Large forces pulling the box into the gripper lead to massive passthrough artifacts in the discrete-time simulation. CENIC resolves the dynamics correctly with a combination of error control (Section V-A ) and implicit treatment of the gripper force model (Section V-C ).

Fig. 15 : Tangential velocities and friction forces for a 1 kg box on an oscillating conveyor belt, demonstrating transitions between static (

\mu_{s}=1.0

) and dynamic (

\mu_{d}=0.5

) friction. An accuracy of

10^{-5}

(left) shows the box “catching” in stiction at low velocities. A tight accuracy (

10^{-8}

, right) is needed to capture the impulsive transition from sliding to stiction. Existing convex time-stepping methods are unable to capture any of these effects, as they only support a single coefficient of friction.

Panel of Figure 18 — Fig. 16 : Undesirable artifacts from simulating the dish rack example (Fig. 1 ) with a too-large time step (fixed-step mode) or a too-loose accuracy (error control mode). Tightening the accuracy is a more efficient means of eliminating these artifacts than reducing the time step, as shown in Table III

Accuracy	$10^{-1}$	$10^{-2}$	$\mathbf{10^{-3}}$	$\mathbf{10^{-4}}$
Real-time rate	218%	463%	342%	132%
Artifacts	Yes ( 16(d) )	Yes ( 16(e) )	No	No
Time Step (s)	$0.02$	$0.01$	$0.002$	$\mathbf{0.001}$
Real-time rate	102%	180%	43%	21%
Artifacts	Yes ( 16(a) )	Yes ( 16(b) )	Yes ( 16(c) )	No

TABLE III : Performance comparison between error-controlled integration (top) and fixed-step integration (bottom) on the dishrack simulation (Fig. 1 ). Error-controlled integration eliminates obvious artifacts (shown in Fig. 16 ) while achieving faster overall simulation speeds.

Panel of Figure 20 — Fig. 17 : Snapshots of the dish rack task on hardware (top) and in simulation with CENIC (bottom, $\varepsilon_{\text{acc}}=10^{-3}$ ). The simulation and real trajectories are not exactly the same, due to model sensitivity and chaotic dynamics. But the overall qualitative behavior in simulation matches the real-world demonstration: there are no drops, missed grasps, passthrough incidents, rattling, or other major numerical artifacts. In contrast, standard discrete time-stepping methods introduce major artifacts, even with relatively small time steps (Fig. 16(a) - 16(c) ).

Fig. 18 : Real-time rate and step size

\delta t

for the complex dishrack simulation (Fig. 1 ) at accuracy

\varepsilon_{\text{acc}}=10^{-3}

. By only taking small time steps when needed to resolve stiff collision dynamics, CENIC achieves faster overall simulation times than discrete time-stepping (Table III ).

Fig. 19 : Computational cost breakdown for each of the examples, run with error control at accuracy

\varepsilon_{\text{acc}}=10^{-3}

Accuracy	$10^{-1}$	$10^{-2}$	$10^{-3}$	$10^{-4}$	$10^{-5}$
Wall time (adaptive)	0.19	0.39	0.65	1.11	2.69
Wall time (fixed)	0.21	0.41	0.71	1.19	2.69
Iterations (adaptive)	4015	9973	19235	39537	95839
Iterations (fixed)	4591	10938	20865	40729	95839

TABLE IV : Impact of adaptive convergence criteria (Sec VI-B ) on hard clutter simulation times (in seconds) and total solver iterations.

Error control	Reuse	Wall time (s)	Steps	Facts.	Iters.
Off	Off	0.26	1000	3980	3980
Off	On	0.27	1000	2293	6068
On	Off	0.65	499	16624	16624
On	On	0.61	510	10557	19235

TABLE V : Impact of Hessian reuse on the hard clutter example. Hessian reuse reduces the number of expensive factorizations (“Facts.”), at the cost of more iterations. With error control enabled (

\varepsilon_{\text{acc}}=10^{-3}

, maximum timestep 0.1 s), this offers a small improvement in overall performance. In fixed step mode (

\delta t=10

ms), Hessian reuse does not improve performance.

Fig. 20 : Comparing different strategies for choosing an initial guess for linesearch, for the hard clutter problem at different accuracies. Our cubic spline method consistently reduces the number of iterations required for exact linesearch, resulting in modest improvements in overall speed.

Contact Stiffness (N/m)	Steps	Iterations	Wall time (s)
$10^{2}$	585	2232	0.22
$10^{4}$	440	11179	0.34
$10^{6}$	664	24765	0.92
$10^{8}$	1378	37183	1.742
$10^{10}$	1621	34647	1.53
$10^{12}$	1641	37495	1.63

TABLE VI : Statistics for hard clutter simulation (10 seconds of simulation time) at accuracy

\varepsilon_{\text{acc}}=10^{-3}

, and varying contact stiffness

k_{i}

. Stiffer contact parameters lead to more time steps, more iterations, and longer wall times. Nonetheless, CENIC is able to simulate this contact-rich scene significantly faster than real-time even with extremely stiff contact parameters (note that steel has stiffness around

k_{i}=10^{7}

N/m [ 17 ] ).

Robotics Simulation

Robotics simulation builds virtual robots and environments to model, train, and evaluate autonomous systems before real-world deployment.

See all 320 papers

GPU-Accelerated Continuous-Time Successive Convexification for Contact-Implicit Legged Locomotion
Apr 11, 2026Samuel C. Buckner, Purnanand ElangoTrajectory Optimization Robot Locomotion
Department of Aeronautics and Astronautics, University of Washington, Seattle, WA 98195. · Mitsubishi Electric Research Laboratories, Cambridge, MA 02139.
Highly-Efficient Differentiable Simulation for Robotics
Sep 11, 2024Quentin Le Lidec, Louis Montaut, Yann de Mont-Marin +3Robotics Simulation Multibody Dynamics
Inria, Ecole normale sup´erieure CNRS, PSL Research University 75005 Paris, France
FlashNeRD: Performance-First Contact-Rich Neural Robot Dynamics
Oct 6, 2026Mohammadmehdi Ataei, Mohammad Amin Nabian, Jie Xu +9Robot Dynamics Learning Robotics Simulation
NVIDIA

Apr 11, 2026J/K move · Enter open · S save

GPU-Accelerated Continuous-Time Successive Convexification for Contact-Implicit Legged Locomotion

Samuel C. Buckner, Purnanand Elango

Department of Aeronautics and Astronautics, University of Washington, Seattle, WA 98195. · Mitsubishi Electric Research Laboratories, Cambridge, MA 02139.

Trajectory Optimization Robot Locomotion Legged Robotics Nonconvex Optimization

Open paper AskPDF arXiv

Contact-implicit trajectory optimization (CITO) enables the automatic discovery of contact sequences, but most methods rely on fine time discretization to capture all contact events accurately, which increases problem size and runtime while tying solution quality to grid resolution. We extend the recently proposed sequential convex programming (SCP) approach for trajectory optimization, continuous-time successive convexification (ct-SCvx), to CITO by introducing integral cross-complementarity constraints, which eliminate the risk of missing contact events between discretization nodes while preserving the flexibility of contact mode changes. The resulting framework, contact-implicit successive convexification (ci-SCvx), models full multibody dynamics in maximal coordinates, including stick-slip friction and partially elastic impacts. To handle complementarity constraints, we embed a backtracking homotopy scheme within SCP for reliable convergence. We implement this framework in a stand-alone Python software, leveraging JAX for GPU acceleration and a custom canonical-form parser for the convex subproblems of SCP to avoid the overhead of general-purpose modeling tools such as CVXPY. We demonstrate ci-SCvx on diverse legged-locomotion tasks. In particular, we validate the approach in MuJoCo with the Gymnasium HalfCheetah model against the MuJoCo MPC baseline, showing that a tracking simulation with the optimized torque profiles from ci-SCvx produces physically consistent trajectories with lesser energy consumption. We also show that the resulting software achieves faster solve times than existing state-of-the-art SCP implementations by over an order of magnitude, thereby demonstrating a practically important contribution to scalable real-time trajectory optimization.

Abstract

Robotics Simulation

Explore similar work

GPU-Accelerated Continuous-Time Successive Convexification for Contact-Implicit Legged Locomotion

Highly-Efficient Differentiable Simulation for Robotics

FlashNeRD: Performance-First Contact-Rich Neural Robot Dynamics