Contact Models and Multibody Dynamics

Contact Models

Solving contact dynamics is NP-hard problem[1] due to its non-convexity and discontinuity. To make the problem tractable, physics engines approximate the contact problem as a relaxed problem that enables using efficient solving methods.

One common method is polytope approximation that relaxes Coulomb's friction cone constraint as pyramid shape. This approximation transforms Non-linear Complementarity Problem (NCP) to Linear Complementarity Problem, which has the efficient solving methods such as Lemke[2], Dantzig[3] and Projected Gauss-Seidel. This approach is mathematically clear but sacrifices reality of simulation significantly, especially when the simulation has many slipping contacts.

ODE and DART apply this approach for contact simulation. They provide LCP Dantzig and LCP Projected Gauss-Seidel (PGS) solver that can be selected by a user. While Dantzig solver finds an exact solution of LCP, PGS solver finds an approximated solution iteratively. Note that DART's Dantzig solver is cloned from ODE Dantzig solver.

Bullet approximates the problem into Mixed Linear Complementary Problem and uses Sequential Impulse solver to find the solution that also approximates the friction cone into polytope. However it post-processes the solution to satisfy round shape friction cone constraint.[4]

RaiSim uses the Bisection method[5] that iteratively finds impulse solution which satisfies both Coulomb's friction cone constraint and zero normal velocity constraint. It uses a geometrical interpretation of the NCP problem to find the feasible solution in a remarkably efficient manner.

While the above engines hold strict complementarity in velocity level, MuJoCo is based on the complementarity-free method that approximately enforces the complementarity condition.[6] It relaxes complementary constraint and makes the hard contact problem into convex. With this relaxation, the penetration between object happens, but the solution can be obtained efficiently.

Multibody Dynamics

Solving dynamics of a robotic articulated system with $n$ links has $O(n^3)$ complexity which is very expensive. Most state-of-the-art physics engines, thus, adopted the Featherstone algorithm that solves the dynamics for open kinematic chain structure in $O(n)$ complexity.[7] It uses minimal (or reduced) coordinate formulation that represents $n$ link robot's configuration with $n$ or $n+6$ variables for fixed based and floating based robot respectively.

RaiSim, Bullet, MuJoCo, and DART uses minimal coordinate representation and the Featherstone algorithm for the multibody dynamics. On the other hand, ODE does not use minimal coordinate, instead use the maximal coordinate representation with $6n$ variables.


  1. D. Kaufman et al., “Staggered projections for frictional contact in multibody systems,” 2008.
  2. R. W. Cottle, J. Pang and R. E. Stone, "The Linear Complementarity Problem," 1992.
  3. E. Coumans, "Exploring MLCP solvers and Featherstone," GDC (2014)
  4. J. Hwangbo, J. Lee, and M. Hutter, "Per-Contact Iteration Method for Solving Contact Dynamics," IEEE Robotics and Automation Letters (2018).
  5. E. Todorov, "Convex and analytically-invertible dynamics with contacts and constraints: Theory and implementation in MuJoCo," ICRA (2014).
  6. R. Featherstone, "Rigid Body Dynamics Algorithms," Springer Verlag (2008)