Non-Smooth Newton Methods for Deformable Multi-Body Dynamics

TL;DR

This paper introduces a non-smooth Newton method for simulating contact and friction in deformable multi-body systems, supporting complex nonlinear dynamics with improved convergence and GPU acceleration for interactive applications.

Contribution

It presents a novel non-smooth Newton framework with a new complementarity preconditioner and GPU-based solver for efficient, robust simulation of deformable and rigid multi-body dynamics.

Findings

Supports nonlinear dynamics including hyperelastic bodies and articulated mechanisms

Achieves faster convergence with the new complementarity preconditioner

Enables real-time simulation in robotics scenarios using GPU acceleration

Abstract

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.

Tables

Table 1. Table 1 . Glossary of terms.

Symbol	Meaning
$𝐪$	Generalized system coordinates
$\dot{𝐪}$	First time derivative of generalized coordinates
$\ddot{𝐪}$	Second time derivative of generalized coordinates
$\widetilde{𝐪}$	Predicted or inertial system configuration
$𝐟$	Generalized force function
$𝐮$	Generalized velocity vector
$𝐌$	System mass matrix
$𝐊$	Geometric stiffness matrix
$𝐄$	Compliance matrix
$𝐃$	Friction force basis
$𝐖$	Friction compliance matrix
$𝐆$	Kinematic mapping transform
$𝐍$	Material stiffness matrix
$\Psi$	Strain energy density function
${𝐜}_b$	Bilateral constraint vector
${𝐜}_n$	Normal contact constraint vector
${𝐜}_f$	Frictional contact constraint vector
$𝐧$	Contact normal
$d$	Separation distance for a contact
${\varphi }_n$	Normal contact NCP constraint vector
${\varphi }_f$	Frictional contact NCP constraint vector
${𝐉}_b$	Jacobian of bilateral constraints
${𝐉}_n$	Jacobian of contact normal NCP constraints
${𝐉}_f$	Jacobian of contact friction NCP constraints
${𝝀}_b$	Lagrange multiplier vector for bilateral constraints
${𝝀}_n$	Lagrange multipliers for normal contact constraints
${𝝀}_f$	Lagrange multipliers for friction contact constraints
$n_d$	Number of system degrees of freedom
$n_b$	Number of bilateral constraints
$n_c$	Number of contacts

Table 2. Table 2 . Parameters and statistics for different examples. All scenarios have used a time-step of h = 0.0083 ℎ 0.0083 h=0.0083 s. Contact counts are for a representative step of the simulation. Performance cost is typically dominated by the linear solver step, while the matrix assembly cost is small. We have reported solver times for a single time-step of our GPU-based PCR solver. We note that for small problems the performance of iterative methods is limited by fixed cost kernel launch overhead and that scaling with problem size is typically sub-linear up to available compute resources. For contacts we also report the maximum number of contacts in a single simulation island in brackets.

Example	Bodies	Joints	Tetra	Contacts	Newton Iters	Linear Iters	Assembly (ms)		Solve (ms)
	#	#	#	# (island)	#	#	Avg	Std. Dev.	Avg	Std. Dev.
Allegro DexNet	21	21	0	104 (104)	8	20	0.4	0.1	10.0	1.5
Allegro Ball	21	21	427	409 (409)	6	50	0.35	0.05	11.4	1.2
PneuNet	1	0	2241	532 (532)	4	80	0.3	0.05	22.25	1.4
Fetch Tomato	29	28	319	371 (371)	4	50	0.45	0.05	12.2	1.25
Fetch Beam	42	42	0	300 (300)	6	50	0.2	0.1	18.6	2.25
Arch	20	0	0	134 (134)	6	20	0.15	0.05	7.75	1.4
Table Pile	54	0	0	936 (736)	4	20	0.1	0.05	3.7	0.95
Humanoid Run	5200	4800	0	10893 (33)	4	10	0.8	0.05	8.9	0.25
Yumi Cabinet	6400	6600	0	4082 (55)	5	25	1.8	0.25	56.0	6.4

Equations

$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}-\mathbf{f}(\mathbf{q},\dot{\mathbf{q}})=0,$$

$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}-\mathbf{f}(\mathbf{q},\dot{\mathbf{q}})-\nabla {\mathbf{c}}_b^T{\lambda }_b$$

$${\mathbf{c}}_b(\mathbf{q})$$

$$U(\mathbf{q})=\frac{k}{2}{C}_b(\mathbf{q}{)}^2.$$

$${\mathbf{f}}_b=-\nabla U^T=-k\nabla C_b^T{C}_b(\mathbf{q})=\nabla C_b^T{\lambda }_b.$$

$${\lambda }_b=-k{C}_b(\mathbf{q}).$$

$$C_b(\mathbf{q})+e{\lambda }_b=0,$$

$${\mathbf{c}}_b(\mathbf{q})+\mathbf{E}{\lambda }_b=0.$$

$$C_n(\mathbf{q})={\mathbf{n}}^T\left[\mathbf{a}(\mathbf{q})-\mathbf{b}(\mathbf{q})\right]-d\ge 0.$$

$$0\le C_n(\mathbf{q})\perp {\lambda }_n\ge 0,$$

$$\underset{{\lambda }_f}{\mathrm{argmin}}$$

$$|{\lambda }_f|\le \mu {\lambda }_n.$$

$$\mathcal{L}({\lambda }_f,{\lambda }_n)={\dot{\mathbf{q}}}^T\mathbf{D}{\lambda }_f+{\lambda }_s(|{\lambda }_f|-\mu {\lambda }_n).$$

$${\mathbf{D}}^T\dot{\mathbf{q}}+{\lambda }_s\text{partialderivative}|{\lambda }_f|{\lambda }_f=0$$

$$0\le {\lambda }_s\perp \mu {\lambda }_n-|{\lambda }_f|\ge 0.$$

$${\mathbf{D}}^T\dot{\mathbf{q}}+\frac{|{\mathbf{D}}^T\dot{\mathbf{q}}|}{|{\lambda }_f|}{\lambda }_f=0$$

$$0\le |{\mathbf{D}}^T\dot{\mathbf{q}}|\perp \mu {\lambda }_n-|{\lambda }_f|\ge 0.$$

$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}-\mathbf{f}(\mathbf{q},\dot{\mathbf{q}})-\nabla {\mathbf{c}}_b^T(\mathbf{q}){\lambda }_b-\nabla {\mathbf{c}}_n^T(\mathbf{q}){\lambda }_n-\nabla {\mathbf{c}}_f^T(\mathbf{q}){\lambda }_f=0$$

$${\mathbf{c}}_b(\mathbf{q})+\mathbf{E}{\lambda }_b=0$$

$$0\le {\mathbf{c}}_n(\mathbf{q})\perp {\lambda }_n\ge 0$$

$$i\in \mathcal{A},{\mathbf{D}}_i^T\dot{\mathbf{q}}+\frac{|{\mathbf{D}}_i^T\dot{\mathbf{q}}|}{|{\lambda }_{f,i}|}{\lambda }_{f,i}=0$$

$$i\in \mathcal{A},0\le |{\mathbf{D}}_i^T\dot{\mathbf{q}}|\perp {\mu }_i{\lambda }_{n,i}-|{\lambda }_{f,i}|\ge 0$$

$$i\in \mathcal{I},{\lambda }_{f,i}=0$$

$$\varphi (a,b)=0⟺0\le a\perp b\ge 0.$$

$${\varphi }_{\text{min}}(a,b)=\text{min}(a,b)=0.$$

$$0\le C_n(\mathbf{q})\perp {\lambda }_n\ge 0.$$

$${\varphi }_n(\mathbf{q},{\lambda }_n)=\text{min}(C_n(\mathbf{q}),{\lambda }_n)=0,$$

$$\text{partialderivative}{\varphi }_n\mathbf{q}$$

$$\text{partialderivative}{\varphi }_n{\lambda }_n$$

$${\varphi }_{\text{FB}}(a,b)=a+b-\sqrt{a^2+b^2}=0.$$

$$\alpha (a,b)=\text{partialderivative}{\varphi }_{FB}a$$

$$\beta (a,b)=\text{partialderivative}{\varphi }_{FB}b$$

$$\alpha (0,0)$$

$$\beta (0,0)$$

$${\varphi }_n(\mathbf{q},{\lambda }_n)={\varphi }_{\text{FB}}(C_n(\mathbf{q}),{\lambda }_n)=0,$$

$$\text{partialderivative}{\varphi }_n\mathbf{q}$$

Full text

numerical optimization, friction, contact, multi-body dynamics, robotics

††copyright: acmlicensed††journal: TOG††journalyear: 2019††journalvolume: 38††journalnumber: 5††article: 20††publicationmonth: 6††price: 15.00††doi: 10.1145/3338695††ccs: Computing methodologies Simulation by animation††ccs: Computing methodologies Interactive simulation††ccs: Computer systems organization Robotics

1. Introduction

