Andrews, Sheldon

Erleben, Kenny

Teichmann, Marek

English

Copenhagen University Research Information System

u n i ve r s i t y  o f  co pe n h ag e n  Constraint reordering for iterative multi-body simulation with contactAndrews, Sheldon; Erleben, Kenny; Teichmann, MarekPublication date:2017Document versionPeer reviewed versionDocument license:UnspecifiedCitation for published version (APA):Andrews, S., Erleben, K., & Teichmann, M. (2017). Constraint reordering for iterative multi-body simulation withcontact. Abstract from ECCOMAS Thematic Conference on Multibody Dynamic, Prag, Czech Republic.Download date: 09. okt.. 2020ECCOMAS Thematic Conference on Multibody DynamicsJune 19-22, 2017, Prague, Czech RepublicConstraint reordering for iterative multi-body simulation with contactSheldon Andrews1,2, Kenny Erleben3, Paul G. Kry1, Marek Teichmann21 McGill UniversityMontreal, Canadasheldon.andrews@mail.mcgill.cakry@cs.mcgill.ca2 CM Labs Simulations, Inc.Montreal, Canadamarek@cm-labs.com3 University of CopenhagenCopenhagen, DKkenny@di.ku.dkAbstractMulti-body simulations with contact are non-smooth systems and wrought with discontinuities which arisedue to non-interpenetration and frictional constraints. Linear systems are used for applications where real-timeperformance is a concern, such as interactive training or video games, which gives rise to a linear complementarityproblem (LCP). A common mathematical formulation [5] of the LCP for the velocity-level equations of motion isJM−1JT︸ ︷︷ ︸A∆tλ︸︷︷︸z+J(v+∆tM−1fext)︸ ︷︷ ︸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 thegeneralized 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 tangentialfrictional 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 pivotingapproach by Judice and Pires [6]. These methods are able to provide exact solutions to the LCP, but are com-putationally infeasible for more than several hundred contact constraints. Iterative methods are more prolific forsimulations involving a large number of contacts since an approximate solution can be found after only small num-ber of iterations. Where performance is a concern, the algorithm can terminate early once specified error tolerancehas 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. No-tably, the projected Gauss-Seidel (PGS) method [1, 2, 5], which has even been extended to handle non-linearcomplementarity problems [7]. These solvers are used extensively in computer graphics applications due to theirspeed, stability, and convergence properties. Briefly, the PGS algorithm works by splitting the lead matrix asA= 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 thatzk+1i ←bi−∑i−1j=1Li, jzk+1j −∑mj=i+1Ui, jzkjDi,i, (2)which is followed by a projection step zk+1i = min(max(zlo,zk+1i ),zhi). It’s obvious by inspection of Eq.(2) thatthe solution at each fixed-point iteration is dependent on the constraint order. This affects the number of iterationsrequired 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 physicalbehavior. Likewise, Erleben [5] observes that correct stacking and shock propagation behavior cannot be simulatedusing standard PGS or Jacobi solvers and proposes a modified algorithm that sequentially solves blocks of contactconstraints spatial ordering (e.g., in a bottom up order for stacking).Figure 1: The examples used in our experiments. A capsule (left) andthree sphere body (right) slide across a planar surface with an instanta-neous angular velocity 1 rad/s and linear velocity of 10 m/s with directionsshown as dashed orange arrows. The capsule and three sphere body aresimulated 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 af-fects rate of convergence. We highlight this by simulatingthe examples shown in Fig. 1 using the complimentarity for-mulation of [5]. Since the simulations involve only two orthree frictional contacts, all possible permutations of con-straint equation ordering can easily be evaluated (i.e. thereare 720 and 362880 permutations, respectively, for the cap-sule and three sphere body example). The default constraintordering is such that if j = (i mod 3) is zero, it correspondsto the non-interpenetration constraint of the jth contact, androws 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 con-verges in less than 20 iterations for both examples. However, for the worst ordering, the error remains high even0 5 10 15 20 25Iteration10-1010-5100LCP errorDefaultBestWorst0 5 10 15 20 25Iteration10-1010-5100LCP errorDefaultBestWorstFigure 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 theconstraint equations (black) is compared versus the optimal order (blue) and worst case order (red).0 5 10 15 20 25Iterations00.050.10.150.20.250.3% permutationsIterations to reach error < 1e-40 5 10 15 20 25Iterations00.020.040.060.080.10.120.140.160.18% permutationsIterations to reach error < 1e-4Figure 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 requiredto reach a reasonable error threshold.Motivated by these results, our work investigates strategies to accelerate the convergence of iterative solversfor multibody simulation by reordering of the constraint equations. We present an analysis of the following strate-gies:• 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 frictionalcontact and develop heuristics that allow automatic re-ordering, and grouping, of constraint equations to improvesolver performance.References[1] Mangasarian, OL. (1977). Solution of symmetric linear complementarity problems by iterative methods.Journal of Optimization Theory and Applications, 22(4), p. 465.[2] Ko, M., and Zowe, J. (1994). An iterative two-step algorithm for linear complementarity problems. Nu-merische Mathematik, 68(1), p. 95.[3] Fratarcangeli, M., Tibaldo, V., and Pellacini, F. (2016). Vivace: a practical Gauss-Seidel method for stablesoft body dynamics. ACM Transactions on Graphics (TOG), 35(6), p. 214.[4] Stam, J. (2009). Nucleus: Towards a unified dynamics solver for computer graphics. In 11th IEEE Interna-tional Conference on Computer-Aided Design and Computer Graphics, (p. 1–11).[5] Erleben, K. (2007). Velocity-based shock propagation for multibody dynamics animation. ACM Transac-tions on Graphics (TOG), 26(2), p. 12.[6] Júdice, J. J., and Pires, F. M. (1994). A block principal pivoting algorithm for large-scale strictly monotonelinear complementarity problems. Computers and Operations Research, 21(5), p. 587-596.[7] Tasora, A., and Anitescu, M. (2009). A fast NCP solver for large rigid-body problems with contacts, friction,and joints. In Multibody Dynamics, pp. 45–55.

Constraint reordering for iterative multi-body simulation with contact

https://curis.ku.dk/portal/files/180633715/Andrews_2017_Constraint_recordering.pdf

Constraint reordering for iterative multi-body simulation with contact

Abstract

Similar works

Full text

Available Versions

Copenhagen University Research Information System