{'fraction_non_alphanumeric': 0.028989258444415685, 'fraction_numerical': 0.009317975928562185, 'mean_word_length': 4.614993338985103, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving variational inequalities and cone complementarity problems in nonsmooth dynamics using the alternating direction method of multipliers

This work presents a numerical method for the solution of variational inequalities arising in nonsmooth flexible multibody problems that involve set‐valued forces. For the special case of hard frictional contacts, the method solves a second order cone complementarity problem. We ground our algorithm on the Alternating Direction Method of Multipliers (ADMM), an efficient and robust optimization method that draws on few computational primitives. In order to improve computational performance, we reformulated the original ADMM scheme in order to exploit the sparsity of constraint jacobians and we added optimizations such as warm starting and adaptive step scaling. The proposed method can be used in scenarios that pose major difficulties to other methods available in literature for complementarity in contact dynamics, namely when using very stiff finite elements and when simulating articulated mechanisms with odd mass ratios. The method can have applications in the fields of robotics, vehicle dynamics, virtual reality, and multiphysics simulation in general.

they also generalize to the case where velocity is assumed to be a (possibly) discontinuous function of bounded variation, therefore they admit impulsive events. 3 MDIs can be solved using special time stepping methods that offer high stability and robustness, however this comes at the cost of solving a Complementarity Problem (CP) per each time step. 4 These CPs correspond to the broader class of Variational Inequalities (VI). 5,6 In general, the solution of VIs or CPs represents the major numerical bottleneck in the entire time stepping process of the MDI, especially when dealing with many parts or many contacts. For instance, similar scenarios can happen when simulating packaging devices, ground-machine interaction, rock dynamics, masonry structures, and granular materials, where the unknowns of the VIs (contact forces) can amount to thousands or millions. This challenge stimulated research on efficient numerical methods in the last three decades.
Among the former strategies for solving this class of problems, we cite the idea of solving a Linear Complementarity Problem (LCP) per each time step. 7 Once cast as a LCP, the VI problem can be solved by direct LCP solvers such as the Lemke algorithm. However, LCP direct solvers often lead to complex subiterations that do not scale well with increasing number of unknowns, hence they are not much used nowadays. Another drawback of casting the problem as a LCP is that they introduce a polyhedral approximation, for example the Coulomb friction cone would be approximated as a faceted pyramid, introducing a numerical source of anisotropy.
At the other end of the spectrum are iterative methods that aim at high robustness and efficiency even at the expense of sacrificing accuracy, since this is often acceptable in fields like real-time simulation, video games and virtual reality. 8 In these fields, the most used approach nowadays consists in fixed point iterations, similarly to Gauss-Seidel, SOR or Jacobi stationary methods, interleaved with simple projections onto the constraint sets at each iteration. [9][10][11] The major issue of these methods is that they exhibit slow convergence when the system contains long kinematic chains or long sequences of contacts, such as when including high stacks of objects. 12 Their convergence is even more challenged if odd mass ratios are encountered: this often happens when simulating industrial robots, earth-moving machines, complex vehicle suspensions, car drivelines. In these cases the convergence might stall soon, hence if the iterations are truncated in order to meet a time budget in real-time applications, the resulting simulation can be affected by low accuracy-usually detectable in the form of objects that interpenetrate at the point of contact and mechanisms that fall apart. Optimizations like warm starting and substepping can alleviate the drawbacks of these methods. 13 A more advanced class of solvers is represented by non-smooth Newton methods, such as the one presented in Reference 14, that assume a generic Nonlinear Complementarity Problem (NCP), again a subcase of a VI, and that enforce the NCP complementarity constraints using nonsmooth functions like the Fisher-Burmeister function. In these methods, a generalized nonsmooth Newton method can be used to find the zero of the functions, but this comes at the cost of solving a linear system per each iteration.
Recent efforts pushed forward the idea of casting the VI of nonsmooth dynamics as a convex optimization problem, namely a Quadratic Program (QP) with convex constraints. This requires a simplifying assumption about the Coulomb friction being associative instead of nonassociative. This artificial associativity of the friction model means that the relative sliding velocity in a contact point will be restricted to the dual cone of the Coulomb friction cone, whereas the original model has a more relaxed condition, assuming that such velocity could be even parallel to the surface: the practical consequence is that grazing motion cause an artificial lift-off-however, this artifact can be attenuated thanks to the introduction of a stabilization term. 15 In this setting, the problem is also a convex Cone Complementarity Problem (CCP). The convexity assumption allows the adoption of numerical methods already available for large-scale optimization problems. In Reference 16, the Barzilai-Borwein spectral projected gradient is used to this end, and in Reference 17 we used the Nesterov accelerated projected gradient descend to simulate large-scale contact problems. These are first-order optimization methods, in the sense that they do not require the computation of the Hessian matrix: in fact they rely only on matrix-by-vector multiplications and inner products (similarly to Krylov subspace methods in the field of linear systems), plus a fast projection operator. Although they offer superior convergence with respect to fixed point iterations, they are not easily applicable to problems involving finite elements. In fact, this class of solvers often perform the matrix-by-vector primitive over a Schur complement matrix, a product of constraint jacobians and inverse of the mass matrix. Since the mass matrix is often (block) diagonal, the Schur complement can be quickly factored at the beginning of the iteration, but if finite elements are used, the Schur matrix would also involve to some extent the inverse of element stiffness matrices and damping matrices, a costly and inefficient operation unless special developments are considered. 18 Extending complementarity approaches in order to encompass finite elements is important because recent research has shown that complementarity-based formulations for large deformation dynamics with contacts often show superior performance when compared to penalty-based approaches. 19,20 Another class of solvers that can draw on the optimization-based approach are Interior Point Methods (IPMs). These solvers share some similarities with the previously mentioned nonsmooth Newton methods, in that they require the solution of a saddle-point linear system at each iteration. IPMs offer the best theoretical convergence for the problem at hand. 21 However, their implementation is complex, and despite the encouraging theoretical properties, in practical scenarios they do not scale well for large problems. One of their drawbacks, also, is that at the present state of knowledge there are no efficient ways to warm-start them. 22 This led us to investigate other methods for constrained optimization problems.
The Alternating Direction Method of Multipliers (ADMM) has been proposed in the 1970s 23,24 as a practical way to solve constrained optimization problems by iterating over the minimization of two (not necessarily differentiable) objectives. ADMM introduces auxiliary variables and casts the original problem as the minimization of a separable objective subject to a linear constraint between the original variables and the auxiliary variables. With a proper choice of auxiliary variables, the minimization of the separable function is obtained iterating on two distinct minimization steps of lower complexity. The method iterates updating primal variables, auxiliary variables and dual variables similarly to a fixed point iteration: as such, they compare less favourably to IPMs if one is interested in convergence to high precision; however, recent developments showed that in practical scenarios they offer superior speed, scalability, and robustness. 25 These positive properties triggered a recent revival of ADMM and similar operator-splitting methods even as alternatives to IPMs, especially in distributed convex optimization with large-scale problems where moderate precision is acceptable. 26 ADMM methods have been used in many fields, with recent and notable examples in computer graphics, 27 computational fluid dynamics, 28 simulation of deformable structures, 29 and contact dynamics. 30,31 In the field of robotics, a sequential unconstrained minimization technique (SUMT) method based on augmented Lagrangian has been discussed in References 32,33, to our knowledge for the first time in a context of nonsmooth dynamics: SUMT leads to an iteration with the same structure of ADMM, for the same type of constraints.
Recent developments addressed the possibility of accelerating the ADMM method. In Reference 34 an acceleration scheme has been proposed for special cases of symmetric ADMMs, and assuming that one of the two objectives is differentiable. In Reference 35 the Nesterov acceleration method has been proposed for problems where at least one of the two objectives is quadratic, and under the assumption that both objectives are strongly convex; if the latter assumption is not guaranteed, restarting strategies are needed. In Reference 36 a similar acceleration method has been proposed, convergence rate, but without the assumption that one of the objectives must be quadratic. Other types of improvements, based on the Anderson fixed point acceleration, have been put forward in References 27,37, showing high performance in computer graphics applications.
Recent developments also focused on the possibility of parallelizing the ADMM iterations on GPU and parallel architectures. 38 In Reference 39 a specialized ADMM method for quadratic problems has been presented: building off this method, we develop a solver that can exploit the nature of nonsmooth dynamical problems. An attractive property of the approach is that it requires the solution of a linear system per each iteration, but unlike IPMs, the linear system most often remains unchanged during the iterations, so a factorization can be reused multiple times with a great benefit in terms of speed.
In the approach proposed herein, unknowns are primal variables (velocity measures) and dual variables (impulses in contact points and joints), the convex cone constraints stem from the Coulomb friction law, and the system matrix includes terms from the mass matrices and from the tangent stiffness/damping matrices (hence it is sparse, but not necessarily diagonal). This paper will present a model for the non-smooth dynamics, it will show how to cast it as a QP with convex conic constraint, it will present the ADMM method to solve it, it will discuss practical implementation details and then it will show benchmarks.

DEFINITIONS
In this section we list some basic definitions that we will use in the rest of the paper.

Definition 1.
A second-order cone, also Lorentz cone, is a self-dual self-scaled symmetric cone defined as Definition 2. The dual cone  * of a set  in a vector space equipped with an inner product, is always convex and is defined as Definition 3. The polar cone is defined as Definition 4. The normal cone to a closed convex set  at the point x ∈  is closed and convex and is defined as Definition 5. The indicator function of a subset  ∈  is a scalar function I ∶   → R defined as: Remark. The indicator function of any closed set is lower semicontinuous.

Definition 6.
The subdifferential f (x 0 ) at x 0 of a convex, possibly nondifferentiable scalar function f ∶ R n  → R, is the closed convex set of all subgradients g at x 0 : Remark. The subdifferential is a set-valued function f (x) ∶ R n ⇉ R n in general, but as a special case, Remark. The subdifferential of an indicator function of a convex set is also its normal cone: Definition 7. The projection of y on a nonempty closed convex set  is a R n → R n mapping defined as: Remark. If f (y) is the indicator function I  (y) to a nonempty closed convex set , defined as in (5), then prox f (y) is the orthogonal projection onto that set:

THE NONSMOOTH MULTIBODY MODEL
Different time stepping schemes have been proposed for MDIs, here we refer to the one discussed in Reference 12 without lack of generality. After regularization and convexification and discretization, the MDIs leads to a major numerical bottleneck to be solved at each time step: a (mixed) CCP with unknowns v and : where • unknown v ∈ R n v is the speed at the end of the time step h, • unknown ∈ R n c is the reaction in contacts and bilateral joints, a vector-signed measure with Lebesgue decomposition in atomic parts (impacts) and continuous parts (continuous reactions), • H ∈ R n v ×n v is a positive definite matrix containing M, the block-diagonal matrix of the masses and inertia tensors of the bodies; if stiff loads are added too (ex. when finite elements come into play) it becomes a block-sparse matrix including the tangent stiffness and damping matrices, for instance H = M − h 2 ∇ q f − h∇ v f in case of the first-order time steppers in References 12,18, • D ∈ R n v ×n c is a sparse matrix, the transpose jacobian of all constraints, • k ∈ R n v is a vector containing terms proportional to applied forces f and to the last-known speed, • b ∈ R n c is a vector containing constraint stabilization terms and rheonomic terms, • Υ is the Cartesian product of all cones of admissible constraint forces, Υ= ⨉ j Υ j , in detail: • if a frictional contact is added, Υ j ⊂ R 3 is a second-order Lorentz cone (see Definition 1) with aperture proportional to friction coefficient, For more details on the model we refer to Reference 40. Since H is positive definite, from (10a) one also gets: and Equation (10b) can be written as where we introduced the Schur complement * and the vector In this way, one can also reformulate the problem as the following CCP: This CCP corresponds exactly to a first-order optimality condition of a convex quadratic program: and in fact, the CCP (14) can be written in the more conventional language of the Karush-Kuhn-Tucker optimality conditions on dual variables y multipliers, for y = −u , and primal variables : After the convex program (15) is solved by ADMM, one can compute v using (11) with an immediate step. The ADMM method in Reference 39 can be used to solve problems in the form * The N matrix, which is the Schur complement of the matrix of the linear system Equations (10a) and (10b), is also called Delassus operator in the contact dynamics community. and recalling Definitions (3) and (4), one can see that (16) is a special case of (17) where A = I, P = N,  = Υ, q = r, x = , b = 0, where some optimizations can take place because of the structure of our problem. We remark that the ADMM method just makes the assumption of Υ being convex, so the problem can be generalized to ∈  where  is a generic convex set, and at the same time we assume an associated † flow u ∈ −  ( ). For instance,  could be a capped friction cone to represent plasticization of contacts, or a cone translated downward to represent possible adhesion in contact up to a threshold, or the Von Mises yield region if represent stresses in finite elements undergoing plasticization.
Also, aiming at highest generality, in presence of finite elements one might need to use tangent stiffness matrices K and damping matrices D to accommodate an implicit integration scheme for stiff elements, and this means that the Schur complement would be computed as N = D T H −1 D . Here H is a linear combination of M, K, D (depending on the integration scheme) and in general is not block-diagonal anymore.

THE ADMM SOLVER
For variables (x, z) ∈ R m × R n , the ADMM methods solve constrained optimization problems with separable structure where f ∶ R n x → R and g ∶ R n z → R are proper convex lower-semi-continuous functions, and A ∈ R n x ×n y , B ∈ R n z ×n y are linear operators. The optimality conditions require that primal feasibility and dual feasibility must hold at the solution (x ⋆ , z ⋆ , y ⋆ ), saddle point of the Lagrangian with dual variables y ∈ R n y . By introducing an augmented Lagrangian the ADMM method converges to the solution (x ⋆ , z ⋆ , y ⋆ ) by iterating over two distinct minimization problems as in the following loop: The convergence of such method has rate  . Although the convergence is fast in the first iterations, it tends to deteriorate later as in most fixed-point iterations; however in practical scenarios where loose tolerances can be accepted, even this basic formulation of ADMM proves to be efficient and robust.
The performance of ADMM depends on the efficiency of the updates in (22) and (23). We will design the method such that the step (23) will correspond to a very efficient projection onto a cone with separable structure, † Not all dissipative constitutive models are associated. For example, in some computational plasticity models, the plastic flow might deviate with respect to the normal to the yield surface. Moreover, the Coulomb contact model itself is not associated-in fact the contact velocity u can range into a wider cone than the polar of the friction cone, but our MDI model compensates for such effect via a stabilization term. We will assume associated models heretofore, and we do not deal with non-associated models to avoid additional complexity for the moment.
whereas the only bottleneck will be (22), corresponding to a quadratic optimization to be solved with a linear system.
In order to rewrite (15) as a sum of two functions as in (20), we introduce ∈ R n c as primal variables, z ∈ R n c as auxiliary variables, we add − z = 0 as the linear constraint, we reformulate the ∈ Υ constraint by adding a nonsmooth penalty function I Υ (z) given by an indicator function (see Definition 5), and finally we have We can write the augmented Lagrangian introducing a step size parameter and a vector of dual variables y, obtaining: that is also, using the property y T r The first step of ADMM requires to compute = arg min = arg min This is an unconstrained convex quadratic program whose optimality conditions lead to the following linear problem: The linear system (31) can be solved as it is, but in our case it would be better to exploit the fact that the Schur matrix N is a product D T H −1 D . We would like to avoid computing H −1 , even storing a precomputed inverse for all iterations would be unpractical because often dense (except when H is a diagonal mass matrix M). So we propose to replace (31) with an equivalent saddle point problem:

Theorem 1. Let k+1 be a solution to (31), then it is also a solution to
[ Proof. Performing the multiplication of the upper part of the saddle point matrix one has Hv k+1 − D k+1 = k, pre-multiplying by Recalling the definition of N in (12) and the definition of r in (13), this can be written also as (N + I) k+1 = z k − r − y k . ▪ Corollary 1. The linear system (32) introduces an auxiliary variable v ∈ R n v , however, the matrix is very sparse and does not require storing or computing any H −1 matrix as we should do if using (31).

Corollary 2.
A positive side effect of this approach is that its auxiliary variable v is also the velocity term in the original mixed CCP (10), so one does not need to compute it as v = H −1 (k + D ) after the iteration converged because it is already a byproduct of the linear solver in (32).

Corollary 3.
If constraint compliance is added to the MDI, a compliance matrix C can be present in the CCP, thus leading to a modified version of (10b): In such case one would have N = D T H −1 D + C , and (32) would change into: The second step of ADMM requires to compute This is a minimization problem too, but it can be rephrased in terms of a projection on the Υ set, which has a separable structure. In fact, recalling (27) and removing constant terms: Using the Definition 7, this leads to: Note that the projection Π Υ can be performed efficiently in our case, because Υ is the Cartesian product of n c Coulomb second-order friction cones of lower dimension, assuming a set of n c contacts: The separable structure allows the projection to be performed as a tuple of simpler projections: In a more general setting, where contains both contact reactions and reactions in bilateral and unilateral joints, the projection operator for the ith bilateral reaction is simply Π R (p i ) = p i and the projection operator for the ith unilateral reaction is Π R + (p i ) = max(0, p i ).
Note that the entire Π Υ projection requires an inexpensive and parallelizable operation. A basic form of ADMM will then iterate over these updates: Before being able to apply this iteration in real problems, we should introduce some optimizations that will make a big difference for the performance of the method, and that will be discussed in the next section.

Residuals and termination
In order to monitor the convergence of the method, a measure of residual quantities must be introduced. Given the definitions (19), (20a), (20b), one can define a primal and a dual residual as: An alternative way to compute r k+1 dual , that does not need to use the N matrix, is Another option (see Reference 25) is to compute it as which is even faster, at the expense of storing the previous value of z.
The iteration is terminated when both r primal and r dual fall under prescribed tolerances.
Here primal and dual residuals have an interesting physical interpretation: a large primal residual means that the solution provides reaction forces that do not satisfy the set inclusions (e.g., the friction might be overestimated), whereas a large dual residual means that the speed (in the contact point metric) is wrong. That is, in the discussed time stepping scheme, the measuring units of r primal is the same of reaction impulses, e.g. N.s, whereas the measuring units of r primal can be considered the same of speeds, for example, m/s.
We remark that, depending on the measuring units, the two residual errors can have different importance: this means that it is useful to set two distinct tolerances primal and dual . Different heuristics can be used to define relative tolerances too, as shown in Reference 39.

5.2
Step size selection The basic method outlined in (36)-(38) uses a fixed step .
In general, small values of cause a fast reduction of the dual residual, but a slow reduction of the primal residual. This has a physical interpretation: truncating iterations based on too small most often lead to solutions where parts do not interpenetrate, but some contact forces do not satisfy the inclusion in the Coulomb cone. That is, some parts might stick like contacts were glued.
On the other hand, large values of cause a fast reduction of the primal residual, but a slow reduction of the dual residual. This means that truncating iterations with too large most often lead to solutions where contact forces satisfy very well the inclusion in the Coulomb friction cones, but parts might interpenetrate.
Clearly, this calls for the introduction of a scheme that adaptively changes the step size in order to achieve a balanced reduction of both residuals. In most cases one can use heuristics and adjust only one or two times during the iterations, hence allowing to reuse the factorization of the linear system (36) as much as possible. In fact, the most time-consuming step is the factorization of the matrix in (36), but this happens only at the first iteration and at all iterations where changes, otherwise one can reuse the last factorization because the matrix is the same, and perform only the back solve step with the different right-hand side. On average, in our tests, the back solve operation took less than 1/10 of the time needed for the factorization.
We tested different policies for scaling the step in our algorithm. All policies share the following concept: starting from an initial value 0 , the step is scaled once each n s iterations such that the primal and dual residuals are well balanced, 41 that is: We define the balanced scaling policy as We define the balanced normalized scaling policy as Another idea is presented in Reference 42, where the the step is changed using the Barzilai-Borwein spectral formula. This is referenced as the spectral scaling policy in the following.
In all the cases above, one needs safeguards to avoid excessive scaling, for example we limit the scaling factor by clamping it in the ] interval by default. Also, in order to avoid an excessive number of factorizations, one could force the scaling factor to = 1 (hence no change in ) if it is not too far from the unity anyway-for example, by default we use a tolerance interval , and if the scaling falls inside that interval, it is forced to = 1.

Custom treatment of bilateral joints
We experimented with the idea of using a value that changes on a per-constraint basis. If one could know in advance which contact will be active at the solution, one could split i values between very high or very low values for the fastest convergence; however, this is not possible because the set of active constraints is known only at the solution. However, we still can use a simple heuristics consisting in forcing a small value b = 1 × 10 −9 for constraints belonging to the set  B of bilateral constraints, where one knows in advance that the corresponding dual variable y i would be null anyway. This produces better convergence results in systems containing both contacts and many bilateral joints, such as when simulating articulated robots that grasp some objects. This means that, instead of using a scalar k , in our code we use a diagonal matrix Θ k defined as:

FINAL ALGORITHM
Using the above results, one can finally obtain the ADDM algorithm for solving the CCP of the nonsmooth dynamics timestepper: The algorithm above contains a single computational bottleneck, that is the solution of the linear system. However, if a direct solver is used, a LU decomposition ‡ could be computed at the beginning and then recomputed only when the step size is adjusted, otherwise in all other iterations with k+1 = k the decomposition is unaltered and only a relatively inexpensive back solve is required.

Preconditioning
We experienced that the convergence of Algorithm 1 is negatively affected by the presence of odd mass ratios in the mechanical system. The immediate effect of uneven mass ratio is a broadening of the spectrum of the N matrix. Note that this can be also a consequence of using different measuring units for different constraint multipliers i , for instance even a system with odd mass ratios can have a badly conditioned N matrix if some constraint reactions are assumed in newtons, other in kilo newtons-clearly, as the choice of measuring units is arbitrary, one should find a way to make the solver as much insensitive as possible to the scaling of the constraints. ‡ In many cases the H matrix is hermitian, for example when just block-diagonal mass matrices are used, hence the decomposition of the symmetric indefinite saddle point matrix can be done via a LDL decomposition, which is faster than the LU decomposition. However there are cased where un-symmetric tangent stiffness matrices might be added to H, as happens in some finite element problems, precluding this approach.

Algorithm 1. ADMM for solving CCPs in nonsmooth dynamics
Input: Initial approximations y 0 , z 0 Outputs: Solution ⋆ , v ⋆ , optionally also y ⋆ , z ⋆ 1: Θ 0 = Θ( 0 ) 2: while k ≤ n maxiters do 3: Θ k+1 = Θ( k+1 ) 13: end while A possible remedy for this difficulty is to perform a simplified form of preconditioning by doing a diagonal scaling of N and r before starting the iteration, thus operating on the scaled form of the problem: with̆= S −1 ,N = SNS,ȓ T = Sr T . We set S ii = . Note that other more advanced (and CPU intensive) preconditioners could be used, for example in Reference 39 a modified Ruitz equilibration is proposed. After the solution is computed, one can quickly recover = S̆. Note that in search of high performance we never build explicitlyN andȓ T : what we need is just a scaling of the saddle point matrix involved in (36), that will simply include the scaled jacobiansD = D S, the scaled complianceC = SC S, the scaled termb = Sb , the scaled setΥ = {̆∶ S̆∈ Υ } . The scaled form of the problem would require a projection ΠΥ, however this can be simplified we use the same S ii values for each triplet corresponding to the jth friction cone constraint, for uniform scaling of the jth cone, thus ΠΥ(z) = Π Υ (z). More in general, for a convex set , again assuming a uniform scaling for each  j set, one can use the simplification The preconditioned form of the ADMM algorithm, that we call P-ADMM heretofore, is not reported here for compactness: it is enough to replace the above scaled quantities in Algorithm 1, wherĕ,y,z will replace , y, z, then a final step is added after convergence, to recover the unscaled solution: = S̆. Moreover, if one is interested in the residuals in the not-scaled original form, for instance to verify convergence, these can be computed as As an alternative, one can see that the preconditioned algorithm can be written with the original variables , y, z if one substitutes the scaled variables and performs some algebraic simplifications, obtaining the following steps: The method above corresponds to using the original ADMM algorithm with a nonuniform step Θ k S −2 , with (S −2 ) ii = N ii , hence it avoids the necessity of scaling all variables and matrices. However, when one needs to compute the residuals (for instance for termination criteria or for step scaling policies), in this simplified algorithm (47) those must be computed as:

Warm starting
One of the nice properties of ADMM methods is that they can be easily warm started, something that for example is difficult to do with IPMs. This is a relevant feature in our context, because we must solve the problem (14) at each time step, where one can expect some degree of temporal coherency in the values of between each step. This is especially true for simulations involving stacked objects, like when simulating environments for robots, granular flows or masonry buildings, because the contact forces often show limited changes over time, and the amount of contacts that change state can be limited. If so, one can reuse the last computed values of in the previous time step to warm-start the ADMM solver at the next time step.
This poses two difficulties. First of all, one must develop algorithms and data structures that are able to convey information from the contact at time t A to a time t B = t A + h. In fact, after ADMM computes the contact forces at step t A , those could be saved into the contact data structures; however, those contacts would be recomputed by the collision detection engine at the next time step t B : since the objects move, also the contact manifold changes, thus it is not trivial to compute which contacts are persistent between the two steps. By using tolerances, one can detect if some contacts at time t B are the same contacts used at time t A , if so, their last reaction forces can be copied to the new contacts and used for the warm starting, whereas completely new contacts are initialized with = 0.
Furthermore, one can see that the ADDM Algorithm 1 is a fixed point , thus one needs {y 0 , z 0 } to warm start it, rather than 0 . This leads to two options: store both y and z values in the contact data (wasting some memory in case of large scale simulations, and complicating both data structures and bookkeeping) otherwise assume that the last ADMM run converged exactly up to r primal = 0, so one could just store values in contacts, and then compute This idea, implemented in Algorithm 2, provides a speedup of 2×-10× in our tests.
Note that also the last k value from the previous time step can be used as the initial 0 for the next step, adding a further optimization. Statistically, given some degree of temporal coherency in the mechanical system, such pretuned value would be better than a default value like 0 = 0.1.

ACCELERATED ADMM
Following the idea presented in Reference 35, one can improve the convergence of ADMM by introducing a predictor-corrector step based on the Nesterov acceleration, with minimal computational overhead. Heretofore we will call this variant as FADDM (Fast ADMM). The convergence of the original Nesterov-accelerated ADMM requires the assumption that both f (x) and g(z) in (20) are strongly convex.

Definition 8.
A function f ∶ R n → R, possibly nondifferentiable, is -strongly convex if for all points x, y in its domain it holds: In our case, f (x) is a quadratic function, hence strongly convex, but g(z) is an indicator function, hence convex but not strongly convex. In the case where just one of the two functions is strongly convex, a restart scheme can be applied as in Reference 35, where the method reverts to the original ADMM if the combined residual r k+1 is not monotonically decreasing. The monotone decrease is assessed with a factor r < 1, as a default value we use r = 0.999 to avoid too frequent restarts. Using the same performance optimization strategies that we used to develop Algorithm 1, we finally obtain Algorithm 3.
We note in passing that the Algorithm 3 is based on a z → → y update ordering whereas the Algorithm 1 is based on a → z → y ordering, leading to a slightly modified warm starting method.
Also, a preconditioned version with diagonal scaling can be obtained with the already discussed transformation of (46), leading to a method that we call P-FADMM in the following.

BENCHMARKS
In the literature there are many examples of iterative solvers for nonsmooth dynamics, but many of them are targeting mostly computer-graphics applications, where efficiency and robustness are more important than fidelity and precision. For instance, many benchmarks in computer graphics involve objects with low stiffness (e.g., tissues, ropes, plastics), whereas we need to test our solver with structures of engineering interest, such as steel structures. Also, the proposed method allows an accurate retrieval of the contact stresses, something that is particularly important for mechanics-related applications and that can be neglected in the computer graphics field. On the basis of these requirements we designed the benchmarks that follow. At each time step the ADMM solver requires an updated set of contacts with the corresponding D jacobians: to this end a contact detection algorithm provides an updated set of pairs of potential colliding points between the geometric primitives (edges, faces, curved surfaces) if such points are within a distance threshold. When a deformable structure is present, the contact detection algorithm finds the contact points using the triangles of a tessellated surface that covers the outer skin of the finite elements.
All benchmarks have been computed on a Intel® Core® i7-8750H CPU, clocked at 2.20 GHz, with six physical cores and 16 GB of RAM. The ADMM method has been implemented in C++ within the Chrono open-source library, 43 using the Eigen library v.3.37 for dense and sparse linear algebra, and Intel® MKL Pardiso 44 as a direct solver for sparse linear systems.

High stack with odd mass ratio
This benchmark represents a typical worst case scenario in non-smooth dynamics, where a heavy object sits on a high stack of rigid bodies of much lower mass, as shown in Figure 1. The stable stacking is already a demanding case for contact dynamics, but the presence of extremely odd mass ratios add a further complication.
To this end we designed two subcases: Test 1.a and Test 1.b. Both feature a vertical stack of 20 spheres with mass m = 10 kg each, except the 10th sphere that has a mass of m = 10, 000 kg. A further sphere with mass m = 1000 kg is placed F I G U R E 2 Convergence of alternating direction method of multipliers methods compared to the Jacobi projected fixed point iteration, for Test 1.b The plot shows the violation of unilateral constraints in terms of penetrating residual velocity

F I G U R E 3 Convergence of primal and dual residuals in Test 1a
F I G U R E 4 Convergence of primal and dual residuals in Test 1.b on the floor. All rigid bodies are subject to a vertical gravitational field g = 9.8 m/s 2 . For Test 1.a one has the analytical solution (all contacts are active). The Test 1.b differs from Test 1.a in that the heavy sphere is pulled upward by a force twice its weight, hence causing a separation at the contact below it.
As expected, in both Test 1.a and Test 1.b the projected Jacobi iteration, one of the most common methods in nonsmooth dynamics, has a very bad convergence rate § . As shown in Figure 2, the ADMM method presented in this paper behave much better than the projected Jacobi iteration and, interesting enough, it also converges faster than the Preconditioned Spectral Projected Gradient with Fall-Back (P-SPG-FB), an efficient method for non-smooth problems presented in Reference 17.
For Test 1.a, Figure 3 shows that the ADMM method using the balanced adaptive step size can converge to the exact solution in less than 20 iterations (the primal residual is not visible in the semi-logarithmic plot as it is immediately zero from the first iteration) (Figures 4-6). § It is a known fact that the simulation of high stacks of objects is a major difficulty in non-smooth dynamics, especially in real-time simulators and in video games where the projected fixed-point iterations are truncated to keep the computational time within a limited frame budget. In those cases, the limited precision of the solution will lead to underestimated contact forces, undesired slow interpenetration of the parts and, ultimately, to unnatural collapses of stacks. Note also that we force the update of the step every n s steps because each update of corresponds of a costly factorization of a large matrix. Figure 7 shows that, in this benchmark, the convergence improves for denser updates, meeting the tolerance in the same amount of updates: this would suggest using a small n s anyway. However, for more complex and randomized scenarios, we found that a good tradeoff is n s = 5. This tradeoff value can change depending on the speed of solver used for the factorization: the higher the performance for the back solve respect to the factorization, the higher is the optimal n s .
The FADMM method with the default fixed step = 0.05, instead, converges slowly. Repeating the FADMM test with different step sizes as shown in Figure 6, however, the convergence is greatly improved, suggesting that the FADMM method is sensitive to the initial choice of the time step.
For Test 1.b, convergence plots in Figure 4 show that the FADMM is very efficient even if starting from a fixed not optimal step size, although the performance deteriorates if no preconditioning is used. The ADMM method with balanced step adaptivity n s = 5 converges quite well, especially with the help of a preconditioner. We note that ADMM converges even if primal and dual residuals might exhibit oscillatory behavior (in such cases, the combined residual would be monotonically decreasing anyway 35 ).
Showing that a solver is capable of handling this class of object-stacking problems is relevant, for instance, to the field of civil engineering, because it shares the same difficulties of simulating tall masonry structures when assessing stability and seismic response via discrete elements. 45

Wrecking ball
We modeled four parallel walls made with 10 rows of 15 bricks each. Bricks length, height, width are, respectively, 3.96, 2, 4 m, and their density is 100 kg/m 3 . A wrecking ball with diameter 8 m and density 8000 kg/m 3 impacts horizontally with the 600 bricks. The density of the bricks respect to the ball is deliberately high in order to generate a badly conditioned problem with odd mass ratios. The friction coefficient is = 0.4 and the time step is h = 0.02 s. Figure 7 shows a snapshot from the simulation. The convergence of the ADMM method, as shown in Figure 8 for a random time step, is noticeably superior to the convergence of fixed point iterations such as projected Jacobi. The method exhibit even better convergence when using the diagonal preconditioning (46), as shown Figure 9, because of the odd scaling of the masses ( Figure 10).

Robotic manipulator
This benchmark features a 6-DOF industrial robot whose end effector interacts with 15 boxes with height 0.3 m and width 0.4 m, moving them in the working area via contacts and collisions, as shown in Figure 11. The friction coefficient between the boxes is = 0.4, the time step is h = 0.01 s. The difficulty highlighted by this benchmark is the simultaneous presence of nonsmooth contacts and bilateral constraints in a critical articulated mechanism. In fact the simulation of a robot arm actuated via motorized joints is a trivial problem, but here we do the opposite: we drive the end effector via a rheonomic constraint defining an imposed trajectory respect to the base, hence the rest of the arm will naturally move to the configuration prescribed by joint constraints; however, since the mass of the end effector is much lighter than the bicept, since some joints are redundant ¶ and since the robot passes close to a singularity, the simulation of the robot alone is a challenging task. That is, even without contacts, most iterative solvers like SOR or Jacobi would converge too slowly, and even Krylov linear solvers like GMRES, CG, MINRES would require many iterations unless equipped by some preconditioner. In our test, convergence to high precision is achieved in 15 iterations at most, where just one or two of them require a refactorization of the saddle-point matrix. Moreover, when there are no contacts, the solution is provided in exactly one step, as the iteration boils down to a single factorization followed by a forward/back solve just like when simulating the robot via a direct solver in the conventional context of smooth dynamics.

Deformable bars
This benchmark aims at estimating the efficiency and robustness of the method in problems involving also finite elements. Few methods already exist in the literature that can simulate deformable structures in the context of nonsmooth dynamics, but most of them target the field of interactive computer graphics, where the stiffness of structures is often very low. In our case, being interested in engineering applications, we require that the method would be able to handle both very deformable structures (ex. rubber-like materials) as well as very stiff structures. In this benchmark, we simulate the fall of 15 deformable bars, each modeled with a mesh of tetrahedral finite elements, for a total of 9360 degrees of freedom. At initial state bars are not touching and are dropped in five groups of three, each group with random yaw rotation and with increasing height in the 0.1 … 0.5 m range. The outer skin of the bars can collide with each other and with rigid bodies-in this case there is just a single rigid body, namely the flat ground. At the contact points there is a friction coefficient = 0.3 and a null restitution coefficient, the time step is h = 0.05 s. In all simulations the material of the bars has a density of 1 × 10 3 kg/m 3 , a 0.01 Rayleigh damping coefficient and a Poisson ratio = 0.3, but in the case of Test 4.a (see Figure 12) the Young modulus is E = 1 × 10 6 Pa, whereas in Test 4.b (see Figure 13) the Young modulus is many orders of magnitude higher: E = 2 × 10 11 Pa. Tests has shown that our ADMM method converges equally well even in the case of extreme stiffness, converging to the required tolerances with comparable number of iterations. At each time step, both Test 4.a and Test 4.b required an average number of three factorizations for the adjustment of the i step and a variable number of forward/back solves in the range of 4 ÷ 50. When bars come into contact, an average ¶ The robot has a gravity compensator modeled via two revolute joints and one cylindrical joint, leading to an overconstrained loop.

RESULTS AND CONCLUSION
We performed benchmarks involving multibody systems with contacts between multiple parts, showing that the performance of the ADMM method is capable of handling problems that would converge too slowly using conventional projected fixed point methods or first-order spectral methods. Our ADMM method requires few computational primitives: basically a projection of dual variables on conic sets, a forward/backward solve of a linear system, and its factorization. The latter is a computational bottleneck, but it can be performed only once per run, as the matrix does not change often during the iterations. We noted that ADMM can be successfully applied to problems that exhibit temporal coherence because, unlike IPMs, it supports warm-starting. A good estimation of the ADMM step size proved to be fundamental in achieving good convergence: using adaptive step scaling we obtained an efficient auto-tuning algorithm. Another optimization that allowed superior performance is the adoption of a diagonal preconditioning, with block scaling for the triplets of lagrangian multipliers relative to the conic constraints.
Further research on this topic might address the acceleration of the ADMM method by means of Anderson acceleration and the solution of nonconvex problems with complex and non-associated frictional models that go beyond the standard Coulomb model.