Detailed Description

template<class T>
class drake::multibody::CenicIntegrator< T >

Convex Error-controlled Numerical Integration for Contact (CENIC) is a specialized error-controlled implicit integrator for contact-rich robotics simulations [Kurtz and Castro, 2025].

CENIC provides variable-step error-controlled integration for multibody systems with stiff contact interactions, while maintaining the high speeds characteristic of discrete-time solvers required for modern robotics workflows.

Benefits of CENIC include:

Guaranteed convergence. Unlike traditional implicit integrators that rely on non-convex Newton-Raphson solves, CENIC's convex formulation eliminates step rejections due to convergence failures.
Guaranteed accuracy. CENIC inherits the well-studied accuracy guarantees associated with error-controlled integration [Hairer and Wanner, 1996], avoiding discretization artifacts common in fixed-step discrete-time methods.
Automatic time step selection. Users specify a desired accuracy rather than a fixed time step, eliminating a common pain point in authoring multibody simulations.
Implicit treatment of external systems. This means that users can connect arbitrary stiff controllers (e.g., a custom LeafSystem) to the MultibodyPlant and have them treated implicitly in CENIC's convex formulation. This allows for larger time steps, leading to faster and more stable simulations.
Principled static/dynamic friction modeling. Unlike discrete solvers, CENIC can simulate frictional contact with different static and dynamic friction coefficients.
Speed. CENIC consistently outperforms general-purpose integrators by orders of magnitude on contact-rich problems. Error-controlled CENIC is often (but not always) faster than discrete-time simulation, depending on the simulation in question and the requested accuracy.

CENIC works by solving a convex Irrotational Contact Fields (ICF) optimization problem [Castro et al., 2023] to advance the system state at each time step. A simple half-stepping strategy provides a second-order error estimate for automatic step-size selection.

Because CENIC is specific to multibody systems, this integrator requires a system diagram with a MultibodyPlant subsystem named "plant".

Running CENIC in fixed-step mode (with error-control disabled) recovers the "Lagged" variant of discrete-time ICF simulation from [Castro et al., 2023].

References:

[Castro et al., 2023] Castro A., Han X., and Masterjohn J., 2023. Irrotational Contact Fields. https://arxiv.org/abs/2312.03908.

[Hairer and Wanner, 1996] Hairer E. and Wanner G., 1996. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Springer Series in Computational Mathematics, Vol. 14. Springer-Verlag, Berlin, 2nd edition.

[Kurtz and Castro, 2025] Kurtz V. and Castro A., 2025. CENIC: Convex Error-controlled Numerical Integration for Contact. https://arxiv.org/abs/2511.08771.

Template Parameters

T	The scalar type, which must be one of the default nonsymbolic scalars.

#include <drake/multibody/cenic/cenic_integrator.h>

Public Member Functions
	CenicIntegrator (const systems::System< T > &system, systems::Context< T > *context=nullptr)
	Constructs the integrator.
	~CenicIntegrator () final
const MultibodyPlant< T > &	plant () const
	Gets a reference to the MultibodyPlant used to formulate the convex optimization problem.
const contact_solvers::icf::IcfSolverParameters &	get_solver_parameters () const
	Gets the current convex solver tolerances and iteration limits.
void	SetSolverParameters (const contact_solvers::icf::IcfSolverParameters &parameters)
	Sets the convex solver tolerances and iteration limits.
bool	supports_error_estimation () const final
	Derived classes must override this function to indicate whether the integrator supports error estimation.
int	get_error_estimate_order () const final
	Derived classes must override this function to return the order of the asymptotic term in the integrator's error estimate.
Does not allow copy, move, or assignment
	CenicIntegrator (const CenicIntegrator &)=delete
CenicIntegrator &	operator= (const CenicIntegrator &)=delete
	CenicIntegrator (CenicIntegrator &&)=delete
CenicIntegrator &	operator= (CenicIntegrator &&)=delete
Public Member Functions inherited from IntegratorBase< T >
	IntegratorBase (const System< T > &system, Context< T > *context=nullptr)
	Maintains references to the system being integrated and the context used to specify the initial conditions for that system (if any).
virtual	~IntegratorBase ()
	IntegratorBase (const IntegratorBase &)=delete
IntegratorBase &	operator= (const IntegratorBase &)=delete
	IntegratorBase (IntegratorBase &&)=delete