Enabling the next-generation of robots that can, for example, enter a kitchen and prepare dinner, requires new control algorithms capable of navigating complex environments and performing dexterous manipulation of real-world objects, as illustrated in Figure 1. Machine-learning based approaches hold the promise of unlocking this capability, but require large amounts of data to work effectively [Levine et al., 2018]. For many tasks, gathering this data from the real-world may be inefficient, or impractical due to safety concerns. In contrast, simulation is relatively inexpensive, safe, and has been used to learn and transfer behaviors such as walking and jumping to real robots [Tan et al., 2018][Sadeghi et al., 2017]. Extending transfer learning beyond locomotion to a wider range of behaviors requires the robust and efficient simulation of richer environments incorporating multiple physical models [Heess et al., 2017].

We believe the computer graphics community is uniquely placed to address the simulation needs of robotics. One area that computer graphics has studied extensively is two-way coupled simulation of rigid and deformable objects. Such algorithms are necessary to simulate tasks involving dexterous manipulation of soft objects, or even soft robots themselves. An example of the latter is found in the PneuNet gripper [Ilievski et al., 2011] as shown in Figure 2. This design incorporates a flexible elastic gripper, an inextensible paper layer, and a chamber that is pressurized to cause the finger to curve. An additional example is shown in Figure 3 where a deformable ball is manipulated by a humanoid robot hand made of rigid parts. This type of multiphysics scenario is challenging for simulation because the internal forces may be large, and they must interact strongly with contact to allow robust grasping.

Methods for simulating multi-body systems in the presence of contact and friction most commonly formulate a linear complementarity problem (LCP) that is solved in one of two ways: relaxation methods such as projected Gauss-Seidel (PGS), or direct methods such as Dantzig’s pivoting algorithm. Relaxation methods are popular due to their simplicity, but suffer from slow convergence for poorly conditioned problems [Erleben, 2013]. In contrast, direct methods can achieve greater accuracy, but are typically serial algorithms that scale poorly with problem size. Finding methods that combine the simplicity and robustness of relaxation methods with the accuracy of direct methods remains a challenge in computer graphics and robotics. While solving LCP problems efficiently is still the subject of active research, as a model they may not capture all of the dynamics we wish to simulate. For example, hyperelastic materials have highly nonlinear forces that significantly affect behavior compared to linear models [Smith et al., 2018]. In addition, contact models themselves may be nonlinear particularly when considering compliance and deformation [Li and Kao, 2001].

In this work we develop a framework based on Newton’s method that solves the underlying nonlinear complementarity problems (NCPs) arising in multi-body dynamics. We combine a nonlinear contact model, articulated rigid-body model, and a hyperelastic material model as a system of semi-smooth equations, and show how it can be solved efficiently. A key advantage of our Newton-based approach is that it allows the use of off-the-shelf linear solvers as the fundamental computational building block. This flexibility means we can choose to use accurate direct solvers, or take advantage of highly-optimized iterative solvers available for parallel architectures such as graphics-processing units (GPUs). In summary, our main contributions are:

• •

A formulation of smooth isotropic Coulomb friction in terms of non-smooth complementarity functions.

• •

A new complementarity preconditioner that significantly improves convergence for contact problems.

• •

A generalized compliance formulation that supports hyperelastic material models.

• •

A simple approximation of geometric stiffness to improve robustness without changing system dynamics.

2. Related Work

2.1. Contact

The seminal work of Jean and Moreau [1992] introduced an implicit time-stepping scheme for contact problems with friction. This work was further popularized by Stewart and Trinkle [1996] who linearize the Coulomb friction cone and solve LCPs using Lemke’s method in a fixed-point iteration to handle nonlinear forces. We also make use of a fixed-point iteration, but in contrast to their work we use non-smooth functions to model nonlinear friction cones. Kaufman et al. [2008] proposed a method for implicit time-stepping of contact problems by solving two separate QPs for normal and frictional forces in a staggered manner. In our work we do not stagger our system updates but solve for contact forces and friction forces in a combined system, and without using a linearized cone. In addition, rather than using QP solves our method requires only the solution of a symmetric linear system as a building block, allowing for the use of off-the-shelf linear solvers.

Relaxation methods such as projected Gauss-Seidel are popular in computer graphics thanks to their simplicity and robustness [Bender et al., 2014]. While robust, these methods suffer from slow convergence for poorly conditioned problems, e.g.: those with high-mass ratios [Erleben, 2013]. In addition, Gauss-Seidel iteration suffers from order dependence and is challenging to parallelize, while Jacobi methods require modifications to ensure convergence [Tonge et al., 2012]. Daviet et al. [2011] used a change of variables to restate the Coulomb friction cone into a self-dual complementarity cone problem followed by a modified Fischer-Burmeister reformulation to obtain a local non-smooth root search problem. In comparison, we work directly with the friction cone as limit surfaces as we believe this provides us with more modeling freedom, e.g.: for anisotropic or non-symmetric friction cones. Another difference is that we solve for all contacts simultaneously rather than one-by-one.

Otaduy et al. [2009] presented an implicit time-stepping scheme for deformable objects that solves a mixed linear complementarity problem (MLCP) using a nested Gauss-Seidel relaxation over primal and dual variables. Prior work on simulating smooth friction models has used proximal-map projection operators [Jourdan et al., 1998; Jean, 1999; Erleben, 2017] which work by projecting contact forces to the friction cone one contact at a time until convergence. There has been considerable work to address the slow convergence of relaxation methods, Mazhar et al. [2015] use a convexification of the frictional contact problem [Anitescu and Hart, 2004] to obtain a cone complementarity problem (CCP) and solve it using an accelerated version of projected gradient descent. Silcowitz et al. [2009; 2010] developed a method for solving LCPs based on non-smooth nonlinear conjugate gradient (NNCG) applied to a PGS iteration. Francu et al. [2015] proposed an improved Jacobi method based on a projected conjugate residual (CR) method. We also make use of Krylov space linear solvers, however our Newton-based iteration is decoupled from the underlying linear backend, allowing the application of matrix-splitting relaxation methods, or even direct solvers.

Early work on Newton methods for contact problems used a formulation based on a generalized projection operator and an augmented Lagrangian approach to unilateral constraints [Alart and Curnier, 1991; Curnier and Alart, 1988]. Newton-based approaches found in libraries such as PATH [Dirkse and Ferris, 1995] have proved successful in practice, and have been applied to smooth Coulomb friction models in fiber assemblies in computer graphics [Bertails-Descoubes et al., 2011; Daviet et al., 2011; Kaufman et al., 2014]. These approaches formulate the complementarity problem in terms of non-smooth functions, and solve them with a generalized version of Newton’s method. This approach can yield quadratic convergence, although with a higher per-iteration cost than relaxation methods.

Todorov [2010] observed that if solving a nonlinear time-stepping problem, e.g.: due to an implicit time-discretization, it makes little sense to perform the friction cone linearization, since the smooth contact model can be treated simply as an additional set of nonlinear equations in the system. This observation is at the heart of our method, but in contrast to their work we formulate friction in terms of arbitrary NCP-functions, and extend our framework to handle deformable bodies. Our approach is based on Newton’s method, and combined with our proposed preconditioner we show that it enables handling considerably more ill-conditioned problems than relaxation methods, while naturally accommodating nonlinear friction models. For a review of numerical methods for linear complementarity problems we refer to the book by Niebe and Erleben [2015]. For a review of non-smooth methods applied to dynamics problems we refer to the book by Acary & Brogliato [2008].

2.2. Coupled Systems

There has been considerable work in computer graphics on coupling between rigid and deformable bodies. Shinar et al. [2008] proposed a method for the coupled simulation of rigid and deformable bodies through a multi-stage update that peforms collisions, contacts, and stabilization in separate passes. In contrast, we formulate a system update that includes elastic and contact dynamics in a single phase. Duriez [2013] showed real-time control of soft-robots using a co-rotational FEM model and Gauss-Seidel based constraint solver. Servin et al. [2006] introduced a compliant version of elasticity that fits naturally inside constrained rigid body simulators. Their work was extended by Tournier et al. [2015] who include a geometric stiffness term to improve stability. They use temporal coherence of Lagrange multipliers to build the system Jacobian and compute a Cholesky decomposition, followed by a projected-Gauss Seidel solve for contact. For smaller problems our approach is compatible with direct solvers, however we avoid the requirement of dense matrix decompositions by using a diagonal geometric stiffness approximation inspired by the work of Andrews et al. [2017], this improves stability and allows us to apply iterative methods.

In this work we propose a generalized view of compliance, and give a recipe for constructing the compliance form of an arbitrary material model given its strain-energy density function in terms of principle stretches, strains, or other parameterization. As an example we show how to formulate the stable Neo-Hookean material proposed by Smith et al. [2018]. Liu et al. [2016] propose a quasi-Newton method for hyperelastic materials based on Projective Dynamics [Bouaziz et al., 2014] and model contact through stiff penalty forces. In contrast we model contact through complementarity constraints that naturally fit into existing multi-body simulations.

