Constraint reordering for iterative multi-body simulation with contact

Multi-body simulations with contact are non-smooth systems and wrought with discontinuities which arise due to non-interpenetration and frictional constraints. Linear systems are used for applications where real-time performance is a concern, such as interactive training or video games, which gives rise to a linear complementarity problem (LCP). A common mathematical formulation [5] of the LCP for the velocity-level equations of motion is JM−1JT } {{ } A ∆tλ }{{} z +J(v+∆tMfext) } {{ } b = w (1) zlo ≤ z≤ zhi, w⊥ z, where J ∈ Rm×n is the Jacobian matrix encoding the non-penetration and friction constraints, M ∈ Rn×n is the generalized mass matrix, v and fext ∈ Rn are the generalized velocities and external forces of simulation bodies, respectively, λ ∈ Rm are Lagrange multipliers representing the non-interpenetration normal forces and tangential frictional forces of each contact. The box constraints defined by zlo and zhi contain the lower and upper bound, respectively, of the normal and frictional impulses. The unilateral and discontinuous nature of the system in Eq.(1) is problematic for many numerical solvers. Previous work has solved the LCP using simplex based pivoting methods, such as Lemke’s or the block pivoting approach by Judice and Pires [6]. These methods are able to provide exact solutions to the LCP, but are computationally infeasible for more than several hundred contact constraints. Iterative methods are more prolific for simulations involving a large number of contacts since an approximate solution can be found after only small number of iterations. Where performance is a concern, the algorithm can terminate early once specified error tolerance has been reached, or computational time budget has been exceeded. Gauss-Seidel and Jacobi based solvers are among the most popular iterative methods for solving LCPs. Notably, the projected Gauss-Seidel (PGS) method [1, 2, 5], which has even been extended to handle non-linear complementarity problems [7]. These solvers are used extensively in computer graphics applications due to their speed, stability, and convergence properties. Briefly, the PGS algorithm works by splitting the lead matrix as A = L+D+U, where D, L, and U are the diagonal, strictly lower triangle, and strictly upper triangle parts of A, respectively. The algorithm loops over each variable i ∈ m and updates its value at iteration k+1, such that zk+1 i ← bi−∑ j=1 Li, jz k+1 j −∑ m j=i+1 Ui, jzj Di,i , (2) which is followed by a projection step zk+1 i = min(max(zlo,z k+1 i ),zhi). It’s obvious by inspection of Eq.(2) that the solution at each fixed-point iteration is dependent on the constraint order. This affects the number of iterations required to find a solution to the LCP, which is important if early termination is required. Furthermore, Stam [4] notes that for iterative solvers, the order in which constraints are solved has an effect on the resulting physical behavior. Likewise, Erleben [5] observes that correct stacking and shock propagation behavior cannot be simulated using standard PGS or Jacobi solvers and proposes a modified algorithm that sequentially solves blocks of contact constraints spatial ordering (e.g., in a bottom up order for stacking). Figure 1: The examples used in our experiments. A capsule (left) and three sphere body (right) slide across a planar surface with an instantaneous angular velocity 1 rad/s and linear velocity of 10 m/s with directions shown as dashed orange arrows. The capsule and three sphere body are simulated using two and three frictional contact constraints, respectively, and coefficient of friction μ = 1.0. Each body has a mass of 10 kg. Our own experiments verify that the constraint order affects rate of convergence. We highlight this by simulating the examples shown in Fig. 1 using the complimentarity formulation of [5]. Since the simulations involve only two or three frictional contacts, all possible permutations of constraint equation ordering can easily be evaluated (i.e. there are 720 and 362880 permutations, respectively, for the capsule and three sphere body example). The default constraint ordering is such that if j = (i mod 3) is zero, it corresponds to the non-interpenetration constraint of the jth contact, and rows i+1 and i+2 are the associated friction equations. The convergence plots for the default, best, and worst orderings are show in Fig. 2. The best ordering converges in less than 20 iterations for both examples. However, for the worst ordering, the error remains high even 0 5 10 15 20 25 Iteration 10-10 10-5 100 LC P e rr or Default Best Worst 0 5 10 15 20 25 Iteration 10-10 10-5 100 LC P e rr or Default Best Worst Figure 2: The convergence plots for the sliding capsule example (left) and three sphere body (right). The error at each iteration using the default order of the constraint equations (black) is compared versus the optimal order (blue) and worst case order (red). 0 5 10 15 20 25 Iterations 0 0.05 0.1 0.15 0.2 0.25 0.3 % p er m ut at io ns Iterations to reach error < 1e-4 0 5 10 15 20 25 Iterations 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 % p er m ut at io ns Iterations to reach error < 1e-4 Figure 3: Histograms showing the number of iterations required to reach an error < 10−4 for the sliding capsule (left) and three sphere body (right). after 25 iterations. Furthermore, as indicated by Fig. 3, there is a large variation in the number of iterations required to reach a reasonable error threshold. Motivated by these results, our work investigates strategies to accelerate the convergence of iterative solvers for multibody simulation by reordering of the constraint equations. We present an analysis of the following strategies: • Solving constraint equations in a randomized order; • Re-ordering constraint equations by heuristics based on the complementarity error and the effective mass; • Grouping constraint equations and solving for several variables at once by a blockwise PGS. We investigate the viability of each strategy for a number of rigid body simulation scenarios involving frictional contact and develop heuristics that allow automatic re-ordering, and grouping, of constraint equations to improve solver performance.