IntegratorBase &	operator= (IntegratorBase &&)=delete
void	set_target_accuracy (double accuracy)
	Request that the integrator attempt to achieve a particular accuracy for the continuous portions of the simulation.
double	get_target_accuracy () const
	Gets the target accuracy.
double	get_accuracy_in_use () const
	Gets the accuracy in use by the integrator.
const ContinuousState< T > *	get_error_estimate () const
	Gets the error estimate (used only for integrators that support error estimation).
const T &	get_ideal_next_step_size () const
	Return the step size the integrator would like to take next, based primarily on the integrator's accuracy prediction.
void	set_fixed_step_mode (bool flag)
	Sets an integrator with error control to fixed step mode.
bool	get_fixed_step_mode () const
	Gets whether an integrator is running in fixed step mode.
const Eigen::VectorXd &	get_generalized_state_weight_vector () const
	Gets the weighting vector (equivalent to a diagonal matrix) applied to weighting both generalized coordinate and velocity state variable errors, as described in the group documentation.
Eigen::VectorBlock< Eigen::VectorXd >	get_mutable_generalized_state_weight_vector ()
	Gets a mutable weighting vector (equivalent to a diagonal matrix) applied to weighting both generalized coordinate and velocity state variable errors, as described in the group documentation.
const Eigen::VectorXd &	get_misc_state_weight_vector () const
	Gets the weighting vector (equivalent to a diagonal matrix) for weighting errors in miscellaneous continuous state variables z.
Eigen::VectorBlock< Eigen::VectorXd >	get_mutable_misc_state_weight_vector ()
	Gets a mutable weighting vector (equivalent to a diagonal matrix) for weighting errors in miscellaneous continuous state variables z.
void	request_initial_step_size_target (const T &step_size)
	Request that the first attempted integration step have a particular size.
const T &	get_initial_step_size_target () const
	Gets the target size of the first integration step.
void	set_maximum_step_size (const T &max_step_size)
	Sets the maximum step size that may be taken by this integrator.
const T &	get_maximum_step_size () const
	Gets the maximum step size that may be taken by this integrator.
double	get_stretch_factor () const
	Gets the stretch factor (> 1), which is multiplied by the maximum (typically user-designated) integration step size to obtain the amount that the integrator is able to stretch the maximum time step toward hitting an upcoming publish or update event in IntegrateNoFurtherThanTime().
void	set_requested_minimum_step_size (const T &min_step_size)
	Sets the requested minimum step size h_min that may be taken by this integrator.
const T &	get_requested_minimum_step_size () const
	Gets the requested minimum step size h_min for this integrator.
void	set_throw_on_minimum_step_size_violation (bool throws)
	Sets whether the integrator should throw a std::exception when the integrator's step size selection algorithm determines that it must take a step smaller than the minimum step size (for, e.g., purposes of error control).
bool	get_throw_on_minimum_step_size_violation () const
	Reports the current setting of the throw_on_minimum_step_size_violation flag.
T	get_working_minimum_step_size () const
	Gets the current value of the working minimum step size h_work(t) for this integrator, which may vary with the current time t as stored in the integrator's context.
void	Reset ()
	Resets the integrator to initial values, i.e., default construction values.
void	Initialize ()
	An integrator must be initialized before being used.
StepResult	IntegrateNoFurtherThanTime (const T &publish_time, const T &update_time, const T &boundary_time)
	(Internal use only) Integrates the system forward in time by a single step with step size subject to integration error tolerances (assuming that the integrator supports error estimation).
void	IntegrateWithMultipleStepsToTime (const T &t_final)
	Stepping function for integrators operating outside of Simulator that advances the continuous state exactly to t_final.
bool	IntegrateWithSingleFixedStepToTime (const T &t_target)
	Stepping function for integrators operating outside of Simulator that advances the continuous state using a single step to t_target.
void	ResetStatistics ()
	Forget accumulated statistics.
int64_t	get_num_substep_failures () const
	Gets the number of failed sub-steps (implying one or more step size reductions was required to permit solving the necessary nonlinear system of equations).
int64_t	get_num_step_shrinkages_from_substep_failures () const
	Gets the number of step size shrinkages due to sub-step failures (e.g., integrator convergence failures) since the last call to ResetStatistics() or Initialize().
int64_t	get_num_step_shrinkages_from_error_control () const
	Gets the number of step size shrinkages due to failure to meet targeted error tolerances, since the last call to ResetStatistics or Initialize().