3. Background

At the most general level, the dynamics of the systems we wish to simulate can be described through the following second order differential equation:

[TABLE]

where for a system with $n_{d}$ degrees of freedom, $\mathbf{M}\in\mathbb{R}^{n_{d}\times n_{d}}$ is the system mass matrix, $\mathbf{f}\in\mathbb{R}^{n_{d}}$ is a generalized force vector, and $\mathbf{q}\in\mathbb{R}^{n_{d}}$ is a vector of generalized coordinates describing the system configuration with $\dot{\mathbf{q}},\ddot{\mathbf{q}}$ the first and second time derivatives respectively. We refer the reader to Table 1 for a list of symbols used in this paper.

3.1. Hard Kinematic Constraints

We impose hard bilateral (equality) constraints on the system configuration through constraint functions of the form $C_{b}(\mathbf{q})=0$. Using D’Alembert’s principle [Lanczos, 1970] the force due to such a constraint is of the form $\mathbf{f}_{b}=\nabla C_{b}(\mathbf{q})^{T}\lambda_{b}$ where $\nabla C_{b}$ is the constraint’s gradient with respect to the system coordinates, and $\lambda_{b}$ is a Lagrange multiplier variable used to enforce the constraint. Combining these algebraic constraints with the differential equation (1) gives the following Differential Algebraic Equation (DAE):

[TABLE]

For a system with $n_{b}$ equality constraints we define the constraint vector $\mathbf{c}_{b}(\mathbf{q})=[C_{b,1},C_{b,2},\cdots,C_{b,n_{b}}]^{T}\in\mathbb{R}^{n_{b}}$, $\nabla\mathbf{c}_{b}\in\mathbb{R}^{n_{b}\times n_{d}}$ its gradient with respect to system coordinates, and ${\bf\it\lambda}_{b}\in\mathbb{R}^{n_{b}}$ the vector of Lagrange multipliers. In general bilateral constraints may depend on velocity and time, we assume scleronomic constraints for the remainder of the paper to simplify the exposition.

3.2. Soft Kinematic Constraints

In addition to hard constraints, we make extensive use of the compliance form of elasticity [Servin et al., 2006]. This may be derived by considering a quadratic energy potential defined in terms of a constraint function $C_{b}(\mathbf{q})$ and stiffness parameter $k>0$

[TABLE]

The force due to this potential is then given by

[TABLE]

Here we have introduced the variable $\lambda_{b}$, defined as

[TABLE]

We can move all terms to one side to write this as a constraint equation,

[TABLE]

where $e=k^{-1}$ is the compliance, or inverse stiffness coefficient. We can incorporate this into our equations of motion by defining the diagonal matrix $\mathbf{E}=\text{diag}(e_{1},e_{2},\cdots,e_{n_{b}})\in\mathbb{R}^{n_{b}\times n_{b}}$, and updating (3) as follows,

[TABLE]

This form is mathematically equivalent to including quadratic energy potentials defined in terms of a stiffness. The benefit of the compliance form is that it remains numerically stable even as $k\to\infty$ [Tournier et al., 2015]. We describe how to extend this model to handle general energy potentials in Section 9.

3.3. Contact

We treat contacts as inelastic and prevent interpenetration between bodies through inequality constraints of the form

[TABLE]

Here $\mathbf{n}\in\mathbb{R}^{3}$ is the contact plane normal. We use the normalized direction vector between closest points of triangle-mesh features as the normal vector. The points $\mathbf{a}$ and $\mathbf{b}\in\mathbb{R}^{3}$ may be points defined in terms of a rigid body frame, or particle positions in the case of a deformable body. The constant $d$ is a thickness parameter that represents the distance along the contact normal we wish to maintain. Non-zero values of $d$ can be used to model shape thickness, as illustrated in Figure 4. The normal force for a contact is given by $\mathbf{f}_{n}=\nabla C_{n}^{T}\lambda_{n}$, with the additional Signorini-Fichera complementarity condition

[TABLE]

which ensures contact forces only act to separate objects [Stewart, 2000]. We treat the contact normal $\mathbf{n}$ as fixed in world-space. For a system with $n_{c}$ contacts, we define the vector of contact constraints as $\mathbf{c}_{n}=[C_{n,1},\cdots,C_{n,n_{c}}]\in\mathbb{R}^{n_{c}}$, their gradient $\nabla\mathbf{c}_{n}\in\mathbb{R}^{n_{c}\times n_{d}}$, and the associated Lagrange multipliers as ${\bf\it\lambda}_{n}\in\mathbb{R}^{n_{c}}$. In our implementation contact constraints are created when body features come within some fixed distance of each other. This approach works well for reasonably small time-steps, but can lead to over-constrained configurations. More sophisticated non-interpenetration constraints can be formulated to avoid this problem [Williams et al., 2017].

3.4. Friction

We derive our friction model from a principle of maximal dissipation that requires the frictional forces remove the maximum amount of energy from the system while having their magnitude bounded by the normal force [Stewart, 2000]. For each contact we define a two-dimensional basis $\mathbf{D}=\mathbf{\Gamma}(\mathbf{q})[\mathbf{d}_{1}\mathbf{d}_{2}]\in\mathbb{R}^{n_{d}\times 2}$. The contact basis vectors $\mathbf{d}_{1},\mathbf{d}_{2}\in\mathbb{R}^{3\times 1}$ are lifted from spatial to generalized coordinates by the transform $\mathbf{\Gamma}\in\mathbb{R}^{n_{d}\times 3}$ using the notation of Kaufman et al. [2014]. The generalized frictional force for a single contact is then $\mathbf{f}_{f}=\mathbf{D}{\bf\it\lambda}_{f}$, where ${\bf\it\lambda}_{f}\in\mathbb{R}^{2}$ is the solution to the following minimization

[TABLE]

Here $\mu$ is the coefficient of friction, and $\lambda_{n}$ is the Lagrange multiplier for the normal force at this contact which for the moment we assume is known. This minimization defines an admissible cone that the total contact force must lie in, as illustrated in Figure 5. The Lagrangian associated with this minimization is

[TABLE]

where $\lambda_{s}$ is a slack variable used to enforce the Coulomb constraint that the friction force magnitude is bounded by the $\mu$ times the normal force magnitude. When $\mu\lambda_{n}>0$ the problem satisfies Slater’s condition [Boyd and Vandenberghe, 2004] and we can use the first-order Karush-Kuhn-Tucker (KKT) conditions for (11) given by

[TABLE]

Equation (13) requires that the frictional force act in a direction opposite to any relative tangential velocity. Equation (14) is a complementarity condition that describes the sticking and slipping behavior characteristic of dry friction. In the case that $\mu\lambda_{n}=0$ Slater’s condition is not satisfied, however in this case the only feasible point is ${\bf\it\lambda}_{f}=\mathbf{0}$, which we handle explicitly. We now make a transformation to remove the slack variable $\lambda_{s}$. To do this we observe that the derivative of a vector’s length is the vector normalized, so we have $\partialderivative{|{\bf\it\lambda}_{f}|}{{\bf\it\lambda}_{f}}={\bf\it\lambda}_{f}/|{\bf\it\lambda}_{f}|$, which by definition is a unit-vector. We can use this fact to express the optimality conditions in terms of $\dot{\mathbf{q}}$ and ${\bf\it\lambda}_{f}$ only:

[TABLE]

Finally we combine the frictional basis vectors and multipliers for all contacts in the system as a single matrix $\nabla\mathbf{c}_{f}=[\mathbf{D}_{1},\cdots,\mathbf{D}_{n_{c}}]^{T}\in\mathbb{R}^{2n_{c}\times n_{d}}$ and vector ${\bf\it\lambda}_{f}=[{\bf\it\lambda}_{f,1},\cdots,{\bf\it\lambda}_{f,n_{c}}]^{T}\in\mathbb{R}^{2n_{c}}$.

3.5. Governing Equations

Assembling the above components, our continuous equations of motion are given by the following nonlinear system of equations:

[TABLE]

where $\mathcal{A}=\{i\in(1,\cdots,n_{c})\mid\mu_{i}\lambda_{n,i}>0\}$ is the set of all contact indices where the normal contact force is active, and $\mathcal{I}=\{i\in(1,\cdots,n_{c})\mid\mu_{i}\lambda_{n,i}\leq 0\}$ is its complement. This combination of a differential equation with equality and complementarity conditions is known as Differential Variational Inequality (DVI) [Stewart, 2000]. The coupling between normal and frictional complementarity problems makes the problem non-convex, and in the case of implicit time-integration leads to an NP-hard optimization problem [Kaufman et al., 2008]. In the next section we develop a practical method to solve this problem.

4. Nonlinear Complementarity

One successful approach [Ferris and Munson, 2000] to solving nonlinear complementarity problems is to reformulate the problem in terms of a NCP-function whose roots satisfy the original complementarity conditions, i.e.: functions where the following equivalence holds

