Second order cone programming for frictional contact mechanics using interior point algorithm

,


Introduction
Contact problems with dry Coulomb friction are present in many design and validation processes in mechanical engineering.As soon as several objects are involved, the question of computing contact forces arises.Examples include multibody systems and mechanisms (Brogliato, 1999;Pfeiffer and Glocker, 1996), robotic systems and grasping problems (Bicchi and Kumar, 2000;Buss et al., 1997Buss et al., , 1996;;Han et al., 2000), deformable solid mechanics (Johnson and Johnson, 1987;Popov et al., 2010;Wriggers, 2006), and granular materials (Cambou et al., 2013;Radjai and Dubois, 2011).A common model of one-sided contact is to consider a set of inequality constraints on the configuration parameters of the mechanical problem (positions, rotations) associated with forces which must themselves be positive.Coulomb's friction, which governs how objects slide relative to each other, imposes a constraint on the contact forces that must remain within the Lorentz cone.These laws lead us to write second order cone constraints on the contact forces.By introducing a non-linear change of variables of the contact relative velocities which defines the so-called modified velocities, the contact laws can be written as a complementarity condition on second order cones between the modified velocities and the contact forces. 1 In solid and structural dynamics, the discrete mechanical model is usually supplemented by equations of motion that relate the velocities of the system to the applied and contact forces.After discretization in space and time, the resulting system can be expressed as nonlinear Second Order Cone Complementarity Problem SOCCP (Acary et al., 2011;Acary and Brogliato, 2008;Acary and Cadoux, 2013).Section 2 details the problem formulation when the discrete equation of motions are assumed to be linear.
To solve these SOCCP problems numerically, a large number of methods have been used in the computational contact mechanics community.Many of these methods are in fact adaptations of well-known mathematical programming methods for solving variational inequalities and complementarity problems.Some of these methods have been jointly developed by optimization specialists.There are two inherent difficulties in frictional contact mechanics problems.The first difficulty is related to the non-monotone nature of the complementarity problem.Identifying an optimization problem for which the complementarity problem would be the optimality conditions is difficult, and leads to non-convex optimization problems.The second difficulty arises from the fact that the constraints are not of full rank and, in the case of rigid multi-body systems, may be of very low rank with respect to the number of constraints.In (Acary et al., 2018), a review of the main literature is made and methods are compared with performance profiles.When the second-order constraints are of full rank, second-order methods such as semi-smooth Newton methods, e.g., Alart-Curnier's method (Alart and Curnier, 1991), are generally robust and accurate.If not, as it is usual in multi-body systems made of rigid parts (robots, granular material), then the first-order methods such as projected successive over-relaxation gradient methods are robust, but slow and with a limited accuracy in practice.Second order techniques cannot solve the problem due to robustness issues.As far as we know, there is no high accuracy second order method capable of solving the friction contact problems with redundant constraints.
One of the objectives of this work is to propose a second-order method, based on an interior point method that is accurate, robust and efficient for problems where the constraints are rank deficient, but, in the first step, on a convex optimization problem.Some applications of interior point methods for contact problems have already been attempted in the literature (Kleinert et al., 2014;Kučera et al., 2013;Mangoni et al., 2018;Melanz et al., 2017).In the precursor work of (Kučera et al., 2013), the contact problem with Tresca friction (purely quadratic cylindrical constraints) is solved with Mehrotra's algorithm.The work in (Kleinert et al., 2014) is the most advanced for the case of conic constraints.The problem considered is the relaxed convex problem as in this work, which is further regularized by adding a constant diagonal term to the Jacobian matrices of the interior point algorithm.To solve large problems, linear systems are solved by iterative methods.The degeneracy of the conditioning as the iterations progress leads to the use of preconditioners and makes it costly to obtain high accuracy solutions.In (Melanz et al., 2017), a comparison is made between an interior point method and first order methods.The interior point method on second-order cone shows a better convergence rate but without reaching higher accuracies than first order methods.
In this article, we consider a relaxed problem where the actual local velocity is complementary to the reaction force at each contact point.This yields a reformulation of the original problem as a convex second order cone optimization (SOCO) problem.A dedicated interior point method based on the primal-dual algorithm of Mehrotra with the Nesterov-Todd scaling strategy is tailored to solve this problem.In particular, we show on a large bench set of examples that a tight accuracy can be achieved at optimum without regularizing the Jacobian matrices, provided that their conditioning is controlled.
Since our main objective is to have a robust and efficient computer code capable of solving our problem with a high degree of accuracy, we have first performed experiments with existing solvers.At first, we have tried to solve the convex relaxation problem by using a differentiable nonlinear optimization as proposed by Benson and Vanderbei (Benson and Vanderbei, 2003).They suggest four ways of reformulating the second-order cone constraint to get round the non-differentiability of the norm at zero.Our experiments were performed by means of the modeling langage AMPL (Gay, 2015) and the optimization solver KNITRO (Byrd et al., 2006).Our conclusion is that, except in a few rare cases, none of these reformulations is able to solve efficiently and accurately the convex model of our problem.The second experiments were with the interior point solver SDPT3 (Tütüncü et al., 2003).This is an implementation of the Mehrotra predictor-corrector algorithm, which we hoped would be effective and robust in solving our problem.From our point of view, the great advantage of SDPT3 is that it is open source and we can therefore hope to modify it to solve the initial problem.Unfortunately, the results of the experiments with SDPT3 were not conclusive.In particular, the accuracy we wanted to achieve was not sufficient and too many failures were observed.But the experiments were interesting for us, because we identified that the main difficulty in implementing this algorithm is in solving the linear system at each iteration.This led us to carry out our own implementation of this algorithm, devoting our efforts in solving the linear system and in the scaling method.We also made experiments with the interior-point solver MOSEK (Andersen and Andersen, 2000) whose algorithm is based on a homogeneous and self-dual model (Andersen et al., 2003).The numerical results of MOSEK are better than those of SDPT3.However, in order to obtain accurate solutions, the number of iterations can become large or failures occur as soon as the stopping tolerance is reduced (below 10 −7 ).Unfortunately, MOSEK is not open source, so it has not been possible to understand the origin of these failures, nor will it be possible to modify this code to solve our initial problem.These first numerical experiments motivated us to make our own implementation of the interior point algorithm.Especially since in the conclusion of (Tütüncü et al., 2003) the authors state that "For problems with second-order (quadratic) cone constraints, experiments indicate that there is room for improvement in SDPT3".
The outline and the contribution of the article are as follows.In Section 2, we recall the basics of mechanical models and how the convex SOCO problem is obtained.The corresponding primal-dual pair is formulated in Section 3. General results are given on the non-emptiness and compactness of the solution set under the Slater hypothesis.In Section 4, the properties of the central path are detailed.In particular, under the assumption of strict complementarity, we establish a characterization of the limit point of the central path, the analytic center of the optimal set, which to our knowledge is new in the SOCO context.This property of the algorithm is important from the mechanical point of view, since the selection of dual variables, that corresponds to the reaction forces, is completely controlled by the interior point method.The details of the numerical implementation of the algorithm are given in Section 5. Several alternative equivalent formulations of the linear system to be solved at each iteration are detailed and the comparison of their conditioning over the iterations is illustrated.Our experiments show that the choice of formulation is fundamental for the robustness of the algorithm.In Section 6, the interior point method is extended to the case of rolling friction, where the cone of constraints is no longer a Lorentz cone and is not self-dual.