int64_t	get_num_derivative_evaluations () const
	Returns the number of ODE function evaluations (calls to CalcTimeDerivatives()) since the last call to ResetStatistics() or Initialize().
const T &	get_actual_initial_step_size_taken () const
	The actual size of the successful first step.
const T &	get_smallest_adapted_step_size_taken () const
	The size of the smallest step taken as the result of a controlled integration step adjustment since the last Initialize() or ResetStatistics() call.
const T &	get_largest_step_size_taken () const
	The size of the largest step taken since the last Initialize() or ResetStatistics() call.
int64_t	get_num_steps_taken () const
	The number of integration steps taken since the last Initialize() or ResetStatistics() call.
std::vector< NamedStatistic >	GetStatisticsSummary () const
	Returns all integrator statistics as a single collection.
void	add_derivative_evaluations (double evals)
	Manually increments the statistic for the number of ODE evaluations.
const Context< T > &	get_context () const
	Returns a const reference to the internally-maintained Context holding the most recent state in the trajectory.
Context< T > *	get_mutable_context ()
	Returns a mutable pointer to the internally-maintained Context holding the most recent state in the trajectory.
void	reset_context (Context< T > *context)
	Replace the pointer to the internally-maintained Context with a different one.
void	reset_context (std::unique_ptr< Context< T > > context)
	Same as above but allows the integrator to take ownership of the context.
void	StartDenseIntegration ()
	Starts dense integration, allocating a new dense output for this integrator to use.
const trajectories::PiecewisePolynomial< T > *	get_dense_output () const
	Returns a const pointer to the integrator's current PiecewisePolynomial instance, holding a representation of the continuous state trajectory since the last StartDenseIntegration() call.
std::unique_ptr< trajectories::PiecewisePolynomial< T > >	StopDenseIntegration ()
	Stops dense integration, yielding ownership of the current dense output to the caller.
const System< T > &	get_system () const
	Gets a constant reference to the system that is being integrated (and was provided to the constructor of the integrator).
bool	is_initialized () const
	Indicates whether the integrator has been initialized.
const T &	get_previous_integration_step_size () const
	Gets the size of the last (previous) integration step.
std::unique_ptr< IntegratorBase< T > >	Clone () const
	Returns a copy of this integrator with reset statistics, reinitialized internal integrator states, and a cloned system context.

Static Public Attributes
static constexpr double	kDefaultAccuracy = 1e-3
	This target accuracy is established in the constructor, but may be changed by integrator-accuracy methods.

Additional Inherited Members
Public Types inherited from IntegratorBase< T >
enum	StepResult { kReachedPublishTime = 1 , kReachedZeroCrossing = 2 , kReachedUpdateTime = 3 , kTimeHasAdvanced = 4 , kReachedBoundaryTime = 5 , kReachedStepLimit = 6 }
	Status returned by IntegrateNoFurtherThanTime(). More...
const ContinuousState< T > &	EvalTimeDerivatives (const Context< T > &context)
	Evaluates the derivative function and updates call statistics.
template<typename U>
const ContinuousState< U > &	EvalTimeDerivatives (const System< U > &system, const Context< U > &context)
	Evaluates the derivative function (and updates call statistics).
void	set_accuracy_in_use (double accuracy)
	Sets the working ("in use") accuracy for this integrator.
bool	StepOnceErrorControlledAtMost (const T &h_max)
	Default code for advancing the continuous state of the system by a single step of `h_max` (or smaller, depending on error control).
std::pair< bool, T >	CalcAdjustedStepSize (const T &err, const T &attempted_step_size, bool *at_minimum_step_size) const
	Calculates adjusted integrator step sizes toward keeping state variables within error bounds on the next integration step.
virtual void	DoReset ()
	Derived classes can override this method to perform routines when Reset() is called.
virtual std::unique_ptr< IntegratorBase< T > >	DoClone () const
	Derived classes must implement this method to return a copy of themselves as an IntegratorBase instance.
trajectories::PiecewisePolynomial< T > *	get_mutable_dense_output ()
	Returns a mutable pointer to the internally-maintained PiecewisePolynomial instance, holding a representation of the continuous state trajectory since the last time StartDenseIntegration() was called.
bool	DoDenseStep (const T &h)
	Calls DoStep(h) while recording the resulting step in the dense output.
ContinuousState< T > *	get_mutable_error_estimate ()
	Gets an error estimate of the state variables recorded by the last call to StepOnceFixedSize().
