High-accuracy Computation of Rolling Friction Contact Problems

Our goal is to numerically solve optimization problems derived from a mechanical model of unilateral contact between solid bodies with rolling friction. The model is an optimization problem with a strictly convex quadratic objective function and a second-order cone of constraints that is not self-dual. The solver is an implementation of a primal-dual interior-point algorithm with the predictor-corrector scheme of Mehrotra extended to the second-order cone problem. We focused on analyzing the limits of numerical computation and proposed some treatments to achieve optimal solutions with ten significant digits of precision.


I. INTRODUCTION
Acary and Bourrier presented in [1] a general model for the problem of unilateral contact between solid bodies with rolling friction.Let n, m ∈ N be the number of contact points and degrees of freedom.Let d = 5 be the dimension of the velocity and reaction vectors at a contact point.In this paper, we consider only a convex relaxation of this model as where M ∈ R m×m is symmetric and positive-definite, f ∈ R m , H ∈ R nd×m , w ∈ R nd .The vector r ∈ R nd is the reaction force at contact, u ∈ R nd is the velocity of the contact point and v ∈ R m is the velocity of the system.The rolling friction cone is defined by The dual cone is then This cone is not self-dual since F = F * , a significant difference with the frictional contact (no rolling) framework [3].The system (1) coincides with the first order optimality conditions of the primaldual pair of convex problems: After the introduction of artificial variables and a reformulation of the constraints, we will show that these problems fall in the class of second-order cone optimization (SOCO) problems.Thus, a primal-dual interior-point method (IPM) is an appropriate choice to efficiently solve them.This paper proposes an adapted IPM for the solution of (2), in order to accurately solve (1).We do it by means of a C/C++ implementation of the algorithm in SICONOS 1 and perform tests on FCLIB [3], a collection of discrete rolling frictional contact problems 2 .The paper is organized as follows.Section II shows the rolling friction contact model and the system of equations solved by means of the interior point algorithm described in Section III.The numerical difficulties are analyzed in Section IV.The robustness of our implementation is highlighted by the numerical results presented in Section V.The paper ends with a conclusion in Section VI and some basics on Euclidean Jordan algebra as well as the notation in Appendix.

II. ROLLING FRICTION CONTACT MODEL
The cone appearing in a standard SOCO problem, usually the Lorentz cone, is self-dual.This allows to take advantage of the Jordan algebra.But for the rolling friction problems, this is not the case, the cone of constraints is not self-dual.There is no Jordan product such that R * i is a cone of squares, an important property to extend IPM from linear optimization to SOCO [2].To deal with this difficulty, we propose to transfom the conical constraints of Problem (2) into the product of Lorentz cones by means of artificial variables.Knowing that, for a triplet (a, b, c) of real numbers, a ≥ b + c if and only if there exists t, t ∈ R such that t ≥ b, t ≥ c and a = t + t .We apply this statement by setting ui,0 = ti + t i for all i ∈ {1, . . ., n}.Then, the primal-dual Problem (2) can be reformulated as where z ∈ R n×(d+1) is a vector with n components of the form (ti; ūi; t i ; ũi), i ∈ {1, . . ., n} and By virtue of [2, Lemma 15], the orthogonality condition of (1) becomes z • J r = 0, with z ∈ (L 2n ) * = L 2n and J r ∈ L 2n .This representation of (1) provides a square system of equalities.
It is assumed that the Slater hypothesis holds, i.e., there exists v ∈ R m such that Hv + w ∈ int(F * ).Under this assumption, (r, u) ∈ R 2nd is a primal-dual optimal solution of Problem (2) if and only if there exists (t, t ) ∈ R 2n such that the Karush-Kuhn-Tucker (KKT) are satisfied: In practice, the Newton method is applied to a perturbation of (4), such that the solution is maintained in the interior of the cone of constraints.The perturbation is to replace the complementarity conditions by z • J r = 2µe, where the perturbation parameter µ is the duality measure defined by The perturbed KKT system associated to Problem (3) is then It can be shown that the system (5) has a unique solution (z(µ), J r(µ)) ∈ int(L 4n ) for all µ > 0. The corresponding curve, parametrized by µ, is called the central path.It can be shown that under the Slater hypothesis, when µ ↓ 0, the central path converges to a relative interior point of the set of primal-dual solutions of (3).