[TABLE]

Combined with an appropriate time-discretization, such a NCP-function turns our DVI problem into a root finding one. In general the functions $\phi$ are non-smooth, but allow us to apply a wider range of numerical methods [Munson et al., 2001].

4.1. Minimum-Map Formulation

The first such function we consider is the minimum-map defined as

[TABLE]

The equivalence of this function to the original NCP can be verified by examining the values associated with each conditional case [Cottle, 2008]. We now show how this reformulation applies to our unilateral contact constraints. Recall that the complementarity condition associated with a contact constraint $C_{n}(\mathbf{q})$ and its associated Lagrange multiplier $\lambda_{n}$ is

[TABLE]

We can write this in the equivalent minimum-map form as

[TABLE]

which has the following derivatives,

[TABLE]

From these cases we can see that the minimum-map gives rise to an active-set style method where a contact is considered active if $C_{n}(\mathbf{q})\leq\lambda_{n}$. Active contacts are treated as equality constraints, while for inactive contacts the minimum-map enforces that the constraint’s Lagrange multiplier is zero.

4.2. Fischer-Burmeister Formulation

An alternative NCP-function is given by Fischer-Burmeister [1992], who observe the roots of the following equation satisfy complementarity:

[TABLE]

The Fischer-Burmeister function is interesting because, unlike the minimum-map, it is smooth everywhere apart from the point $(a,b)=(0,0)$. For $(a,b)\neq(0,0)$ the partial derivatives of the Fisher-Burmeister function are given by:

[TABLE]

At the point $(a,b)=(0,0)$ the derivative is set-valued (see Figure 8). For Newton methods it suffices to choose any value from this subgradient. Erleben [2013] compared how the choice of derivative at the non-smooth point affects convergence for LCP problems and found no overall best strategy. Thus, for simplicity we make the arbitrary choice of

[TABLE]

For a contact constraint $C_{n}$, with Lagrange multiplier $\lambda_{n}$ we may then define our contact function alternatively as,

[TABLE]

with derivatives given by

[TABLE]

4.3. Frictional Constraints

We now formulate our friction model in terms of non-smooth functions. Our first step is to rewrite the optimality conditions (15)-(16) in terms of an NCP-function. For a single contact with $\mu\lambda_{n}>0$ the following must hold:

[TABLE]

where $\psi_{f}$ may be either the minimum-map or Fischer-Burmeister function. We now introduce a new quantity with the goal of linearizing the relationship between $\mathbf{D}^{T}\dot{\mathbf{q}}$ and ${\bf\it\lambda}_{f}$ such that upon convergence of a fixed-point iteration the initial complementarity problem is solved. We do this by defining the following fixed-point iteration based on the relationship between $|\mathbf{D}^{T}\dot{\mathbf{q}}|$ and $|{\bf\it\lambda}_{f}|$ given by $\psi_{f}$,

[TABLE]

By construction, a fixed-point of this iteration satisfies the original complementarity condition. We use these expressions as replacements for $\frac{|\mathbf{D}^{T}\dot{\mathbf{q}}|}{|{\bf\it\lambda}_{f}|}$ in (31), by defining the following coefficient,

[TABLE]

Inserting this into (31) allows us to write it in the compact form,

[TABLE]

Here $W$ can be interpreted as acting as a compliance term that works to project the frictional force back onto the friction cone as illustrated in Figure 6. The derivatives of ${\bf\it\phi}_{f}$ are then given by,

[TABLE]

where we have treated $W$ as a constant. Replacing the complementarity condition by a fixed-point iteration ensures that we obtain a symmetric system of equations in the following section. In Appendix A we derive the exact form of $W$ for both minimum-map and Fischer-Burmeister functions.

We note that conical equivalents of the Fischer-Burmeister function exist and have been used to model smooth isotropic friction [Fukushima et al., 2002; Daviet et al., 2011]. One advantage of our method being based on a fixed-point iteration is that it can be extended to arbitrary friction surfaces for e.g.: anisotropic or even non-symmetric friction cones.

5. Implicit Time-Stepping

The continuous equations of motion are expressed purely in terms of our generalized coordinates $\mathbf{q}$, their derivatives, and the Lagrange multipliers. Although it is possible to discretize and solve these equations for $\mathbf{q}$ directly, it is convenient to introduce an additional re-parameterization in terms of a generalized velocity vector $\mathbf{u}$:

[TABLE]

Here $\mathbf{G}$ is referred to as the kinematic map, and its components depend on the choice of system coordinates. For simple particles $\mathbf{G}$ is an identity transform, however in Section 9.4 we discuss the mapping of a rigid body’s angular velocity to the corresponding quaternion time-derivative. We use the the kinematic map $\mathbf{G}$ to obtain the following mass matrix,

[TABLE]

and define the bilateral constraint Jacobians with respect to the generalized velocity $\mathbf{u}$ as follows:

[TABLE]

We group the normal and friction NCP-functions for all contacts into two vectors,

[TABLE]

and likewise define their Jacobians with respect to the generalized velocity as

[TABLE]

Using a first-order backwards time-discretization of $\dot{\mathbf{q}}$ gives the following update for the system’s generalized coordinates in terms of generalized velocities over a time-step of length $h$,

[TABLE]

The superscripts $+,-$ represent a variable’s value at the end and beginning of the time-step respectively. Discretizing our continuous equations of motion from the previous section gives the implicit time-stepping equations,

[TABLE]

Here $\mathbf{u}^{+},{\bf\it\lambda}^{+}$ are the unknown velocities and multipliers at the end of the time-step. The constant $\tilde{\mathbf{u}}=\mathbf{u}^{-}+h\mathbf{G}^{T}\mathbf{M}^{-1}\mathbf{f}(\mathbf{q}^{-},\dot{\mathbf{q}}^{-})$ is the unconstrained velocity that includes the external and gyroscopic forces integrated explicitly. Observe that through this time-discretization the original force level model of friction has changed into an impulsive model.

We highlight a few differences from common formulations. First, the equality and inequality constraints have not been linearized through an index reduction step [Hairer and Wanner, 2010]. Index reduction is a common practice that reduces the order of the DAE by solving only for the constraint time-derivatives, e.g.: $\dot{\mathbf{c}}_{b}=\mathbf{0}$. Index reduction results in a linear problem, but requires additional stabilization terms to combat drift and move the system back to the constraint manifold. These stabilization terms are often applied as the equivalent of penalty forces [Ascher et al., 1995] and are known as a source of instability and tuning issues. Work has been done to add additional post-stabilization passes [Cline and Pai, 2003], however these require solving nonlinear projection problems, which gives up the primary benefit of performing the index reduction step. Our discrete equations of motion are also nonlinear, but they require no additional stabilization terms or projection passes.

A second point we highlight is that the friction cone defined through (11) has not been linearized into a faceted cone as is commonly done [Stewart and Trinkle, 2000; Kaufman et al., 2008]. The faceted cone approximation provides a simpler linear complementarity problem (LCP), but increases the number of Lagrange multipliers required (one per-facet), and introduces an approximation bias where frictional effects are stronger in some directions than others. In the next section we propose a method that solves the NCP corresponding to smooth isotropic friction using a fixed number of Lagrange multipliers.

The implicit time-discretization above treats contact forces as impulsive, meaning they are able to prevent sliding and interpenetration instantaneously (over a single time-step in the discrete setting). This avoids the inconsistency raised by Painlevé in the continuous setting.

6. Non-Smooth Newton

We now develop a method to solve the discretized equations of motion (46)-(50). Our starting point is Newton’s method which we briefly review here. Given a set of nonlinear equations $\mathbf{r}(\mathbf{x})=\mathbf{0}$ in terms of a vector valued variable $\mathbf{x}$, Newton’s method gives a fixed point iteration in the following form:

[TABLE]

where $n$ is the Newton iteration index. Newton’s method chooses $\mathbf{A}$ specifically to be the matrix of partial derivatives evaluated at the current system state, i.e.:

[TABLE]

Although Newton’s method is most commonly used for solving systems of smooth functions it may also be applied to non-smooth functions by generalizing the concept of a derivative [Pang, 1990]. Qi and Sun [1993] showed that Newton will converge for non-smooth functions provided $\mathbf{A}\in\partial\mathbf{r}(\mathbf{x}^{k})$, where $\partial\mathbf{r}$ is the generalized Jacobian of $\mathbf{r}$ at $\mathbf{x}^{k}$ defined by Clarke [1990]. Intuitively, this is the convex hull of all directional derivatives at the non-smooth point, as illustrated in Figure 8. The derivatives of our NCP-functions given in the previous section belong to this subgradient, and we can use them in the fixed-point iteration of (51) directly. Algorithms of this type are sometimes referred to as semismooth methods, we refer the reader to the article by Hintermüller [2010] for a more mathematical introduction and the precise definition of semismoothness. In many cases $\mathbf{A}$ is not inverted directly, and the linear system for $\Delta\mathbf{x}$ may only be solved approximately. When this is the case we refer to the method as an inexact Newton method. Additionally, when the partial derivatives in $\mathbf{A}$ are also approximated we refer to it as a quasi-Newton method. In our method we make use of both approximations.