Notation
Let x and y be two vectors of R n .The Euclidean scalar product of x and y is denoted x ⊤ y.The associated ℓ 2 norm is ∥x∥ = √ x ⊤ x.The perp notation x ⊥ y means that x ⊤ y = 0. Let E be a subset of R n .The interior, the relative interior and the boundary of E are respectively denoted by int(E), ri(E) and bd(E).The dual of E is denoted and defined by The definitions and notations related to the Jordan algebra are described in Appendix A.

Mechanical models
Let us first introduce the original mechanical models that we are interested in.Let d be the dimension of the space of the mechanical system.Two classes of problems are considered.The frictional contact (FC) problems (Acary et al., 2018) for which d = 3 and the problems with rolling friction (RF) at contact (Acary and Bourrier, 2021) for which d = 5.Let n ∈ N be the number of contact points and let m ∈ N be the number of degrees of freedom.The mechanical system is described by means of three vectors: the global velocity v ∈ R m , the velocity u ∈ R nd and the reaction r ∈ R nd .After discretizing the dynamical system in time and space, the general model to be considered is a conic complementarity problem of the form where the data of the problem are M ∈ R m×m a symmetric and positive-definite matrix, f ∈ R m , H ∈ R nd×m , w ∈ R nd , K is a cone and K * := {u ∈ R nd : u ⊤ r ≥ 0, for all r ∈ K} is its dual cone.The cone is of the form K = n i=1 K i , each K i being a cone whose definition depends on the model.For a FC problem, we have This is the second order cone of Coulomb friction at the ith contact point.The scalar r i,N (resp.u i,N ) and the vector r i,T (resp.u i,T ) are the normal and tangential components of the reaction (resp.velocity) vector of the ith contact point and µ ∈ R n is the vector of friction coefficients.The dual cone of K, is given by For a RF problem, the cones are of the form The vector r i,R is the rolling friction reaction moment at contact and µ R ∈ R n is the vector of rolling friction coefficients.The dual cone of K i is The vector u i,R is the relative angular velocity at contact.For this model, the components of the function Φ are defined by Φ By observing that Φ(u) = Φ(u + Φ(u)) and by making the change of variable u ← u + Φ(u), the problem (1) is reformulated as (2) In (Acary et al., 2011), for the case of a FC problem, an iterative solution of ( 2) is proposed as follows.The second equation is rewritten as where ϕ : R n → R nd is defined by n components of the form (s i , 0 ⊤ ) ⊤ ∈ R d .The idea is to solve the nonlinear system as a parametric system of the form where s is periodically updated according to There is no guarantee of global convergence of this procedure, but the advantage of this formulation is that, for a fixed s, the parametric system corresponds to the first order optimality conditions of the following convex optimization problem: min This problem belongs to the class of second order cone optimization (SOCO) problems.Our study focuses on the numerical solution of this problem.Obviously, the same parametric strategy can be applied to the solution of an RF problem.