III. PRIMAL-DUAL INTERIOR-POINT APPROACH
Basically, an interior point algorithm consists of successive linearizations of (5), while driving µ to zero and keeping the iterates z and J r in the interior of the Lorentz cones.The linearization of equations ( 5) leads to the following linear system where Z = Arw(z) and R = Arw(J r).Unfortunately, this system can be singular, even if the matrices Z and R are positive definite, see [8].This is due to the fact that the arrow matrices do not commute, except if they share a common set of eigenvectors.To overcome this drawback, the idea is to make a local change of variables, which keeps invariant the Lorentz cone as well as its interior and so that, in the new space, the arrow matrices commute.An extensive literature has discussed several choices of variable changes and it is commonly agreed that the best one for solving a SOCO problem is the Nesterov and Todd (NT) scaling scheme [7].It is defined as follows.Let p, p ∈ R nd be the NT directions defined with n conic components of the following form: for i ∈ {1, . . ., n}, Let us define the block-diagonal matrix Qp = n i=1 (Qp i ⊕ Qp i ).The change of variables is then ẑ := Qpz and y := Q p −1 J r.
A fundamental property of the NT scaling is that ẑ = y, see [2, p.42], which implies that the corresponding arrow matrices commute.
After the change of variables, we obtain the corresponding primaldual problems where Note that KJ = I.After some obvious simplifications, the new pertubed KKT system associated to ( 8) is given by The linearization with respect to the variables v, r and z of this new system leads to the following linear system: where Ẑ := Arw(ẑ) = Arw(y) =: Y .By left-multiplying the third equation by Qp Ẑ−1 and after some simplications, we get the following equivalent symmetric linear system: Algorithm 1 is the implemented algorithm that was used to perform our numerical tests.This is the primal-dual interior point method of Mehrotra [6], extended to the solution of a SOCO problem [7].At each iteration of Algorithm 1, two linear systems are solved, with two different right-hand-sides (RHS).The first one corresponds to a pure Newton direction, called the affine scaling direction and denoted by (d v a , d r a , d z a ).It satisfies the linear system (6) (resp.(11)) with µ = 0.The second one corresponds to the computation of a linear combination of the affine scaling direction and a correction step, that is used to force the iterates to get close to the central path.See [5, §14.2] for details.As a result, when solving (6), the third equation is replaced by the following one: where σ ∈ (0, 1) is the centralization parameter, i.e., the desired reduction factor of the duality measure.When solving (11), this equation becomes where dz a = Qpd z a and ďy a = Q p −1 J d r a .Once a linear system is solved, a line search is performed so that the next iterate remains within the cone of constraints.Let x ∈ L 4n be the current iterate and let d ∈ L 4n be the corresponding solution of the linear system that was solved.The function get_step_length(x, d) returns the largest α ≥ 0 such that x + αd ∈ L 4n .The step length along the affine scaling direction is the largest value α ∈ (0, 1] such that x + αd ∈ L 4n .This value is used to compute the centralization parameter σ.Next, the step length along the full direction is the largest value α ∈ (0, 1] such that x + αd ∈ (1 − γ)x + L 4n .See the detailed formulas for this computation in [7].At the end of the iteration, the variables are updated at Step 26.

26:
Set (v, r, z) = (v, r, z) + α(d v , d r , d z ) 27: end while IV.COMPUTATIONAL ISSUES Although the system (6) may be singular, we have performed numerical experiments with it.The singularity of the matrix is a rather rare event in practice.These first experiments serve as a kind of reference for future experiments, both in terms of number of iterations and CPU time.Since the matrix in ( 6) is not symmetric, an LU factorization is performed.The advantage of this "no-scaling" strategy is that it is simple to implement and computationally fast, because there is not a bunch of calculations due to the scaling.However, since an LU factorization is more expensive than a symmetric matrix factorization, this advantage should be lowered for large problems.
Our first experiments with the NT scaling scheme have been disastrous.This is due to the ill-conditioning of the matrix Qp near an optimal solution.This behavior is well known in the community of SOCO and previous experiments with the software SDPT3 [7] warned us about it.We have been searching for a long time for a formulation of the linear system that allows us to obtain suitable results.The form (11) seems to us not too bad, but not enough robust.We propose an alternative equivalent 3 × 3 symmetric linear system, by left-multiplying the last equation of (11) by Q p −1 and setting dz = Qpd z : The matrix of ( 14) is less ill-conditioned than the one of (11).
Although the ill-conditioned matrix Qp is introduced in the second equation, which leads to a loss of accuracy in the primal feasibility, the use of ( 14) allowed to solve 97% of the problems with a tolerance = 10 −10 .
Since the matrices in ( 11) and ( 14) are symmetric, an LDL factorization is done by means of the Harwell subroutine MA57 [9], an efficient code for the solution of sparse symmetric systems.Furthermore, in order to avoid the propagation of computational errors due to the ill-conditioning of Qp, all calculations of the form Qxy are performed without explicitly building the vector x and the matrix Qx.This is done, for example, to compute the vectors y and ẑ in the RHS of ( 14).We do it by using the formula Qxy = 2(x y)x − (det x)Ry and the spectral decompositions of the vectors x and y.In addition, to improve the accuracy, all these calculations are done with the floating-point data type long double instead of double.Indeed, a difficulty due to the computational limits of processors led us to use this type of data.In experiments with tight tolerances, when the algorithm returns an error, it is usually stopped by a computation involving a NaN (Not a Number).This comes from an invalid arithmetic operation, like a division by zero or the calculation of the square root of a negative number.Despite the fact that z and J r are kept in the interior of the Lorentz cones, when approaching the boundary of the constraints, the second eigenvalue (λ2) of some conic components becomes almost zero or even equal to zero.This further justifies our choice to use the data type long double.However, the use of the long double data type may produce a negative λ2.This is due to the fact that the algorithm is executed on 64-bit processors whose epsilon machine is in double precision, i.e., 16 digit precision [10].
We also tried some experiments with a solution of the linear system in a reduced form, which might have been preferred because of its smaller dimension.By eliminating the variable d z in the system (11), a reduced version is given by But, like as (11), the matrix of ( 15) is very ill-conditioned.Another reduction with a smaller condition number is given by where Ĥ = P −1 H and d r = P d r .The matrix P is such that JQ p −2 J = P P .This factorization can be done in several ways, such as Cholesky, LDU, etc. Whatever the formulation used, in addition to the numerical computation problems mentioned above, we again encounter numerous computations to perform at each iteration, which leads to propagating too many errors.Any reduction to a 2 × 2 form by eliminating the variable d z always returned worse results.