6.1. System Assembly

Linearizing our discrete equations of motion (46)-(50) each Newton iteration provides the following linear system in terms of the change in system variables:

[TABLE]

where $\mathbf{K}$ is the geometric stiffness matrix arising from the spatial derivatives of constraint forces discussed in Section 8.3. The lower-diagonal blocks are the derivatives of our contact functions with respect to their Lagrange multipliers

[TABLE]

Grouping the sub-block components such that

[TABLE]

we can write the system compactly as,

[TABLE]

The right-hand side is given by evaluating our discrete equations of motion (46)-(50) at the current Newton iterate. Here, $\mathbf{g}$ is our momentum balance equation,

[TABLE]

and $\mathbf{h}$ is a vector of constraint errors,

[TABLE]

Here the left-hand side matrix should be considered acting block-wise on the constraint error. Since the frictional constraints are measured at the velocity level they are not scaled by $\frac{1}{h}$ like the positional constraints.

The mass block matrix $\mathbf{M}$ is evaluated at the beginning of the time-step using $\mathbf{q}^{-}$, while the $\mathbf{H}$ matrix is evaluated each Newton iteration as described in Section 8.3. The friction compliance block, $\mathbf{W}$, is evaluated at each Newton iteration using the current Lagrange multipliers. This means the frictional forces are bounded by the normal force from the previous Newton iteration. This is similar to the fixed point iteration described by Stewart and Trinkle [2000] but with a smooth friction cone. It is also similar in principle to the Kaufman et al.’s Staggered Projections [2008], with the difference that we update both normal and frictional forces during a single Newton iteration, rather than solving two separate minimizations. For inactive contacts with $\lambda_{n}\geq 0$ we conceptually disable their frictional constraint equation rows by removing them from the system. In practice this can be performed by zeroing the corresponding rows to avoid changing the system matrix structure.

6.2. Schur-Complement

The system (56) is a saddle-point problem [Benzi et al., 2005] that is indefinite and possibly singular. To obtain a reduced positive semi-definite system, we take the Schur-complement with respect to the mass-block to obtain

[TABLE]

After solving this system for $\Delta{\bf\it\lambda}$ the velocity change $\Delta\mathbf{u}$ may be evaluated directly,

[TABLE]

and the system updated accordingly,

[TABLE]

Here $t$ is a step-size determined by line-search or other means (see Section 8). We refer to (59) as the Newton compliance formulation. Under certain conditions we could alternatively have taken the Schur-complement of (56) with respect to the compliance block $\mathbf{C}$, instead of the mass block $\mathbf{H}$. This transformation is only possible if $\mathbf{C}$ is non-singular, but it leads to what we call the Newton stiffness formulation,

[TABLE]

This form is closely related to that of Projective Dynamics [Bouaziz et al., 2014], and arises from a linearization of the elastic forces due to a quadratic energy potential. Having $\mathbf{C}$ be invertible corresponds to having compliance everywhere, or in other words, no perfectly hard constraints. For stiff materials, where $\mathbf{C}$ is poorly conditioned, this approach leads to numerical problems in calculating $\mathbf{C}^{-1}$. However, if the system has fewer degrees of freedom than constraints this transformation can result in a smaller system. In this work we are interested in methods that combine stiff constraints with deformable bodies, so we use the compliance form which accommodates both.

7. Complementarity Precondtioning

An interesting property of the non-smooth complementarity formulations is that their solutions are invariant up to a constant positive scale $r$ applied to either of the arguments. Specifically, a solution to $\phi(a,b)=0$, is also a solution to the scaled problem, $\phi(a,rb)=0$. Alart [1997] presented an analysis of the optimal choice of $r$ in the context of linear elasticity in a quasistatic setting. Erleben [2017] also explored this free parameter in the context of proximal-map operators for rigid bodies.

In this section we propose a new complementarity preconditioning strategy to improve convergence by choosing $r$ based on the system properties and discrete time-stepping equations. To motivate our method, we make the observation that the two sides of the complementarity condition typically have different physical units. For example, a contact NCP function

[TABLE]

has the units of meters for the first parameter, and units of Newtons for the second parameter. This mismatch can lead to poor convergence in a manner similar to the effect of row-scaling in traditional iterative linear solvers.

Inspired by the use of diagonal preconditioners for linear solvers, our idea is to use the effective system mass and time-stepping equations to put both sides of the complementarity condition into the same units. The appropriate $r$-factor to perform this scaling comes from the relation between ${\bf\it\lambda}$ and our constraint functions, given through the equations of motion. Specifically, for a unilateral constraint with index $i$ we set $r_{i}$ to be,

[TABLE]

Intuitively, this is the time-step scaled effective mass of the system, and relates how a change in Lagrange multiplier affects the corresponding constraint value. This has the effect of making both sides of the complementarity equation have the same slope, as illustrated in Figure 9. For a friction constraint $\phi_{f}$ the correct scaling factor is $r_{i}=h\left[\mathbf{J}\tilde{\mathbf{M}}^{-1}\mathbf{J}^{T}\right]_{ii}$ since this is a velocity-force relationship. When using a Jacobi preconditioner the diagonal of the system matrix is already computed, so applying $r_{i}$ to the NCP function incurs little overhead. We discuss the effect of preconditioning strategies in Section 10.2.

8. Robustness

8.1. Line Search and Starting Iterate

We implemented a back-tracking line-search based on a merit function defined as the $L_{2}$ norm of our residual vector. For frictionless contact problems we found this worked well to globalize the solution. However, for frictional problems we found line search would often cause the iteration to stall and make no progress. We believe this is related to the fact that the frictional problem is non-convex and our search direction is not necessarily a descent direction. Instead, all our results use a simple damped Newton strategy that accepts a constant fraction of the full Newton step with a factor $t\in(0,1)$. This strategy may cause a temporary increase in the merit function, but can result in overall better convergence [Maratos, 1978]. Watchdog strategies [Nocedal and Wright, 2006] may be employed to make stronger convergence guarantees, however we did not find them necessary. We have used a value of $t=0.75$ for all examples unless otherwise specified. This value is ad-hoc, but we have not found our method to be particularly sensitive to its setting.

Starting iterates have a strong effect on most optimization methods, and ours is no different. Although it is common to initialize solvers with the unconstrained velocity at the start of the time-step we found the most robust method was to use the zero-velocity solution as a starting iterate, i.e.:

[TABLE]

This choice is robust in the case of reinforcement-learning, where large random external torques are applied to bodies that can lead to initial points far from the solution. Warm-starting for constraint forces is possible, and our tests indicate this can give a good improvement in efficiency, however due to the additional book-keeping we have not used warm-starting in our reported results.

8.2. Preconditioned Conjugate Residual

A key advantage of our non-smooth formulation is that it allows black-box linear solvers to be used for nonlinear complementarity problems. Nevertheless, the choice of solver is still an important factor that affects performance and robustness. A common issue we encountered is that even simple contact problems may lead to singular Newton systems. Consider for example a tessellated cylinder resting flat on the ground. In a typical simulator a ring of contact points would be created, leading to an under-determined system, i.e.: there are infinitely many possible contact force magnitudes that are valid solutions. This indeterminacy results in a singular Newton system and causes problems for many common linear solvers. The ideal numerical method should be insensitive to these poorly conditioned problems.

In addition, we seek a method that allows solving each Newton system only inexactly. Since the linearized system is only an approximation it would be wasteful to solve it accurately far from the solution. Thus, another trait we seek is a method that smoothly and monotonically decreases the residual so that it can be terminated early e.g.: after a fixed computational time budget. Finally, since our application often requires real-time updates we also look for a method that is amenable to parallelization.

The conjugate gradient (CG) method [Hestenes and Stiefel, 1952] is a popular method for solving linear systems in computer graphics. It is a Krylov space method that minimizes the quadratic energy $e=\frac{1}{2}\mathbf{x}^{T}\mathbf{A}\mathbf{x}-\mathbf{b}^{T}\mathbf{x}$. One side effect of this structure is that the residual $r=|\mathbf{A}\mathbf{x}-\mathbf{b}|$ is not monotonically decreasing. This behavior leads to problems with early termination since, although a given iterate may be closer to the solution, the true residual may actually be larger. This manifests as constraint error changing unpredictably between iterations.

A related method that does not suffer from this problem is the conjugate residual (CR) algorithm [Saad, 2003]. It is a Krylov space method similar to conjugate gradient (CG), however, unlike CG each CR iteration $k$ minimizes the residual

[TABLE]