Basic properties
To simplify the models, let us remove the friction coefficients by means of a simple change of variables.Let us define the diagonal matrix P µ whose ith diagonal block is diag(1, µ i , µ i ) for the FC model and diag(1, µ i , µ i , µ R,i , µ R,i ) for the RF model.The parameter s is assumed to be fixed.Let us make the change of variables: u ← P µ u, H ← P µ H, w ← P µ (w + ϕ(s)) and r ← P −1 µ r.
We keep the same names for the variables, but change the notation of the components.For the RF model we define The same notation is used for the FC model, except that the components ũi and ri do not appear.The system (3) becomes where F is the product of cones.The definition of F depends on the model.For the FC model (d = 3) we have This cone is self-dual (i.e., L * i = L i ), therefore For the RF model (d = 5), the cone of constraints is This cone is not self-dual.The dual cone is then The system (4) can be viewed as the first order optimality conditions of the following second order cone optimization problem: min Let r ∈ R nd be the dual variable associated to the linear constraint of (9).The Lagrangian function associated to problem (9) is defined by Since L is separable in v and u, convex in v and linear in u, and since F * * = F, the dual function is readily obtained and defined by The dual problem is then min Let us denote by R the optimal set of this problem.It is a convex set, characterized by the first order optimality conditions: Note that these conditions are equivalent to (4), thanks to the definitions (10) and the equalities By weak duality, the duality gap (i.e., the difference between the value of the primal and the dual function) is nonnegative.Indeed, let (v, u, r) be a primal-dual feasible solution, we have ≥ 0, the inequality arising from the fact that r ∈ F and u ∈ F * .The strong duality and the existence of an optimal solution for (11), we need a constraints qualification assumption, known as the Slater hypothesis.
With respect to the reduced problem ( 11), the Slater hypothesis can be equivalently formulated as follows: The following result summarizes the connection between the solutions of ( 4) and ( 11).They are well known and are derived from basic properties of convex optimization.We provide a proof in Appendix B to be complete.Proposition 3.2.Consider the primal-dual pair of SOCO problems (9)-( 11), where M ∈ R m×m is symmetric and positive-definite, f ∈ R m , H ∈ R nd×m , w ∈ R nd and F is the product of second order cones of the form (5) or (7).
(ii) If the problem (9) is feasible, then it has a unique optimal solution (v, û).(iii) Assumption 3.1 is satisfied if and only if the optimal set of (11) is non-empty and compact.

Central path related to a frictional contact problem
In this section, we are interested in FC problems.We recall some known results about the convergence of the central path in SOCO and put them in our context.For the sake of completeness we also provide the proofs.
For FC problems, the friction contact cone is equal to the product of n Lorentz cones, i.e., F = L n , where L is the Lorentz cone of dimension d.We recall that this cone is self-dual, i.e., F * = F.There is an algebra, the so-called Euclidean Jordan algebra, associated to symmetric cones, which allows an almost direct extension of the interior-point algorithms for linear optimization to the case of SOCO (Alizadeh and Goldfarb, 2003).A summary of Euclidean Jordan algebra is given in Appendix A. Thus, the orthogonality condition in (12) becomes u • r = 0, see (Alizadeh and Goldfarb, 2003, Lemma 15).This leads to a square system of equations with the additional conical constraints.Under Assumption 3.1, (r, u) is a primal-dual optimal solution of the problem (11) if and only if It is important to keep in mind that the matrix W is not an arbitrary positive semidefinite matrix, but has a structure given by ( 10).This structure will be useful for solving the linear system done at each iteration.In interior-point methods, a perturbation is introduced in the complementarity equation, so that the system to be solved becomes where e is the unit vector associated with the Jordan product and µ > 0 is a parameter that is driven progressively to zero along the iterations, so that in the end a solution of ( 14) is found.The parameter µ is called barrier parameter, because (15) can be interpreted as the optimality condition of the barrier problem min The first order optimality condition of ( 16) is W r + q − 2µr −1 = 0.By introducing the variable u = 2µr −1 , we get (15).This system (15) defines a curve, called the central path.The following result states that under Slater's hypothesis, the central path is well defined and remains bounded for bounded values of µ.The proof is given in Appendix C Proposition 4.1.Under Assumption 3.1, for all µ > 0, the perturbed KKT system (15) has a unique solution (r(µ), u(µ)) ∈ int(L 2n ).For all μ > 0, the set {(r(µ), u(µ)) : 0 ≤ µ ≤ μ} is bounded.
From Proposition 3.2-(ii), the optimal solution of ( 9) is unique, therefore lim µ→0 u(µ) = û.For the curve r(•), it can be shown that lim µ→0 r(µ) = r, where r ∈ ri( R).Such a solution is called a maximally complementary optimal solution of (11).It can be characterized as follows.Regarding to the problem (11), the index set {1, . . ., n} of the Lorentz cones is partitioned into six sets (Mohammad-Nezhad and Terlaky, 2021): An optimal solution r ∈ R is maximally complementary if and only if for all i ∈ B, r i ∈ int(L) and for all The next result shows the convergence of the central path to a maximally complementary solution.A proof is given in Appendix C for self-completeness.
In linear optimization, it is known that the central path converges to the analytic center of the optimal set (Monteiro and Zou, 1998).In semidefinite optimization, it is also known that this result holds under strict complementarity (De Klerk et al., 1998;Halická et al., 2005).But in SOCO, as far as we know, we have not found in the literature any characterization of the analytic center, even under the strict complementarity assumption.In order to complete the theory, we provide such a characterization.
This assumption is equivalent to T i = ∅, for i = 1, 2, 3.In that case, (B, N, R) is a partition of {1, . . ., n}.Under this assumption, we can define the analytic center of the optimal set R. If R is a singleton, then the analytic center of R is this point.Otherwise, it is the unique optimal solution of the problem The next proposition shows that the analytic center is well defined.
Proposition 4.4.Under Assumptions 3.1 and 4.3, if R is not reduced to a singleton, then the problem (17) has a unique optimal solution r ∈ ri( R), characterized by the following property: for all r ∈ ri( R), i∈B Proof.Let r ∈ ri( R).For all i ∈ B, det r i > 0, and for all i ∈ R, To show that problem (17) has at least one solution, it suffices to show that the function ψ + δ ri( R) is coercive.This last property is a direct consequence of the compactness of the set R.
The uniqueness of the minimum comes from the strict convexity of ψ on ri( R).Indeed, for a conic component i ∈ B, we have ∇ 2 r (− log det r) r=ri = 2Q r −1 i , which is positive definite for all r i ∈ int(L), see Appendix A. For all conic components i ∈ R, there exists h i ∈ R d−1 , with ∥h i ∥ = 1, such that for all r ∈ ri( R), we have r i = r i,0 (1, h ⊤ i ) ⊤ .Indeed, suppose that for i ∈ R, there exist r and r ′ in ri( R), such that r i and r ′ i are not collinear.By the triangle inequality, it is easy to see that 1 2 (r i + r ′ i ) ∈ int(L), and thus i ∈ B, which would contradict i ∈ R. Finally, for all r, r ′ ∈ ri( R), if r ̸ = r ′ then there exists i ∈ B such that r i ̸ = r ′ i or there exists i ∈ R such that r i,0 ̸ = r ′ i,0 .In both cases, we have ψ(r) ̸ = ψ(r ′ ), which implies that ψ is strictly convex on ri( R).Since ψ is convex on ri( R), r is optimal if and only if −∇ψ(r) is in the normal cone to R at r, i.e., for all r ∈ ri( R), ∇ψ(r) ⊤ (r − r) ≥ 0. For r ∈ ri( R), we have from which we deduce that ( 18) is satisfied.
Let us state and prove the main result of this section.
Theorem 4.5.Under Assumptions 3.1 and 4.3, the central path r(•) converges to the analytic center of R.
Proof.Suppose that R is not reduced to a singleton, otherwise the result is a direct consequence of Proposition 4.2.Let r ∈ ri( R).As in the proof of Proposition 4.2, for all µ > 0, (44) is satisfied.By using r(µ) • u(µ) = 2µe and the definition of the partition (B, N, R), we have i∈B For any index i, by using the Cauchy-Schwarz inequality and the fact that r i,0 ≥ ∥r i ∥, we have In the same manner, for all i ∈ {1, . . ., n} we have From ( 19), ( 20) and ( 21), we deduce that i∈B When µ tends to zero, the first sum tends to i∈B r ⊤ i r−1 i and the second one to |N|.By using the fact that ri,0 = ∥ ri ∥ and ûi,0 = ∥ ūi ∥ for i ∈ R, each term of the third sum tends to ri,0 2ri,0 + 1 2 .By taking the limit as µ → 0 and using n which shows by Proposition 4.4 that r is the optimal solution of ( 17), the analytic center of R.
We would also like to emphasise an important point concerning the characterization of the analytic center from the point of view of computational mechanics.If the matrix H is not of full row rank, the reaction forces r are not unique.The interior point algorithm allows the optimal reaction forces to be uniquely defined.In addition, the analytic center is of mechanical interest and interpretation.By choosing a solution that is furthest from the boundaries of the solution set and therefore from the boundaries of the cones, the chosen solution maximizes the sum of the normal stresses r 0 and the distance to the edges of the cones.For simple cases, this provides a better distribution of reactions for active constraints rather than particular solutions where reaction forces are concentrated on a few constraints.