V. NUMERICAL RESULTS
Introduced in [3], three collections of 523 rolling friction problem tests have been carried out by Algorithm 1. Detailed descriptions of these problems are given in Table I.The performance profiles were introduced by Dolan and Moré [11] to get the benchmarking optimization solvers through a cumulative distribution function ρs(τ ).It is the probability that the performance ratio of the solver s to the best solver does not exceed a factor 2 τ .The higher the probability is, the better the solver is.Figures 1, 2 and  3 show the performance profiles of Algorithm 1, with a linear system under the form (6) "No scaling", (11) "Qp2" and ( 14) "JQinv", with tolerances = 10 −8 , = 10 −9 and = 10 −10 respectively.In Figure 1, the linear systems have been solved for the tolerance = 10 −8 .Within this tolerance, the three solvers are able to solve all problems successfully.The "no-scaling" system (6), represented by the blue curve, has shown it to be faster (better CPU time) but less efficient (larger number of iterations) than the other two.The curves in red and green, corresponding to (11) and (14), do not differ significantly in this assessment.However, when the tolerance is decreased to 10 −9 and 10 −10 , the difference becomes more significant between these curves.In Figure 2 and Figure 3, despite still being the fastest one, the blue curve has shown the worst efficiency as the number of results is only approximately 82% and 74% of the total number of experiments.Meanwhile, the red curves in Figure 2 and Figure 3 show that (11) has started to suffer from problems of NaN, so its result has stopped at 98% and 85%.Finally, through the greens, we can see that (14) still retains its robustness.There is no failure with the tolerances 10 −8 , 10 −9 .
When the tolerance is 10 −10 , the algorithm is stopped by NaN for only 14/523 tests, thus 97% performance.

VI. CONCLUSION
The reformation of the non self-dual cone of constraints of the original problem into a self-dual one of the SOCO problem is an important contribution to the proposed interior-point algorithm.First, a scaling strategy by NT directions to produce a symmetric linear system is performed to overcome the disadvantages that come from the asymmetric linear system.Then, because of the considerable computational volume of this strategy, it is essential to choose suitable alternative formulas, for example the calculation of Qxy.This not only creates efficiency in terms of execution time, but also is to achieve better accuracy of optimal solutions.Finally, the experiments indicate that Algorithm 1 is efficient and robust for solving the rolling friction problem to obtain the optimal solution of 10 digit precision.

APPENDIX
We will take a glance at the Euclidean Jordan algebra which is a well-known framework associated with SOCO.See the reference paper [2] for a detailed presentation.We consider column vectors of the form (x0; x) ∈ R × R N −1 .For each such vector x ∈ R N , we define an Arrow-shaped matrix given by For any arbitrary power p ∈ Q, we define x p = λ p 1 c1 + λ p 2 c2.For example, a useful formula for the gradient of the logarithm of the determinant is given by ).The quadratic-representation associated to x is defined by Qx = 2Arw 2 (x) − Arw(x 2 ) = 2xx − det(x)R.
For vectors x, y of the form (v1; . . .; v ), such that vi ∈ R n i belongs to the second order cone of dimension ni, we define
Note that x ∈ L := {v = (v0; v) ∈ R × R N −1 : v ≤ v0} (known as Lorentz cone) if and only if Arw(x) is positive definite.In SOCO, it is useful to introduce the identity vector eN = (1;0) ∈ R × R N −1 , the reflection matrix RN = 1 0 0 −I ∈ R N ×N ,the direct sum of matricesA ⊕ B := A 0 0 B and p i=1 Ai = A1 ⊕ . . .⊕ Ap.The Jordan product of two vectors x, y is given byx • y := x y x0 ȳ + y0 x = Arw(x)y = Arw(x)Arw(y)e.The Jordan product is commutative, but not associative.The spectral decomposition of any element of this Jordan algebra is given byx = λ1c1 + λ2c2,where λ1,2 = x0 ± x and c1,2 = and eigenvectors of x.The trace and determinant of x are defined by tr(x) = λ1 + λ2 and det(x) = λ1λ2.