within the Krylov space $\mathcal{K}_{k}=\text{span}\{\mathbf{b},\mathbf{A}\mathbf{b},\mathbf{A}^{2}\mathbf{b},\cdots,\mathbf{A}^{k-1}\mathbf{b}\}$. A remarkable side effect of this minimization property is that CR is monotonically decreasing in the residual norm $r=|\mathbf{A}\mathbf{x}-\mathbf{b}|$, and monotonically increasing in the solution variable norm [Fong and Saunders, 2012]. These are exactly the properties we would like for an inexact Newton method. From a computational viewpoint it is nearly identical to CG, requiring one matrix–vector multiply, and two vector inner products per-iteration. CR requires an additional vector multiply-add per-iteration, but each stage is fully parallelizable, and in our experience we did not observe any performance difference to standard CG.

In Figure 11 we compare the convergence of four common linear solvers, Jacobi, Gauss-Seidel, preconditioned conjugate gradient (PCG), and preconditioned conjugate residual (PCR). We use a diagonal preconditioner for both PCR and PCG. Our surprising finding is that PCR often has an order of magnitude lower residual for the same iteration count compared to other methods. A similar result was reported by Fong et al. [2012] in the numerical optimization literature. For symmetric positive definite (SPD) systems CR will generate the same iterates as MINRES, however it does not handle semidefinite or indefinite problems in general. In practice we observed that CR will converge for close to singular contact systems when CG and even direct Cholesky solvers may fail. We evaluate and discuss the effect of linear solver on a number of test cases in Section 10.4.

8.3. Geometric Stiffness

The upper-left block of the system matrix, $\mathbf{H}=\mathbf{M}-h^{2}\mathbf{K}$, consists of the mass matrix and a second term $\mathbf{K}$, referred to as the geometric stiffness of the system. This is defined as:

[TABLE]

$\mathbf{K}$ is not a physical quantity in the sense that it does not appear in the continuous or discrete equations of motion. It appears only as a side-effect of the numerical method, in this case the linearization performed at each Newton iteration. In practice $\mathbf{K}$ is often dropped from the system matrix, leading to a quasi-Newton method. Tournier et al. [2015] observed that ignoring geometric stiffness can lead to instability, as the Newton search direction is free to move the iterate away from the constraint manifold in the transverse direction, eventually causing the iteration to diverge. However, including $\mathbf{K}$ directly is also problematic because the mass block is no longer simple block-diagonal, and in addition, the possibly indefinite constraint Hessian restricts the numerical methods that can be applied. Inspired by the work of Andrews et al. [2017] we introduce an approximation of geometric stiffness that improves the robustness of our Newton method without these drawbacks.

We now present a method for approximating $\mathbf{H}$ that is related to the method of Broyden [1965]. This method builds the Jacobian iteratively through successive finite differences. We apply this idea to build a diagonal approximation to the geometric stiffness matrix. Using the last two Newton iterates, we can define the following first-order approximation:

[TABLE]

where $\mathbf{H}$ is the unknown matrix we wish to find, and $\Delta\mathbf{u}=\mathbf{u}^{n+1}-\mathbf{u}^{n}$ is the known difference between the last two iterates. This problem is under-determined since $\mathbf{H}$ has $n_{d}\times n_{d}$ entries, while (71) only provides $n_{d}$ equations, however if we assume $\mathbf{H}$ has the special form:

[TABLE]

where $\tilde{\mathbf{K}}=\text{diag}[c_{1},c_{2},\cdots,c_{n}]$ is a diagonal approximation to $\mathbf{K}$, then the individual entries of $\tilde{\mathbf{K}}$ are easily determined from (71) by examining each entry:

[TABLE]

Each Newton iteration we then update the mass block’s diagonal entries to form our approximate $\tilde{\mathbf{H}}$ as

[TABLE]

For most models $\mathbf{M}$ is block-diagonal, e.g.: 3x3 blocks for rigid body inertias, and so $\mathbf{H}^{-1}$ is may be easily computed to form the Schur complement. To ensure $\mathbf{H}$ remains positive-definite we clamp the shift to be positive. (59). This diagonal approximation is inspired by the method presented by Andrews et al. [2017] however this approach has several advantages. First, we do not have to derive or explicitly evaluate the constraint Hessians for different constraint types. This means our method works automatically for contacts and complex deformable models. Second, we do not need to track the Lagrange multipliers from the previous frame, as our method updates the geometric stiffness at each Newton iteration. Finally, since our method does not change the solution to the problem it does not introduce any additional damping provided the solver is run to sufficient convergence. Our approach bears considerable resemblance to a L-BFGS method [Nocedal, 1980], however the use of a diagonal update means the system matrix structure is constant, allowing us to use the efficient block-diagonal inverse for the Schur complement.

8.4. Regularization

Due to contact and constraint redundancy it is often the case that the system matrix is singular, e.g.: for a chair resting with four legs on the ground the contact constraints are linearly dependent, meaning the system is underdetermined, and the Newton matrix becomes singular. To improve conditioning we optionally apply an $\epsilon$ identity shift to the system matrix (59). This type of regularization bears some resemblance to implicit penalty methods, however unlike penalty-based approaches it does not change the underlying physical model. The regularization only applies to the error at each Newton iteration, and the solution approaches the original one as the Newton iterations progress. For the examples shown here we set $\epsilon=10^{-6}$, which is equivalent to adding a small amount of compliance to the system matrix.

9. Modeling

In the following sections we show how to generalize the compliance form of elasticity and discuss some of the physical models we have used in our examples.

9.1. Generalized Compliance

To extend our compliance form of elasticity to arbitrary material models and forces we make a new interpretation of the compliance transformation as a force factorization. Specifically, given a force $\mathbf{f}(\mathbf{q})$ we factor it as:

[TABLE]

where $\mathbf{J}$ and $\mathbf{c}$ may be scalar, vector, or matrix valued functions. The only requirements we place on the factors are that

[TABLE]

where $\mathbf{N}(\mathbf{q})$ is an invertible linear mapping. The goal in choosing a factorization is to separate the force Jacobian into components where one is easily, even analytically, invertible. This splitting moves part of the force Jacobian out of the mass matrix, and onto the compliance matrix, where it becomes a regularization term. To perform the splitting we introduce a variable ${\bf\it\lambda}=-\mathbf{c}$ which allows us to write our equations of motion (simplified for illustration) as,

[TABLE]

In the context of a Newton solve, we require the linearization of (78) which is given by

[TABLE]

Here we have used the condition that $\partialderivative{\mathbf{c}}{\mathbf{q}}=\mathbf{N}(\mathbf{q})\mathbf{J}(\mathbf{q})$. The advantage of this transformation occurs when we use the invertibility of $\mathbf{N}$ to rewrite the linearization as,

[TABLE]

where $\mathbf{E}=\mathbf{N}^{-1}$ is our compliance matrix. Essentially we have factored our force into two components, one of which has an analytically invertible component. This allows us to isolate a poorly conditioned components from the system matrix, e.g.: the material parameters, which may include large or even infinite values.

9.2. Continuum Materials

To illustrate the above transformation we show how to apply it to the case of a strain energy density $\Psi(\mathbf{s})$, where $\mathbf{s}$ is some parameterization of the density function, e.g.: strains, strain rates, or principal stretches. Given such a function the elastic potential energy is

[TABLE]

where $V$ is the volume of integration, and the resulting force on the system is

[TABLE]

To obtain the compliance form of this force we factorize it as $\mathbf{c}=\partialderivative{U}{\mathbf{s}}^{T}$ and $\mathbf{J}=\partialderivative{\mathbf{s}}{\mathbf{q}}^{T}$. The derivative of $\mathbf{c}$ is then given by

[TABLE]

from which we can identify the compliance matrix

[TABLE]

In terms of matrix assembly, each element of a finite element discretization contributes $n_{s}=\text{dim}(\mathbf{s})$ constraints and Lagrange multiplier variables, and a $n_{s}\times n_{s}$ block to the system compliance matrix $\mathbf{C}$.

9.2.1. Linear Materials

For a constant strain element with general linear Hookean constitutive model we have the strain energy

[TABLE]

where $\mathbf{s}=[\epsilon_{xx},\epsilon_{yy},\epsilon_{zz},\epsilon_{yz},\epsilon_{xz},\epsilon_{xy}]$ is a vector of strain elements, $\mathbf{K}_{e}\in\mathbb{R}^{6\times 6}$ is the element’s material stiffness matrix, and $V_{e}$ the element rest volume. The corresponding compliance matrix is $\mathbf{E}=\left(V_{e}\mathbf{K}_{e}\right)^{-1}$ which is a constant that may be precomputed.

9.2.2. Hyperelastic Materials