Numerical solution of the friction contact problems
We have implemented the primal-dual interior point algorithm developed by Tütüncü, Toh and Todd (Tütüncü et al., 2003) and adapted it to our context.The algorithm is an extension of the predictor-corrector algorithm of Mehrotra (Mehrotra, 1992) to the solution of a SOCO problem.The differences with the algorithm implemented in SDPT3 are as follows: • The convex quadratic objective function is taken into account directly, without reformulating the problem with a linear objective function and a new conic constraint.One consequence is that the primal and dual steps must be equal, which is not necessarily the case with a linear objective function.This is explained below.
• The linear system to solve at each iteration, resulting from the linearization of the perturbed KKT system (15), is transformed via the Nesterov-Todd scaling strategy as in SDPT3, but in a different manner.
The scaling matrix (i.e., the quadratic representation of the vector p defined by formula (25) below), which introduces bad conditioning when the iterations are close to an optimal solution, is never explicitly computed and stored in memory.Only matrix-vector products are performed.This method allows us to achieve a high degree of accuracy in the computation of the directions and thus in the solution of the quadratic convex problem.We examine in detail equivalent formulations of the linear system to justify our choice and detail all the computations performed.
• We also exploit the data structure of the problems, by using the matrices H and M instead of the reduced matrix W (see formula (10)), in order to preserve the sparsity of the linear system.This transforms the initial 2 × 2 block linear system into a 3 × 3 block system.We also show, by numerical experiments, that it is better to reduce again this system back to a 2 × 2 block system, while preserving the sparse structure, but avoiding some of the calculations with the scaling matrix.
• Finally, the experiments show that using the C data type long double for all the calculations related to the Nesterov-Todd scaling scheme, improves the robustness of the implementation, while maintaining reasonable computation times.
The main part of the algorithm is the solution of two linear systems that result from the linearization of the equation ( 15) at the current iterate (u, r) ∈ int(L 2n ).They only differ on the right-hand side and are of the form: where U = Arw(u) and R = Arw(r) are the arrow-shaped matrix associated to u and r.See Appendix A for all definitions related to the Euclidean Jordan algebra that are used in this paper.The first direction, denoted (d u a , d r a ) and called the affine scaling direction, is the solution of ( 22) without the square bracketed term in the right-hand side.It then satisfies the linear equation The affine scaling direction is a Newton step on the original optimality system (14).The barrier parameter is set to µ = u ⊤ r n .The centralization parameter σ ∈ (0, 1] is fixed by comparing the current value of µ with its expected reduction obtained along the affine step.The second direction is a linear combination of the affine scaling direction and of a corrector step, to keep the iterates near the central path.It then satisfies the following linear equation u , where α is the largest value in (0, 1] such that for some value τ ∈ (0, 1] (typically τ = 0.995).Contrary to common practice, different primal and dual steplengths are not taken, because of the first equation in ( 22), which is linear and includes both primal and dual variables.Indeed, suppose that the current iterate is such that W r+q−u = 0 and that different steplengths α ̸ = α ′ are taken with d u ̸ = 0. We then have If d u ̸ = 0, then at the next iteration the residual of the linear equation is no longer zero.
The algorithm will be well defined if the matrix in ( 22) is nonsingular at each iteration.Its determinant is equal to det(W + R −1 U ) det R.Although the vectors r and u are kept inside the interior of the second order cones, and thus R and U are positive definite, the matrix W +R −1 U can be singular.This is because the matrix R −1 U is not necessarily symmetric, since in general r and u do not commute.The following example is given by (Peng et al., 2002, p.143 In addition to this singularity problem, there is the symmetry problem.In the context of interior-point algorithms for linear or nonlinear optimization, where the cone of constraints is the nonnegative orthant, the matrices U and R are diagonal and so the matrix of ( 22) can be symmetrized, for example by left-multiplying the second row by −U −1 .See, e.g., (Ghannad et al., 2022) for several symmetrization techniques in interior point methods.The main advantages of a symmetric system are lower factorization costs and an effective control of the inertia of the factorized matrix.In addition, very efficient codes such as MA57 (Duff, 2004) or MUMPS (Amestoy et al., 2001) can be used for this task.
To overcome these problems of singularity and symmetry, a change of variables, called a scaling scheme, is used to obtain a symmetric nonsingular system.The idea is to make a change of variables that leaves the Lorentz cone invariant and such that in the new space the vectors u and r commute.However, this change of variable depends on the current iterate and must be done at each iteration.Let p ∈ int K and let Q p be the associated quadratic representation (see Appendix A).Since p is in the interior of the cone K, the matrix Q p is positive definite.From (Alizadeh and Goldfarb, 2003, Theorem 9), we have The problem (11) becomes min The corresponding perturbed KKT system is W r + q = u and û • ř = 2µe, with (û, ř) ∈ int L 2n .The linearization of these equations leads the following linear system: where Û = Arw(û) and Ř = Arw(ř).The choice of the vector p ∈ int K is made so that û and ř commute, which implies that the matrix of the linear system (24) is nonsingular.Indeed, this matrix is nonsingular if and only if det(W + ( ŘQ p ) −1 Û Q p −1 ) ̸ = 0. Since Û and Ř are positive definite, û and ř commute, and which is symmetric and positive definite.
Several choices for the vector p are possible, see (Alizadeh and Goldfarb, 2003).As mentioned in (Tütüncü et al., 2003), the most efficient scaling technique for the solution of a SOCO problem, is the one using the Nesterov-Todd (NT) direction (Nesterov and Todd, 1997): The main property of the NT direction is that û = ř, which implies that The symmetrization of the system ( 24) is done by left-multiplying the last row by −Q p Û −1 , leading to a symmetric matrix of the form To take advantage of the sparsity of the matrices M and H, the system ( 22) is considered in the following equivalent augmented form: The system (26) can be interpreted as the linearization of the perturbed KKT conditions of the problem (9).Applying the scaling scheme, the system becomes A reduction of ( 27) can be done by eliminating the variable d u , while keeping the sparse structure.This leads to the reduced symmetric system The big flaw of the scaling strategy is the ill-conditioning of the matrix Q p 2 when the solution pair (u, r) approaches an optimal solution.Indeed, suppose that (u * , r * ) is a primal-dual optimal solution of (11), which satisfies the strict complementarity condition.Let (u, r) be an interior point iterate near the optimal solution and let p be defined by (25).For i ∈ {1, . . ., n}, three situations can occur (Cai and Toh, 2006): • u * i ∈ int(L) and r * i = 0, then all the eigenvalues of Q p 2 i are of order µ := u ⊤ r; • u * i = 0 and r * i ∈ int(L), then all the eigenvalues of Q p 2 i are of order 1/µ; • u * i ∈ bd(L), r * i ∈ bd(L) and (u * i , r * i ) ̸ = (0, 0), then the largest eigenvalue of Q p 2 i is of order 1/µ and the smallest is of order µ.
To overcome the difficulties due to ill-conditioning, we propose to solve the linear systems ( 27) and ( 28) under the following equivalent form: where In our numerical experiments, we also consider the reduced system for which the matrix is positive definite.
Figure 1 shows the behavior of the condition number of the matrices of the six linear systems ( 26)-( 31) along the iterations of the interior-point algorithm for two examples.The first example (left figure) has a single contact point: . Since H is of full rank, the primal-dual solution is unique, the optimal reaction and relative velocity vectors are non-zero and on the boundary of the Lorentz cone.
The second example is a "Box Stacks" problem from FCLIB, a collection of discrete three-dimensional frictional contact problems (Acary et al., 2014) 2 .There are n = 69 contact points and m = 450 degrees of freedom, H ∈ R 207×450 and rank(H) = 157.The optimal solution satisfies the strict complementarity condition and (|B|, |N|, |R|) = (18,5,46).The matrix of (26) at the endpoint of the minimization procedure, is nearly rank deficient, there are 15 singular values less than √ ϵ, where ϵ is the machine epsilon.These two examples show that the matrices in ( 29) and ( 30) remain the least ill-conditioned.
An advantage of the systems ( 29) and ( 30) is that the quadratic representation matrices are never explicitly built in memory for the entire computation.Only the product of these matrices times a vector needs to be performed.Indeed, for a pair of vectors (x, y) ∈ L 2 , the product of a vector y by the quadratic representation of x can be done via the formula Q x y = 2(x ⊤ y)x − (det x)R d y.Therefore, the product of a vector by a matrix Q p , where p is the , can be done by performing only three products of a quadratic representation matrix by a vector.Similarly, the computation of the inverse or the square root of a vector in the Jordan algebra, is done by using the spectral decomposition of that vector.Moreover, even if the number of cones can be large, the computational cost of a spectral decomposition per cone is very small, since the dimension of a Lorentz cone is only three.32) is satisfied then return (v, u, r) as primal-dual solution of ( 9);
The linear system ( 26) is solved by means of a LU factorization, the symmetric ones with a LDL ⊤ factorization with MA57 (Duff, 2004).Even for the solution of the positive definite system (31) MA57 is used.Two types of failures are returned during a run: • Failure 1: the stopping criterion (32) is not satisfied after a maximum of 100 iterations.
• Failure 2: A NaN (Not a Number) is detected during the computation of the new iterate.
Table 2 indicates the number of successes and failures when solving the 1091 problems of the FCLIB collection, with a tolerance fixed to tol = 10 −10 .Each row corresponds to a run of Algorithm 1 with the numerical solution of the indicated linear system.Figure 2 shows the corresponding performance profiles (Dolan and Moré, 2002).With the system (26) the failures are due to a nearly singular system.In these cases, either the algorithm stalls to a spurious solution (13 out of 22 cases) or the convergence becomes very slow (9 out of 22 cases).It should be noted, however, that the "no-scaling" strategy yields an optimal solution for almost 98% of the problems.The systems ( 27) and ( 28) return a large number of failures of type 2. This is mainly due to the ill-conditioning of the matrix Q p 2 when approaching an optimal solution.The reduction of the system worsens the results.Surprisingly, the worst results are obtained with the system (29).A deterioration of the residual of the second linear equation in (29) over the iterations is observed, when the matrix Q p becomes increasingly ill-conditioned.This leads to a loss of the primal feasibility of the iterates.This is mainly due to the scaling of the linear equation −Hd v + d u = u − Hv − w.To overcome this problem, the refinement procedure described in the MA57 documentation can be used to solve (29).We performed a run with a refinement tolerance fixed to the tolerance tol and a maximum of 10 refinement iterations.This results in only six type 2 failures and no more failure of type 1, but it takes longer to run than with (30) as shown in Figure 2. The best performance in terms of robustness is obtained with the system (30).No failures were detected.The average number of iterations to solve one of the 1091 problems is equal to 18, while the minimum and maximum numbers are equal to 8 and 34.The positive definite system (31) gives a good performance in terms of efficiency, but its robustness is not sufficient, even when refinement is applied.2: Performance profiles for seven different linear system choices.For τ ≥ 0, ρ s (τ ) is the fraction of problems for which the performance of a given version of the algorithm is within a factor 2 τ of the best one.
At last, it should be mentioned that all the computations related to scaling are done in C using the floating point datatype long double.This data type is also used in the computation of the steplengths.This always results in a better accuracy, although the computing time is slightly longer.Table 3 shows the results of solving the FCLIB problems with Algorithm 1, the linear system under the form (30), with a stopping tolerance tol = 10 −11 .The comparison is between the use of the datatype long double versus double.We can see that the number of failures of type 2 is more than doubled with the type double.For these runs, even in case of a failure, the residual (left term in (32)) is small, which means that an optimal solution has been found.The column max res indicates the maximum value of the residual norm (32) of all the 1091 residuals and shows that the type long double allows to obtain a better precision.Even with the data type double, all problems are successfully solved for tol = 10 −10 .The last column of this table shows the total computational time needed to solve all the FCLIB problems with this tolerance.It can be seen that the increase in computational time is less than 10% with long double.Table 3: Performance comparison double versus long double when solving the FCLIB problems with system (30) and tol = 10 −11 .