void	set_actual_initial_step_size_taken (const T &h)
void	set_smallest_adapted_step_size_taken (const T &h)
	Sets the size of the smallest-step-taken statistic as the result of a controlled integration step adjustment.
void	set_largest_step_size_taken (const T &h)
void	set_ideal_next_step_size (const T &h)

Constructor & Destructor Documentation

◆ CenicIntegrator() [1/3]

template<class T>

CenicIntegrator ( const CenicIntegrator< T > & )

delete

◆ CenicIntegrator() [2/3]

template<class T>

CenicIntegrator ( CenicIntegrator< T > && )

delete

◆ CenicIntegrator() [3/3]

template<class T>

CenicIntegrator	(	const systems::System< T > &	system,
		systems::Context< T > *	context = nullptr )

explicit

Constructs the integrator.

Parameters

system	The overall system diagram to simulate. Must include a MultibodyPlant and associated SceneGraph, with the plant found as a direct child of the system diagram using the subsystem name "plant". The plant must be a continuous-time plant. This system is aliased by this object so must remain alive longer than the integrator.
context	context for the overall system.

◆ ~CenicIntegrator()

template<class T>

~CenicIntegrator ( )

final

Member Function Documentation

◆ get_error_estimate_order()

template<class T>

int get_error_estimate_order ( ) const

finalvirtual

Derived classes must override this function to return the order of the asymptotic term in the integrator's error estimate.

An error estimator approximates the truncation error in an integrator's solution. That truncation error e(.) is approximated by a Taylor Series expansion in the neighborhood around t:

e(t+h) ≈ e(t) + he(t) + he'(t) + ½h²e''(t) + ...
       ≈ e(t) + he(t) + he'(t) + ½h²e''(t) + O(h³)

where we have replaced the "..." with the asymptotic error of all terms truncated from the series.

Implementions should return the order of the asymptotic term in the Taylor Series expansion around the expression for the error. For an integrator that propagates a second-order solution and provides an estimate of the error using an embedded first-order method, this method should return "2", as can be seen in the derivation below, using y* as the true solution:

y̅ = y* + O(h³)   [second order solution]
ŷ = y* + O(h²)   [embedded first-order method]
e = (y̅ - ŷ) = O(h²)

If the integrator does not provide an error estimate, the derived class implementation should return 0.

Implements IntegratorBase< T >.

◆ get_solver_parameters()

template<class T>

const contact_solvers::icf::IcfSolverParameters & get_solver_parameters ( ) const

Gets the current convex solver tolerances and iteration limits.

◆ operator=() [1/2]

template<class T>

CenicIntegrator & operator= ( CenicIntegrator< T > && )

delete

◆ operator=() [2/2]

template<class T>

CenicIntegrator & operator= ( const CenicIntegrator< T > & )

delete

◆ plant()

template<class T>

const MultibodyPlant< T > & plant ( ) const

Gets a reference to the MultibodyPlant used to formulate the convex optimization problem.

◆ SetSolverParameters()

template<class T>

void SetSolverParameters ( const contact_solvers::icf::IcfSolverParameters & parameters )

Sets the convex solver tolerances and iteration limits.

◆ supports_error_estimation()

template<class T>

bool supports_error_estimation ( ) const

finalvirtual

Derived classes must override this function to indicate whether the integrator supports error estimation.

Without error estimation, the target accuracy setting (see accuracy settings) will be unused.

Implements IntegratorBase< T >.

Member Data Documentation

◆ kDefaultAccuracy

template<class T>

double kDefaultAccuracy = 1e-3

staticconstexpr

This target accuracy is established in the constructor, but may be changed by integrator-accuracy methods.

CENIC works best at loose accuracy.

The documentation for this class was generated from the following file:

drake/multibody/cenic/cenic_integrator.h

Detailed Description

Public Member Functions

Static Public Attributes

Additional Inherited Members

Constructor & Destructor Documentation

◆ CenicIntegrator() [1/3]

◆ CenicIntegrator() [2/3]

◆ CenicIntegrator() [3/3]

◆ ~CenicIntegrator()

Member Function Documentation

◆ get_error_estimate_order()

◆ get_solver_parameters()

◆ operator=() [1/2]

◆ operator=() [2/2]

◆ plant()

◆ SetSolverParameters()

◆ supports_error_estimation()

Member Data Documentation

◆ kDefaultAccuracy