Linear material models may exhibit significant volume loss during large deformations. Hyperelastic models, where stiffness increases as a function of strain, are less prone to this artifact, as illustrated in Figure 13. To incorporate isotropic hyperelastic materials we may parameterize the strain energy density in terms of principal stretches $\mathbf{s}(\mathbf{q})=[s_{1},s_{2},s_{3}]$ of the deformation gradient $\mathbf{F}(\mathbf{q})$. The corresponding material stiffness matrix $\mathbf{N}(\mathbf{q})\in\mathbb{R}^{3\times 3}$ is no longer constant, but is a function of the system configuration. During the Newton solve we evaluate $\mathbf{N}$ at each iteration and perform a direct inversion of it to obtain $\mathbf{E}$. A practical consideration for hyperelastic materials is that, as $\mathbf{E}$ is the inverse of a Hessian it may be indefinite, in which case it may be necessary to project it back to the positive definite cone before including it in our system matrices [Teran et al., 2005], although we have also found that a simple diagonal approximation to $\mathbf{E}$ is often sufficient. In Appendix B we give the derivation of $\mathbf{E}=(\partialderivative[2]{U}{\mathbf{s}})^{-1}$ for the stable Neo-Hookean model presented by Smith et al [2018].

9.3. Particles

For deformable bodies we use tetrahedral finite elements defined over a set of particles where a particle with index $i$ contributes 3 degrees of freedom,

[TABLE]

The kinematic mapping for particles is the identity transform, $\mathbf{G}_{i}=\mathbf{1}_{3x3}$ and the lumped mass matrix is $\tilde{\mathbf{M}}_{i}=m\mathbf{1}_{3x3}$, where $m$ is the particle mass.

9.4. Rigid Bodies

We describe the state of a rigid body with index $i$ using a maximal coordinate representation consisting of the position of the body’s center of mass, $\mathbf{x}_{i}$, its orientation expressed as a quaternion ${\bf\it\theta}_{i}=[\theta_{1},\theta_{2},\theta_{3},\theta_{4}]$, and a generalized velocity vector $\mathbf{u}_{i}$

[TABLE]

where ${\bf\it\omega}$ is the body’s angular velocity. The kinematic map from generalized velocities to the system coordinate time-derivatives is then given by $\dot{\mathbf{q}}_{i}=\mathbf{G}\mathbf{u}_{i}$ where

[TABLE]

where $\mathbf{Q}$ is the matrix that performs a truncated quaternion rotation,

[TABLE]

The corresponding mass matrix is then

[TABLE]

where $m$ is the body mass and $\mathbf{I}$ is the inertia tensor. When using quaternions to represent rotations there is an implicit constraint that $\|{\bf\it\theta}\|=1$, rather than include this directly in the Newton solve we periodically normalize the quaternion to unit-length.

10. Results

We implemented our algorithm in CUDA and run it on an NVIDIA GTX 1070 GPU. The assembly of the system matrix is performed on the GPU in compressed sparse row (CSR) format for maximum flexibility. We build $\mathbf{H},\mathbf{M},\mathbf{J}$, and $\mathbf{C}$ separately and perform the matrix multiply operations necessary for Krylov methods through successive multiplications of individual matrices.

Collision detection is performed between triangle-mesh features once per-time step using the system’s unconstrained velocity to generate candidate pairs. We define contact normals as the normalized vector between each feature pair’s closest points using the configuration at the start of the time-step. Constraint manifold refinement (CMR) could be used to further improve robustness [Otaduy et al., 2009].

10.1. Experimental Setup

In this section we discuss the test scenarios we have used for our method. We report scene statistics and performance numbers in Table 2. When running in an interactive setting we use a display rate of 60hz which corresponds to a frame time of 16.6ms. For robust collision detection we use two simulation time-steps per visual frame, each with $h=8.3$ms. If the simulation computation takes longer than this the effect for the user is a slightly slower than real-time update rate.

10.1.1. Fetch Tomato

We test our method on a pick-and-place task using the Fetch robot as shown in Figure 1. The robot consists of rigid bodies connected by joints as defined by the Unified Robot Description Format (URDF) file. The tomato is modeled using tetrahedral FEM with Young’s modulus of $Y=0.1$MPa, Poisson’s ratio of $\nu=0.45$, density of $\rho=1000\text{kg}/\text{m}^{3}$, and a coefficient of friction of $\mu=0.75$ between the grippers and the tomato. The robot is controlled by a human operator who directs the arm and the grippers to grasp the object and transfer it to the mechanical scales. The scales are modeled through rigid bodies connected by prismatic and revolute joints that drive the needle and accurately reads the weight of the tomato. The most challenging part of this scenario is the grasp and transfer of the tomato, which requires tight coupling between frictional contact and the internal dynamics of the tomato. Our method forms stable grasps while undergoing large translational and rotational motion.

10.1.2. Fetch Beam

We test our method on a flexible beam insertion task by modeling the beam as a series of connected rigid bodies with finite bending stiffness as shown in Figure 7. The lightweight beam is modeled by 16 rigid bodies each with mass of $m=0.003$kg, and connected through joints with a bending stiffness of $250\text{N.m}$. This example highlights the limitations of traditional relaxation-based approaches that cannot achieve the desired stiffness on the beam’s joints even with hundreds of iterations as shown in the convergence plot of Figure 18.

10.1.3. DexNet

We evaluate our method on the problem of dexterous grasping using the DexNet adversarial object database [Mahler et al., 2017]. These models are highly irregular with many concave areas that make forming stable grasps difficult. The underlying triangle meshes are high resolution and non-manifold which tends to generate many redundant contacts, a challenge for most complementarity solvers. We use the Allegro hand with a coefficient of friction $\mu=0.8$ to perform grasping using a human control interface, and find stable grasps for the objects in collection as shown in Figure 10.

10.1.4. PneuNet

To test coupling between deformable and rigid bodies we simulate a three-fingered gripper based on the PneuNet design [Ilievski et al., 2011]. We model the deformable finger using tetrahedral FEM with a linear isotropic material model and parameters for silicone rubber of $Y=0.01\text{GPa},\nu=0.47,\rho=1200\text{kg}/\text{m}^{3}$. To model inflation we use an activation function that uniformly adds an internal volumetric stress to each tetrahedra in the finger arches. We do not model the chamber cavity explicitly, however we found this simple activation model was sufficient to reproduce the characteristic curvature of the gripper. We observe robust coupling through contact by picking up a ball with mass $m=0.32\text{kg}$, using a friction coefficient of $\mu=0.7$ as shown in Figure 2.

10.1.5. Rigid Body Contact

We test our method on a variety of rigid-body contact problems as illustrated in Figure 14. Our method achieves stable configurations for difficult problems including self-supported structures, and stacks with extreme mass ratios. For the self-supported arch we use $\mu=0.6$ with masses in the range $m=[15,\cdots,110]$. For the heavy stack of boxes we use $\mu=0.5$, with masses that increase geometrically as $m=[8,64,\cdots,32768]$kg. For the table piling scene each table has a mass of $m=4.7$kg with $\mu=0.7$. Particularly on the scenes with high-mass ratios we observe that relaxation methods struggle to reduce error, while Krylov methods form stable structures and successfully prevent interpenetration.

10.1.6. Material Extension

We perform a material extension test and compare the behavior between a linear co-rotational model and the hyperelastic model of Smith et al. [2018] with Young’s modulus set to $E=10^{5}\text{Pa}$ and Poisson’s ratio of $\nu=0.45$. We visualize volumetric strain and observe high volume loss for the linear model as shown in Figure 13. In the supplementary video we show the impact of geometric stiffness on this test and find it is essential to obtain a stable simulation when strains are large. This is also illustrated in in Figure 12 which shows the iterate oscillating around the solution.

10.1.7. Reinforcement Learning

Reinforcement learning (RL) is a good test of robustness for a simulator since it generates many random inputs in the form of forces, torques, and constraint configurations. We apply our simulator to two scenarios using reinforcement learning. The first is the problem of training the ABB Yumi robot to grasp a cabinet handle and open a drawer as shown in Figure 15 (left). We use the PPO algorithm [Schulman et al., 2017] to train the network and find it quickly (less than 100 training iterations) learns a policy to extend, grasp and open the drawer through implicit PD controls applied through the joints. Our second RL example is the Humanoid Flagrun Harder scene adapted from the OpenAI Roboschool [2017]. In this task a humanoid model must learn to stand up and run in a randomly assigned direction that changes periodically as shown in Figure 15 (right). The learned actions are torques, applied as external forces. Using our simulator we are able to achieve good running results and note the agent taking advantage of stick-slip transitions on the feet during fast turns.

10.2. Effect of Complementarity Preconditiong

In Figure 16 we look at the distribution of convergence behavior over 132 simulation steps during the flexible beam insertion scenario. We use our proposed complementarity preconditioner with 40 PCR iterations per-Newton iteration and observe super-linear convergence in significantly more steps using our preconditioning strategy than with identity scaling. In addition, the median error with our strategy is often an order of magnitude lower for the same iteration count. We use the rigid body contact test scenarios to evaluate the effectiveness of our complementarity preconditioner, and compare the convergence of different strategies in Figure 18. We observe that identity scaling with $r=1$ may fail badly when there are large masses, or large contact sliding velocities. Constant global scaling by $r=h^{2}$ improves convergence, but still suffers from problems with large mass-ratios. The combination of time-step and effective mass leads to the most reliable observed convergence. Based on the poor performance of identity scaling we believe that a preconditioning strategy such as ours is a necessity to make such non-smooth formulations practical.

