Source-linked AI summary
Non-Smooth Newton Methods for Deformable Multi-Body Dynamics
Miles Macklin, Kenny Erleben, Matthias Müller, Nuttapong Chentanez, Stefan Jeschke, Viktor Makoviychuk
TL;DR
Robotics simulation for dexterous manipulation must couple richer physical models, while LCP-based approaches can struggle with nonlinear dynamics and poorly conditioned problems. This paper solves nonlinear complementarity problems with a Newton-based framework, and evaluates it on robotics scenarios where it performs well on deformable-object manipulation and ill-conditioned rigid contact.
Problem
Dexterous manipulation and soft robots require robust simulation of multiple physical models, while hyperelastic forces and nonlinear contact challenge LCP formulations.
Source: §1 · p. 1Method
A Newton-based framework solves multi-body NCPs combining nonlinear contact, articulated rigid-body dynamics, and hyperelastic materials, with symmetric linear systems as its solver building block.
Source: §1 · p. 2Results
Across robotics scenarios, the method performs well on deformable-object grasping and dexterous manipulation, as well as ill-conditioned rigid-body contact problems.
Source: §12 · p. 17Takeaways & Limitations
The authors conclude that their complementarity preconditioner makes non-smooth contact formulations practical for interactive applications.
Source: §12 · p. 17Krylov solvers’ global line search lets residual error in one row affect reduction in others; choosing the complementarity preconditioner’s r factor also leaves room for better heuristics.
Source: §11 · p. 17
Abstract
from arXiv · showhide
We present a framework for the simulation of rigid and deformable bodies in the presence of contact and friction. Our method is based on a non-smooth Newton iteration that solves the underlying nonlinear complementarity problems (NCPs) directly. This approach allows us to support nonlinear dynamics models, including hyperelastic deformable bodies and articulated rigid mechanisms, coupled through a smooth isotropic friction model. The fixed-point nature of our method means it requires only the solution of a symmetric linear system as a building block. We propose a new complementarity preconditioner for NCP functions that improves convergence, and we develop an efficient GPU-based solver based on the conjugate residual (CR) method that is suitable for interactive simulations. We show how to improve robustness using a new geometric stiffness approximation and evaluate our method's performance on a number of robotics simulation scenarios, including dexterous manipulation and training using reinforcement learning.
1 INTRODUCTION
Robotics simulation can make learning safer and less costly, but dexterous manipulation needs robust coupling of multiple physical models. The paper addresses this with a Newton-based NCP framework for nonlinear dynamics and contact.
1 INTRODUCTION
- Simulation offers a safer, less costly alternative to real-world data gathering for learning robot behaviors, but extending transfer beyond locomotion requires richer multi-physics simulation.The motivating applications include dexterous manipulation and soft robots.Source: §1 · pp. 1–2
- Hyperelastic internal forces and nonlinear contact challenge LCP-based methods, whose relaxation solvers can converge slowly and direct solvers scale poorly.The paper identifies combining relaxation methods’ simplicity and robustness with direct methods’ accuracy as an open challenge.Source: §1 · p. 2
- The framework solves multi-body NCPs with Newton’s method, combining nonlinear contact, articulated rigid-body dynamics, and hyperelastic materials as semi-smooth equations.Off-the-shelf linear solvers serve as the computational building block.Source: §1 · p. 2
- The contributions include smooth isotropic Coulomb friction as complementarity functions, a convergence-improving preconditioner, generalized compliance for hyperelastic materials, and an approximate geometric stiffness for robustness.Source: §1 · p. 2
2 RELATED WORK
Prior work addressed contact, friction, and rigid–deformable coupling with complementarity, relaxation, and Newton-based methods. This paper distinguishes its approach through combined contact-force updates, smooth friction-cone treatment, and a generalized compliance formulation.
2.1 Contact
- Earlier frictional contact methods used linearized cones and LCP or staggered QP solves; this work instead solves contact and friction together using nonlinear complementarity functions and symmetric linear systems.Jean and Moreau [1992] introduced implicit frictional contact time-stepping, while Stewart and Trinkle [1996] used a linearized cone and Lemke iterations.Source: §2.1 · p. 2
- Relaxation methods can converge slowly on poorly conditioned problems, while Gauss-Seidel is order-dependent and difficult to parallelize; the authors solve all contacts simultaneously.Their formulation works directly with friction limit surfaces, which they say offers modeling freedom for anisotropic or non-symmetric cones.Source: §2.1 · p. 2
- Newton methods can yield quadratic convergence at higher per-iteration cost, and this work extends the approach to deformable bodies and nonlinear friction models.The authors contrast their approach with relaxation methods and report handling more ill-conditioned problems when combined with their preconditioner.Source: §2.1 · pp. 2–3
2.2 Coupled Systems
- For rigid–deformable coupling, the authors put elastic and contact dynamics in one system update, unlike prior multi-stage updates; a diagonal geometric-stiffness approximation avoids dense decompositions and supports iterative solvers.The approximation is described as improving stability.Source: §2.2 · p. 3
- The paper proposes a recipe for constructing compliance forms from arbitrary material strain-energy density functions and demonstrates it for stable Neo-Hookean material.The authors contrast complementarity contact constraints with prior work using stiff penalty forces.Source: §2.2 · p. 3
3 BACKGROUND
The background formulates dynamics with equality constraints, compliant elasticity, unilateral contact, and friction. These pieces combine into a differential variational inequality with coupled normal and frictional complementarity conditions.
3 BACKGROUND
- The background begins with a second-order differential equation for system dynamics, with generalized coordinates and their time derivatives describing configuration and motion.Source: §3 · p. 3
3.1 Hard Kinematic Constraints
- A system’s second-order dynamics use its mass matrix, generalized force, and configuration coordinates; bilateral constraints add Lagrange-multiplier forces and yield a differential-algebraic equation.The exposition assumes scleronomic bilateral constraints for simplicity.Source: §3 · p. 3, §3.1 · p. 4
3.2 Soft Kinematic Constraints
- Compliance, the inverse stiffness coefficient, expresses elastic constraints in a form mathematically equivalent to quadratic energy potentials and numerically stable as stiffness tends to infinity.The compliance matrix is diagonal in the stated formulation.Source: §3.2 · p. 4
3.3 Contact
- Contact is an inelastic non-interpenetration inequality along a fixed world-space normal; creating constraints within a fixed feature distance works at small time-steps but can over-constrain configurations.The contact points may be rigid-body frame points or deformable-body particles; more sophisticated constraints are cited as an alternative [Williams et al. 2017].Source: §3.3 · pp. 4–5
3.4 Friction
- Friction follows maximal dissipation, with a two-dimensional tangential basis and force magnitude bounded by the normal force through the Coulomb constraint.The model’s complementarity condition describes sticking and slipping.Source: §3.4 · p. 5
3.5 Governing Equations
- Combining the motion equations with normal and frictional complementarity conditions gives a differential variational inequality; coupling those conditions makes the problem non-convex.The passage states that implicit time-integration leads to an NP-hard optimization problem [Kaufman et al. 2008].Source: §3.5 · pp. 5–6
4 NONLINEAR COMPLEMENTARITY
The section reformulates contact and friction conditions as nonlinear complementarity functions for root-finding. It uses minimum-map and Fischer-Burmeister formulations for contact and a fixed-point formulation for smooth friction.
4 NONLINEAR COMPLEMENTARITY
- A non-smooth NCP function converts complementarity into a root-finding problem when combined with time-discretization.The function's roots satisfy the original complementarity conditions.Source: §4 · p. 6
4.1 Minimum-Map Formulation
- Minimum-map reformulates unilateral contact as an active-set method: contacts with Cn(q) ≤ λn are active, while inactive contacts have zero multiplier.At the friction cone limit, the minimum-map and Fischer-Burmeister friction error functions are both zero, indicating sliding.Source: §4.1 · p. 6, Figure 6 · p. 6
4.2 Fischer-Burmeister Formulation
- The Fischer-Burmeister function is smooth except at (a,b) = (0,0), where its derivative is set-valued and Newton methods may choose a subgradient.The authors choose a derivative at that point for simplicity.Source: §4.2 · p. 6
4.3 Frictional Constraints
- The friction fixed-point iteration uses either minimum-map or Fischer-Burmeister functions and has fixed points satisfying the original complementarity condition.Its compliance term W projects frictional force back onto the friction cone, and treating W as constant yields a symmetric system.Source: §4.3 · p. 7
5 IMPLICIT TIME-STEPPING
The section discretizes the equations of motion implicitly in generalized velocity and coordinates, with contact and friction represented through constraint functions. The resulting equations are nonlinear and require no additional stabilization terms or projection passes.
5 IMPLICIT TIME-STEPPING
- Backward time-discretization updates generalized coordinates from generalized velocity, and the implicit equations solve for end-of-step velocities and multipliers.The unconstrained velocity includes external and gyroscopic forces integrated explicitly; the discretization makes friction impulsive.Source: §5 · p. 8
- The implicit discrete equations remain nonlinear but require no additional stabilization terms or projection passes.Unlike index reduction, this formulation does not add stabilization terms to combat drift or project the system back to the constraint manifold.Source: §5 · p. 8
- The friction cone remains smooth rather than being faceted, avoiding extra per-facet multipliers and directional approximation bias.The discrete model treats contact forces as impulsive, preventing sliding and interpenetration over a single time-step.Source: §5 · p. 8
6 NON-SMOOTH NEWTON
The method applies generalized Newton iterations to the discretized equations, using generalized Jacobians for non-smooth functions and approximate linear solves and derivatives. System assembly leads to a saddle-point problem that is reduced with a Schur complement.
6 NON-SMOOTH NEWTON
- Newton's method applies to the non-smooth equations using a matrix from the generalized Jacobian, with approximate solves and derivatives making the method inexact and quasi-Newton.The generalized Jacobian is the convex hull of directional derivatives at a non-smooth point.Source: §6 · p. 8
6.1 System Assembly
- The assembled system includes geometric stiffness and contact-function derivatives, while friction compliance is evaluated each Newton iteration using current multipliers.The frictional forces are bounded by the normal force from the previous Newton iteration.Source: §6.1 · p. 9
6.2 Schur-Complement
- Linearizing the discrete equations gives a saddle-point system that is indefinite and possibly singular; a mass-block Schur complement yields a reduced positive-semidefinite system.After solving for the multiplier change, the velocity change is evaluated directly.Source: §6.2 · p. 9
- The alternative compliance-block Schur complement requires an invertible compliance matrix, which can be poorly conditioned for stiff materials.The paper uses the mass-block compliance formulation to accommodate stiff constraints with deformable bodies.Source: §6.2 · pp. 9–10
7 COMPLEMENTARITY PRECONDTIONING
The section introduces complementarity preconditioning by scaling the arguments according to the system and discrete time-stepping equations. The scaling addresses differing physical units and aligns the slopes of the complementarity terms.
7 COMPLEMENTARITY PRECONDTIONING
- The preconditioner chooses a scaling factor from system properties and discrete time-stepping equations to improve convergence.The motivation is that complementarity arguments can have different physical units, which can impair convergence.Source: §7 · p. 10
- Scaling uses effective system mass and time-stepping equations to put both complementarity arguments in the same units.For unilateral constraints, the factor is the time-step-scaled effective mass and makes the terms have the same slope.Source: §7 · p. 10
- For friction constraints, the scaling factor is h times the corresponding diagonal entry, and a Jacobi preconditioner makes this scaling inexpensive.The diagonal of the system matrix is already computed when using Jacobi preconditioning.Source: §7 · p. 10
8 ROBUSTNESS
The section improves robustness through solver choices and Newton-system approximations. It also describes practical choices for line search, initialization, geometric stiffness, and regularization.
8.1 Line Search and Starting Iterate
- For frictional problems, merit-function backtracking often stalls, so the reported results use a damped Newton step instead.The method accepts a constant fraction of the full Newton step, which may temporarily increase the merit function.Source: §8.1 · pp. 10–11
- Zero velocity is the most robust starting iterate, including under the large random torques used in reinforcement learning.Warm-starting constraint forces can improve efficiency, but it was omitted from reported results because of bookkeeping.Source: §8.1 · p. 11
8.2 Preconditioned Conjugate Residual
- PCR often yields an order-of-magnitude lower residual at equal iterations and decreases residual monotonically, supporting useful early termination.The authors also report convergence on close-to-singular contact systems where CG and Cholesky may fail.Source: §8.2 · p. 11, Figure 11 · p. 11
8.3 Geometric Stiffness
- A diagonal geometric-stiffness approximation improves Newton robustness without changing the converged solution; omitting it causes overshooting and oscillations at large strains.The approximation updates the mass block diagonally and avoids explicitly evaluating constraint Hessians.Source: Figure 12 · p. 12, §8.3 · p. 12
8.4 Regularization
- An optional ϵ identity shift addresses singular systems from redundant constraints, regularizing Newton-iteration error without changing the physical model.The examples use ϵ = 10−6, equivalent to adding a small amount of compliance.Source: §8.4 · p. 12
9 MODELING
The modeling section extends compliance-based elasticity to general material models and forces. It describes continuum material choices alongside particle and rigid-body representations.
9.1 Generalized Compliance
- Generalized compliance factors forces to move part of the force Jacobian out of the mass matrix and into a compliance regularization term.The approach isolates poorly conditioned components, such as large or infinite material parameters.Source: §9.1 · pp. 12–13
9.2 Continuum Materials
- Continuum materials are specified through a strain-energy density Ψ(s), where s can represent strains, strain rates, or principal stretches.Each finite element contributes constraint and multiplier variables and a block to the system compliance matrix.Source: §9.2 · pp. 13–14
9.2.1 Linear Materials.
- For linear Hookean constant-strain elements, the compliance matrix E = (VeKe)−1 is constant and can be precomputed.Here Ke is the element material stiffness matrix and Ve its rest volume.Source: §9.2.1 · p. 14
9.2.2 Hyperelastic Materials.
- Linear materials may lose significant volume under large deformation, while hyperelastic models are less prone to this artifact.For isotropic hyperelasticity, stiffness depends on configuration and is evaluated and inverted at each Newton iteration; E may require projection to remain positive definite.Source: Figure 13 · p. 13, §9.2.2 · p. 14
9.3 Particles
- Deformable bodies use tetrahedral finite elements over particles, with three degrees of freedom per particle.Source: §9.3 · p. 14
9.4 Rigid Bodies
- Rigid bodies use maximal coordinates for center-of-mass position, quaternion orientation, and generalized velocity.The quaternion is periodically normalized to unit length rather than constrained directly in the Newton solve.Source: §9.4 · p. 14
10 RESULTS
The experiments cover rigid and deformable contact, reinforcement-learning tasks, solver behavior, and numerical accuracy. They report stable grasps and structures, improved convergence with complementarity preconditioning, and an accuracy model that generally predicts observed limits.
10.1.1 Fetch Tomato.
- Stable grasps were formed while transferring a tomato, grasping DexNet objects, and lifting a ball with the deformable PneuNet gripper.The tomato task couples frictional contact with internal dynamics; the PneuNet finger uses an activation model rather than an explicit chamber cavity.Source: §10.1.1 · p. 15, §10.1.3 · p. 15, §10.1.4 · p. 15
10.1.2 Fetch Beam.
- Relaxation methods could not achieve the beam’s desired joint stiffness even after hundreds of iterations.The beam is modeled as 16 rigid bodies connected by joints with 250 N.m bending stiffness.Source: §10.1.2 · p. 15
10.1.5 Rigid Body Contact.
- Krylov methods form stable structures and prevent interpenetration in difficult rigid-contact scenes, while relaxation methods struggle with high mass ratios.The tests include a self-supporting arch, a heavy box stack, and table piling.Source: §10.1.5 · p. 15, Figure 14 · p. 15, §10.4 · p. 17
10.1.6 Material Extension.
- Large-strain simulation was stable with geometric stiffness, while the linear co-rotational model showed high volume loss relative to the hyperelastic model.The comparison uses Young’s modulus E = 105Pa and Poisson’s ratio ν = 0.45.Source: §10.1.6 · p. 15
10.1.7 Reinforcement Learning.
- The simulator supported learned policies for opening a drawer and for humanoid running, with the latter using stick-slip transitions during fast turns.The drawer-opening policy was learned in fewer than 100 training iterations.Source: §10.1.7 · p. 16, Figure 15 · p. 16
10.2 Effect of Complementarity Preconditiong
- Across 132 beam-insertion steps, the proposed complementarity preconditioner produced super-linear convergence in more cases and often an order-of-magnitude lower median error than identity scaling.For rigid contact, combining time-step and effective-mass scaling gave the most reliable observed convergence.Source: §10.2 · p. 16, Figure 16 · p. 16
10.3 Effect of NCP-Function
- Fischer-Burmeister produced smoother contact-force distributions than the minimum-map for redundant contact sets, and was used for the examples because of its improved interactive stability.The authors note that the proposed explanation based on the functions’ non-smooth points needs further investigation.Source: §10.3 · p. 16
10.5 Error Analysis
- With 32-bit floating point, the error model generally predicted observed accuracy limits, except that it underestimated residual accuracy for the flexible beam.The linear-solver comparison fixes the iteration count at 40; early termination can cause erratic convergence for PCG, while relaxation and PCR are generally unaffected.Source: §10.4 · pp. 16–17, §10.5 · p. 17
11 LIMITATIONS AND FUTURE WORK
The formulation is intended to work with different generalized-coordinate parameterizations, but several solver and modeling limits remain. The authors identify possible improvements in linear solves and complementarity scaling.
11 LIMITATIONS AND FUTURE WORK
- The formulation can use different generalized-coordinate parameterizations by supplying the corresponding mass matrix, kinematic map, and constraint Jacobians.This supports reduced-coordinate articulations, which encode hard constraints directly and reduce parameters and explicit constraint equations.Source: §11 · p. 17
- Krylov solvers’ global step can let error in one row slow residual reduction in other independent simulation parts.The authors suggest detecting islands and solving independent systems separately.Source: §11 · p. 17
- Elastic collisions and energy-preserving integrators are not specifically addressed, though the authors believe Smith et al. [2012] can be combined with their formulation.Source: §11 · p. 17
- The geometric-stiffness approximation is less effective on rigid-articulated mechanisms and can cause jitter at joint limits.The authors therefore apply it only to particle degrees of freedom and identify symmetric rank-1 updates as future work.Source: §11 · p. 17
- A symmetric linear solve is the method’s primary computational cost; matrix-free iterative solvers and algebraic multigrid are proposed as ways to improve performance.The authors note that algebraic multigrid may accelerate larger-scale linear subproblems.Source: §11 · p. 17
- The complementarity preconditioner’s r-factor has a large design space, and global-information heuristics may improve performance.Source: §11 · p. 17
12 CONCLUSION
The paper presents a framework for rigid and deformable contact problems that can use off-the-shelf linear solvers. Across robotics scenarios, it performs well on deformable-object manipulation and ill-conditioned rigid contact problems.
12 CONCLUSION
- The framework supports off-the-shelf linear solvers for rigid and deformable contact problems and performs well on deformable-object manipulation and ill-conditioned rigid contact scenarios.The authors believe their complementarity preconditioner makes non-smooth contact formulations practical for interactive applications.Source: §12 · p. 17
A CONTACT FORMULATION
The appendix supplies contact and friction NCP forms and derivatives, then gives the compliance-form derivatives for stable Neo-Hookean materials.
A CONTACT FORMULATION
- The appendix gives full contact constraints and derivatives for the minimum-map and Fischer-Burmeister NCP functions, including r-factor scaling.Source: A CONTACT FORMULATION · p. 20
A.1 Contact
- The Fischer-Burmeister contact constraint and its derivatives are provided in closed form.Source: §A.1 · p. 20
- The minimum-map contact formulation and its derivatives are also given.Source: §A.1 · p. 20
A.2 Friction
- For frictional conditions, the appendix gives the full Fischer-Burmeister NCP function, including r-factor scaling.Source: §A.2 · p. 20
- The Fischer-Burmeister friction equation yields a fixed-point iteration for |DT Ûq| and |λf|.Source: §A.2 · p. 20
- The frictional compliance quantities are given for both NCP functions, and used to form the compliance sub-block for one contact.Source: §A.2 · p. 20
B COMPLIANCE FORM OF STABLE NEO-HOOKEAN MATERIALS
- For stable Neo-Hookean materials, the appendix gives the elastic strain-energy density and the Jacobian and Hessian with respect to principal stretches.Source: B COMPLIANCE FORM OF STABLE… · p. 20
- For a constant-strain tetrahedron, elastic potential energy is rest volume times strain-energy density; compliance is obtained by a 3 × 3 inverse, using singular-value derivatives from Perez et al. 2013.Source: B COMPLIANCE FORM OF STABLE… · p. 21