Rolling friction contact problem
We now consider the solution of (4) in the framework of the RF model (d = 5) defined by the rolling friction cones (7) and ( 8).The main difficulty is that an elementary cone R i is no more self-dual and therefore not symmetric.There is no Jordan product such that R i is a cone of squares with respect to this product.A potential approach is to transform the cone of constraints related to problem (9) into the product of Lorentz cones by means of the following usual trick of introducing artificial variables.For real numbers, a ≥ b + c if and only if there exist two real numbers t ≥ b and t ′ ≥ c such that a = t + t ′ .By setting u i,0 = ti + ti for all i ∈ {1, . . ., n}, the primal-dual pair of problems ( 9)-( 11) can be rewritten as and min where d+1) is a block diagonal matrix with n blocks of the form d+1) .
For all i ∈ {1, . . ., n}, we have ȷz i = ( ti + ti , ū⊤ i , ũ⊤ i ) ⊤ and ȷ ⊤ r i = (r i0 , r⊤ i , r i0 , r⊤ i ) ⊤ .The perturbed KKT system associated with the problem (34) is then This optimality system is associated with the barrier problem min r φ µ (r) := 1 2 r ⊤ W r + q ⊤ r − µ n i=1 log(det(r i,0 , ri ) det(r i,0 , ri )). ( 36) As for Proposition 4.1 and by coercivity of the barrier function, under Assumption 3.1, for all µ > 0 the system (35) has a unique solution such that (z(µ), J ⊤ r(µ)) ∈ int(L 4n ).Although the optimal solution of ( 33) is not unique, because J is non-injective, it can be shown, like for Proposition 4.2, that the central path converges to a relative interior point of the primal-dual optimal set.Under the hypothesis of strict complementarity, it can also be shown that the central path r(•) converges to the analytic center of the dual optimal set.For the sake of completeness, we state the result, but without proof in order to lighten the paper.
As in Section 4, under the strict complementarity assumption, the index set {1, . . ., n} can be partitioned into two partitions ( B, N, R) and ( B, Ñ, R), whose definition is a direct extension of the previous one to the current framework.Let Z and R be the primal and dual optimal sets of problems ( 33) and (34).The first partition is defined as follows, the second in a similar way: We can then define the analytic center of the optimal set R as follows.If R is reduced to a singleton, then it is this point, otherwise it is the minimum of the problem min r∈ri( R) The analytic center can be characterized like in Proposition 4.4, by which it can be shown that Theorem 4.5 still holds for the rolling friction framework.
Algorithm 1 is modified in order to solve a RF problem.This is described by Algorithm 2.
The linear system solved at each iteration is obtained by linearizing the system (35).It is reformulated under the form of the following augmented system where Z = Arw(z), R = Arw(J ⊤ r) and µ = r ⊤ Jz n .As in Algorithm 1, the affine scaling direction (d v a , d r a , d z a ) is the solution of (38) without the square bracketed term in the right-hand side, while the full step (d v , d r , d z ) is the solution of the complete system.
The scaling strategy is similar to that described in Section 5.The NT direction p is defined by the formula (25) where u and r are respectively replaced by z and J ⊤ r.The change of variables is done by setting ẑ := Q p z and y := Q p −1 J ⊤ r.
Recall that ẑ = y, which allows to symmetrize the linear system under the form Because of the ill-conditioning of the matrix Q p 2 , an equivalent form of (39) has been considered: where dz = Q p d z .Figure 3 shows the condition number of the three matrices ( 38)-( 40) along the iterations of the numerical resolution of a RF problem.It can be seen that the system ( 40) is better conditioned than (39).We also tried several reductions to a 2 × 2 form as in ( 28) or (30), leading to matrices of the form where Ĥ = P −1 H and P P ⊤ = JQ p −2 J ⊤ .We also tried several ways to compute the matrix P , by performing a Cholesky factorization or by directly exploiting the structure of a quadratic representation matrix.Despite such a reduction and in contrast to the results obtained with the RF problems, the numerical performance has not been improved.Moreover, the computation of the matrix P increases the overall computational cost, without any real improvement.We also observe the same for the 1 × 1 system like (31).The numerical tests were carried out on 526 RF problems of the FCLIB family (Acary et al., 2014), whose characteristics are in Table 4.The numerical results and performances of Algorithm 2 with the three linear systems described previously, are reported in Table 5 and Figure 4.These results show that, with a tolerance tol = 10 −10 , the choice of system (40) gives the best performance.But the performance gap between the systems with the NT scaling is smaller than those observed for the FC problems.It can also be observed that, as for the RF problems, without scaling the failures are of type 1, while with NT scaling the failures only occur when a NaN is detected.Moreover, in the latter case the maximum value of the residual (32) (right column of Table 5), shows that the stopping point of the algorithm is nearly optimal.
Figure 4: Performance profiles for the three linear systems for the RF model primal-dual algorithm, we were not satisfied with its efficiency and robustness.The Nesterov-Todd scaling strategy is a wonderful theoretical tool, but numerically very painful.As the iterates approach the boundary of the second order cones, the conditioning of the linear system explodes, the iterations get stuck on the boundary and divisions by zero occur, producing NaN and thus an emergency stop.We have therefore examined a large number of equivalent formulations of the linear system and found that the one which gives the best results, or shall we say the least bad, is the one in which the quadratic representation matrix which allows scaling is not a direct component of the matrix of the linear system.In addition, special attention must be paid to the way in which the matrix-vector products are performed in order to construct the system to be factorized.Unfortunately, we have not been able to find a formulation for rolling friction problems that is as efficient and robust as for friction cones.The accuracy we have achieved is slightly lower.Nevertheless, the accuracies and computation times we have achieved for both models seem to us to be quite adequate and can be used for real applications..The next step in this research work is to extend the primal-dual algorithm to solve the original model (1).A natural idea is to solve a sequence of systems parametrized by µ > 0: The main issues are that this system is not the optimality system of an optimization problem and that there is a non-smooth equation.
One theoretical question remains from this study.It concerns the characterisation of the limit point of the central path, like the formula (17) which defines the analytic center, but without the assumption of strict complementarity.Wriggers, P. (2006).Computational contact mechanics.Vol. 2. Springer Berlin Heidelberg.[doi], [oa].