10.3. Effect of NCP-Function

We found that although both the minimum-map and Fischer-Burmeister functions can achieve high accuracy given enough iterations, the minimum-map tends to produce noisy results where the contact forces are not distributed evenly around a contact area. This is primarily a problem when the contact set is redundant, and a small change in the problem may lead to a large change in the active-set. We found the Fischer-Burmeister function was less sensitive to this problem, and would produce smooth contact force distributions even for redundant contact sets. We suspect this is because the minimum-map has many non-smooth points while Fischer-Burmeister has only one, however further investigation to verify this is needed. Due to the improved stability in interactive environments we have used the Fischer-Burmeister function for all examples unless otherwise stated. For scenarios where force distributions are critical, e.g.: force feedback based controllers, it may be appropriate to use a combination of warm-starting to provide temporal coherence, and redistribution of contact forces as a post-process after the contact solve, as proposed by Zheng & James [2011]. Yet another option is to change the contact model itself by introducing compliance. This makes the problem well-posed by allowing some interpenetration, and may be supported in our formulation by augmenting the compliance block on the contact constraints..

10.4. Effect of Linear Solver

In Figure 18 we compare convergence for different linear solvers over the course of a Newton solve during a single time-step. For performance sensitive applications it is typically not practical or desirable to run each Newton solve to convergence so for this test we fix the number of linear solver iterations per-step to 40. This early termination is generally no problem for relaxation and PCR methods, but it can cause problems for solvers like PCG which decrease the error non-monotonically, leading to erratic convergence. In Figure 11 we plot the behavior of each iterative method on a single linear subproblem. Please see the supplementary material for our test matrices and reference PCR implementation in MATLAB format.

The convergence of Jacobi and Gauss-Seidel with our contact formulation is consistent with the behavior observed in other engines such as Bullet [Coumans, 2015] and XPBD [Macklin et al., 2016]. While relaxation methods perform quite well for reasonably well-conditioned problems, they are very slow to converge for situations involving high mass-ratios as shown in Figure 18. This is highlighted in the heavy-stack example which shows catastrophic interpenetration.

10.5. Error Analysis

To better understand our results we perform an error analysis to establish a baseline accuracy limit given finite precision floating point. Our analysis is based on that given by Tisseur [2001] who shows that Newton’s method applied to the problem of solving $\mathbf{r}(\mathbf{x})=\mathbf{0}$ will have a limiting step-size, or solution accuracy of

[TABLE]

Here $\mathbf{x}_{*}$ is the true solution, $\mathbf{A}_{*}$ is the system Jacobian at the solution, $\nu$ is an upper bound on the residual error, and $\delta$ is the machine epsilon. In general we do not know the true solution $\mathbf{x}^{*}$ so we use the lowest error solution as an approximation. Likewise we use $\nu=\delta\|\mathbf{r}\|$ as the residual error bound, where $\mathbf{r}$ is the lowest achieved error. Likewise, the minimum residual is limited by available precision and Tisseur show that its predicted limiting magnitude is

[TABLE]

We use 32-bit floating point for all calculations, and plot these limiting values for the $L_{\infty}$ norm as dashed lines in Figures 18-18. We find that this error model does a good job of predicting the observed accuracy, with the exception of the flexible beam example where the predicted residual accuracy is under-estimated.

11. Limitations and Future Work

In robotics it is common to use a reduced coordinate representation to model rigid articulations. This establishes hard constraints directly in the system degrees of freedom, and reduces the number of parameters and explicit constraint equations. We have derived our formulation in terms of generalized coordinates, so our method should be applicable to any parameterization simply by using the appropriate mass matrix $\mathbf{M}$, kinematic map $\mathbf{G}$, and constraint Jacobians.

Krylov methods such as PCG and PCR use a globally optimal line-search step at each iteration. This has the sometimes undesirable effect that error in one row can affect the rate of residual reduction in other rows, leading to slower than necessary convergence for independent parts of the simulation. A natural solution would be to include an island-detection step and solve these independent systems separately.

We have not dealt specifically with elastic collisions, or energy preserving integrators, however we believe that the approach of Smith et al. [2012] can be combined with non-smooth complementarity formulations such as ours.

We found our geometric stiffness approximation to be quite effective on particle-based objects, but less so on rigid-articulated mechanisms, and in the worst case can cause some jitter at joint limits. To avoid this we apply geometric stiffness only to the particle degrees of freedom. In the future we plan to explore more advanced quasi-Newton methods e.g.: those based on symmetric rank-1 (SR1) updates.

Our method’s primary computational cost is the solution of a symmetric linear system at each step. All the iterative solvers considered here may be implemented in a matrix-free way which would likely improve performance. Because the linear solver is treated as a black box it would be interesting to apply more advanced methods, such as algebraic multi-grid (AMG), to accelerate convergence for larger-scale linear subproblems.

Finally, the design space for choosing the $r$ factor in the complementarity preconditioner is large, and we think that heuristics based on global information could provide significant performance improvements.

12. Conclusion

We have presented a framework for multi-body dynamics that allows off-the-shelf linear solvers to be used for rigid and deformable contact problems. We evaluate our method on a variety of scenarios in robotics and found that it performs well on grasping and dexterous manipulation of deformable objects, as well as ill-conditioned rigid body contact problems. We believe our complementarity preconditioner makes non-smooth formulations of contact practical for interactive applications for the first time, and hope that this opens the door for future works to apply new contact models, preconditioners, and linear solver methods to multi-body problems.

Acknowledgements.

The authors would like to thank Ken Goldberg and Jeff Mahler at Berkeley AUTOLAB for their help with the DexNet models. Many thanks to Jacky Liang for creating the Fetch flexible beam insertion environment and the Allegro hand interface. We also thank Yevgen Chebotar for contributing the RL Yumi scene, Jan Issac for creating the cabinet model, and Richard Tonge for creating the table model and piling scene. The Allegro hand URDF and models are provided courtesy SimLab Co., Ltd. The Fetch robot URDF and model is provided courtesy of Fetch Robotics, Inc. The Yumi robot URDF is provided courtesy of ABB robotics. The humanoid model is adapted from the MuJoCo model created by Vikash Kumar.

Appendix A Contact Formulation

In this appendix we give the full form of the contact constraints and derivatives, including the $r$-factor scaling parameter, for the minimum-map and Fischer-Burmeister NCP-functions.

A.1. Contact

The contact constraint given in terms of the Fischer-Burmeister function can be written as,

[TABLE]

which has the following derivatives,

[TABLE]

Likewise, for the minimum-map,

[TABLE]

with derivatives given by,

[TABLE]

A.2. Friction

We now give the closed form of the frictional compliance terms including $r$-factor scaling. First, we write in full the Fischer-Burmeister NCP-function for the frictional conditions (32),

[TABLE]

We use the equation ${\psi_{f}}_{\text{FB}}=0$ to obtain a fixed point iteration for the quantities $|\mathbf{D}^{T}\dot{\mathbf{q}}|$ and $|{\bf\it\lambda}_{f}|$ as follows,

[TABLE]

[TABLE]

Using these we can write our frictional compliance quantity $W_{FB}$,

[TABLE]

Following the same pattern we find that $W_{\text{min}}$ for the minimum-map is given by,

[TABLE]

Using the $W$ for one or the other NCP-function, the compliance matrix sub-block for a single contact with index $i$ is

[TABLE]

Appendix B Compliance Form of Stable Neo-Hookean Materials

Here we give the derivatives required for the compliance form of the stable Neo-Hookean material model introduced by Smith et al [2018]. The elastic strain-energy density is given by

[TABLE]

The material constants $C_{1}$, $D_{1}$, and $\alpha$ can be related to the Lamé parameters according to the original paper. Defining $\mathbf{s}=[s_{1},s_{2},s_{3}]^{T}$ as the vector of principal stretches of the deformation gradient $\mathbf{F}$ we have $I_{C}=s_{1}^{2}+s_{2}^{2}+s_{3}^{2}$ the first invariant of stretch, and $J=s_{1}s_{2}s_{3}$, the relative volume change induced by $\mathbf{F}$. To perform the compliance transformation we require the Jacobian and Hessian of $\Psi$ with respect to $\mathbf{s}$, which we provide below:

[TABLE]

For a constant strain tetrahedron with rest volume $V_{e}$ the elastic potential energy is $U(\mathbf{q})=V_{e}\Psi_{E}(\mathbf{s})$, and the compliance matrix is $\mathbf{E}=\left(V_{e}\partialderivative[2]{\Psi_{E}}{\mathbf{s}}\right)^{-1}$ which can be obtained through a $3\times 3$ matrix inverse. The derivatives of the singular values with respect to the vertices, $\mathbf{J}=\partialderivative{\mathbf{s}}{\mathbf{q}}$, are given in [Perez et al., 2013].