A Euclidean Jordan algebra
Let us consider the set K = n i=1 K i where K i is an n i -dimensional Lorentz cone defined by Let N = n i=1 n i .For x ∈ R N , we denote x = (x 1 , . . ., x n ), where x i = (x i0 , xi ).For two vectors x and y in R N , the Jordan product is defined by , for i = 1, . . ., n.
Let x ∈ R N and x 2 = x • x.A fundamental property of the Jordan algebra for interior-point algorithms, is that the Lorentz cone is the cone of squares, that is K = {x 2 : x ∈ R N }, see (Alizadeh and Goldfarb, 2003, pp. 18-19).
For matrices A and B, we define For i ∈ {1, . . ., n}, let det(x i ) = x 2 i0 −∥x i ∥ 2 be the determinant of x i ∈ R ni .If x i is nonsingular, i.e., det(x i ) ̸ = 0, the inverse of x i is the unique vector of R ni such that x i • x −1 i = e i and is given by x −1 i = det(x i ) −1 R ni x i , where R ni is the reflexion matrix defined by If for all i ∈ {1, . . ., n} x i is nonsingular, we have For i ∈ {1, . . ., n}, the spectral decomposition of a vector x i ∈ R ni is defined by x i = λ i c i + λ ′ i c ′ i , where The scalars λ i and λ ′ i are eigenvalues of Arw(x i ), with the corresponding eigenvectors c i and c ′ i .This pair of eigenvectors is called the Jordan frame of x i .It is said that two vectors commute if they share the same Jordan frame.In that case the corresponding arrow matrices commute.From these definitions, it follows that x i ∈ K i (resp.x i ∈ int(K i )) if and only if λ ′ i ≥ 0 (resp.λ ′ i > 0).If x i is nonsingular, then More generally, for a continuous function f we can define At last, the quadratic representation of a vector x i ∈ K i is defined as The scalars λ 2 i and λ ′ 2 i are eigenvalues of Q xi with the eigenvectors c i and c ′ i .In particular, two vectors commute if and only if their quadratic representation matrices commute.For x ∈ K, we define Q x = Q x1 ⊕ . . .⊕ Q xn .For i ∈ {1, . . ., n} and x ∈ int(K i ), we have ∇ x (log det(x i )) = 2x −1 i and ∇ 2 x (log det(x i )) = −2Q x −1 i .These operators Arw and Q • are useful for the design of interior algorithms.See (Alizadeh and Goldfarb, 2003, §4) for a review of their interesting properties.

B Proof of Proposition 3.2
Before giving a proof of Proposition 3.2, we propose some equivalent formulations of Assumption 3.1.They are proved thanks to the following well-known lemma (Skarpness and Sposito, 1982, Lemma 2), called Tucker's theorem of the alternative in the case where K is the nonnegative orthant, see (Mangasarian, 1969, p. 29).The proof can be made by applying a separation theorem of convex sets, see, e.g., (Rockafellar, 1970, Theorem 11.3).
Lemma B.1.Let K ⊂ R n be a closed pointed convex cone, A ∈ R m×n and B ∈ R p×n .One and only one of the following statements is true.
(i) There exists a non-zero x ∈ K such that Ax = 0 and Bx ≤ 0.
The following lemma provides equivalent formulations of Assumption 3.1.
Lemma B.2.Let F be the product of second order cones of the form (5) or (7).Let M ∈ R m×m be symmetric and positive-definite, f ∈ R m , H ∈ R nd×m , w ∈ R nd and W = HM −1 H ⊤ .The following four assertions are equivalent.
(i) There exists v ∈ R m such that Hv + w ∈ int(F * ).
(iii) There does not exist a nonzero d ∈ F, such that H ⊤ d = 0 and w ⊤ d ≤ 0.
(iv) There does not exist a nonzero d ∈ F, such that W d = 0 and q ⊤ d ≤ 0.
Proof.The implication (i) ⇒ (ii) is obvious.The equivalence (ii) ⇔ (iii) is a direct consequence of the fact that F * is a closed pointed convex cone and of Lemma B.1.The equivalence (iii) ⇔ (iv) follows from the definitions of W and q in (10), and of the positive definiteness of M .It remains to prove that (ii) ⇒ (i).Let (v, t) ∈ R m × R + such that Hv + tw ∈ int(F * ).If t > 0, set v ′ = 1 t v. Since Hv + tw = t(Hv ′ + w) and int(F * ) is a cone, Hv ′ + w ∈ int(F * ).Suppose now that t = 0.There exists ϵ > 0 such that Hv + B ϵ ⊂ F * , where B ϵ is the open ball centered at 0 with radius ϵ.Since F * is a cone, for all t > 0, tHv + B tϵ ⊂ F * .Let us choose t > 1 ϵ ∥w∥.We then have tHv + w ∈ int(F * ), which proves (i).We can prove now that the Slater's assumption is equivalent to the non-emptiness and compactness of the optimal set of the reduced problem (11).
Proof of Proposition 3.2.The assertions (i) and (ii) are direct consequences of the weak duality and of the strong convexity of the objective function of (9).The outcome (iii) can be proved by means of some useful tools from asymptotic analysis in convex optimization (Auslender and Teboulle, 2003).Let us define the function

Figure 1 :
Figure 1: Condition number (κ) of the matrix of the linear systems along the iterations when solving two different problems.The left figure is related to a problem with only one contact point, (nd, m) = (3, 3), while the right figure is with sixty-nine contact points, (nd, m) = (207, 450).

Figure 3 :
Figure 3: Condition number (κ) of the matrix of the linear systems along the iterations of the Algorithm 2 when solving the problem PrimitiveSoup-ndof-6000-nc-37-0 with n = 37 contact points.

Table 1 :
Sizes of problems of the FCLIB collectionAlgorithm 1 One iteration of Mehrotra primal-dual algorithm for solving a FC problem Parameters:

Table 2 :
Number of successes and failures when solving the FCLIB problems with tol = 10 −10 .

Table 4 :
Sizes of rolling friction problems of the FCLIB collection

Table 5 :
Number of successes and failures when solving the RF problems with tol = 10 −10