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

Vincent Acary, Paul Armand, Hoang Minh Nguyen, Maksym Shpakovych

► To cite this version:

Vincent Acary, Paul Armand, Hoang Minh Nguyen, Maksym Shpakovych. Second order cone programming for frictional contact mechanics using interior point algorithm. Optimization Methods and Software, 2024, pp.1-31. 10.1080/10556788.2023.2296438. hal-03913568v3

HAL Id: hal-03913568

https://hal.science/hal-03913568v3

Submitted on 16 Jan 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Second order cone programming for frictional contact mechanics using interior point algorithm

Vincent Acary, Paul Armand, Hoang Minh Nguyen, Maksym Shpakovych

► To cite this version:

Vincent Acary, Paul Armand, Hoang Minh Nguyen, Maksym Shpakovych. Second order cone programming for frictional contact mechanics using interior point algorithm. Optimization Methods and Software, 2024, Special issue honoring Oleg Burdakov, 39 (3), pp.634-663. 10.1080/10556788.2023.2296438. hal-03913568v3

HAL Id: hal-03913568

https://hal.science/hal-03913568v3

Submitted on 16 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Second order cone programming for frictional contact mechanics
using interior point algorithm

Vincent Acary* Paul Armand† Hoang Minh Nguyen‡ Maksym Shpakovych§

January 16, 2024

Abstract

We report experiments of an implementation of a primal-dual interior point algorithm for solving me-
chanical models of one-sided contact problems with Coulomb friction. The objective is to recover an optimal
solution with high accuracy and as quickly as possible. These developments are part of the design of Siconos<fnref n="1" />,
an open-source software for modeling and simulating non-smooth dynamical systems. Currently, Siconos uses
mainly first order methods for the numerical solution of these systems. These methods are very robust, but
suffer from a linear rate of convergence and are therefore too slow to find accurate solutions in a reasonable
time. Our main objective is to apply second-order methods to speed up convergence. In this paper, we focus
on solving a relaxed formulation of the initial mechanical model, corresponding to the optimality conditions
of a convex quadratic minimization problem with second-order cone constraints. We will present in detail a
primal-dual interior point algorithm for solving this type of problem. The main difficulty in implementing
this algorithm arises from the fact that at each iteration of the algorithm, a change of variable, called a
scaling, has to be performed in order to guarantee the non-singularity of the linear system to be solved, as
well as to recover a symmetric system. Although this scaling strategy is very nice from a theoretical point
of view, it leads to an enormous deterioration of the conditioning of the linear system when approaching the
optimal solution and therefore to all the numerical difficulties that result from it. We will describe in detail
the numerical algebra that we have developed in our implementation, in order to overcome these problems of
numerical instability. We will also present the solution of the models resulting from the problems of rolling
friction, for which the cone of constraints is no longer self-dual like the Lorentz cone.

**keywords** Second order cone optimization, Interior-point algorithm, Contact mechanics, Coulomb friction

# 1 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., 1997, 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.

<footnote>*INRIA Rhône-Alpes; e-mail: vincent.acary@inria.fr</footnote>

<footnote>†Laboratoire XLIM, Université de Limoges; e-mail: paul.armand@unilim.fr</footnote>

<footnote>‡INRIA Rhône-Alpes; e-mail: hoang-minh.nguyen@inria.fr</footnote>

<footnote>§INRIA Bordeaux; e-mail: maksym.shpakovych@unilim.fr</footnote>

<footnote n="1">¹https://nonsmooth.gricad-pages.univ-grenoble-alpes.fr/siconos</footnote>
Second order cone programming for frictional contact mechanics
using interior point algorithm

Vincent Acary* Paul Armand† Hoang Minh Nguyen‡ Maksym Shpakovych§

January 16, 2024

Abstract

We report experiments of an implementation of a primal-dual interior point algorithm for solving me-
chanical models of one-sided contact problems with Coulomb friction. The objective is to recover an optimal
solution with high accuracy and as quickly as possible. These developments are part of the design of Siconos<fnref n="1" />,
an open-source software for modeling and simulating non-smooth dynamical systems. Currently, Siconos uses
mainly first order methods for the numerical solution of these systems. These methods are very robust, but
suffer from a linear rate of convergence and are therefore too slow to find accurate solutions in a reasonable
time. Our main objective is to apply second-order methods to speed up convergence. In this paper, we focus
on solving a relaxed formulation of the initial mechanical model, corresponding to the optimality conditions
of a convex quadratic minimization problem with second-order cone constraints. We will present in detail a
primal-dual interior point algorithm for solving this type of problem. The main difficulty in implementing
this algorithm arises from the fact that at each iteration of the algorithm, a change of variable, called a
scaling, has to be performed in order to guarantee the non-singularity of the linear system to be solved, as
well as to recover a symmetric system. Although this scaling strategy is very nice from a theoretical point
of view, it leads to an enormous deterioration of the conditioning of the linear system when approaching the
optimal solution and therefore to all the numerical difficulties that result from it. We will describe in detail
the numerical algebra that we have developed in our implementation, in order to overcome these problems of
numerical instability. We will also present the solution of the models resulting from the problems of rolling
friction, for which the cone of constraints is no longer self-dual like the Lorentz cone.

**keywords** Second order cone optimization, Interior-point algorithm, Contact mechanics, Coulomb friction

# 1 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., 1997, 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.

<footnote>*INRIA Rhône-Alpes; e-mail: vincent.acary@inria.fr</footnote>

<footnote>†Laboratoire XLIM, Université de Limoges; e-mail: paul.armand@unilim.fr</footnote>

<footnote>‡INRIA Rhône-Alpes; e-mail: hoang-minh.nguyen@inria.fr</footnote>

<footnote>§INRIA Bordeaux; e-mail: maksym.shpakovych@unilim.fr</footnote>

<footnote n="1">¹https://nonsmooth.gricad-pages.univ-grenoble-alpes.fr/siconos</footnote>
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 language 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
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 language 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 $\mathbb{R}^n$. The Euclidean scalar product of $x$ and $y$ is denoted $x^\top y$. The associated $\ell_2$ norm is $\|x\| = \sqrt{x^\top x}$. The perp notation $x \perp y$ means that $x^\top y = 0$. Let $E$ be a subset of $\mathbb{R}^n$. The interior, the relative interior and the boundary of $E$ are respectively denoted by $\text{int}(E)$, $\text{ri}(E)$ and $\text{bd}(E)$. The dual of $E$ is denoted and defined by $E^* = \{y \in \mathbb{R}^n : \forall x \in E, x^\top y \ge 0\}$. The indicator function of $E$ is $\delta_E(x) = 0$ if $x \in E$ and $\delta_E(x) = +\infty$ if $x \notin E$. The domain of an extended-valued function $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ is $\text{dom}f = \{x \in \mathbb{R}^n : f(x) \in \mathbb{R}\}$. The definitions and notations related to the Jordan algebra are described in Appendix A.

# 2 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 \in \mathbb{N}$ be the number of contact points and let $m \in \mathbb{N}$ be the number of degrees of freedom. The mechanical system is described by means of three vectors: the global velocity $v \in \mathbb{R}^m$, the velocity $u \in \mathbb{R}^{nd}$ and the reaction $r \in \mathbb{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

$$ 
\begin{aligned}
Mv + f &= H^\top r, \\
Hv + w &= u, \\
K^* &\ni u + \Phi(u) \perp r \in K,
\end{aligned} 
\quad (1) 
$$

where the data of the problem are $M \in \mathbb{R}^{m \times m}$ a symmetric and positive-definite matrix, $f \in \mathbb{R}^m$, $H \in \mathbb{R}^{nd \times m}$, $w \in \mathbb{R}^{nd}$, $K$ is a cone and $K^* := \{u \in \mathbb{R}^{nd} : u^\top r \ge 0, \text{ for all } r \in K\}$ is its dual cone. The cone is of the form $K = \prod_{i=1}^n K_i$, each $K_i$ being a cone whose definition depends on the model.

For a FC problem, we have

$$ K_i = \{r_i = (r_{i,N}, r_{i,T}^\top)^\top \in \mathbb{R} \times \mathbb{R}^{d-1} : \mu_i r_{i,N} \ge \|r_{i,T}\|\}. $$
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 $\mathbb{R}^n$. The Euclidean scalar product of $x$ and $y$ is denoted $x^\top y$. The associated $\ell_2$ norm is $\|x\| = \sqrt{x^\top x}$. The perp notation $x \perp y$ means that $x^\top y = 0$. Let $E$ be a subset of $\mathbb{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 $E^* = \{y \in \mathbb{R}^n : \forall x \in E, x^\top y \ge 0\}$. The indicator function of $E$ is $\delta_E(x) = 0$ if $x \in E$ and $\delta_E(x) = +\infty$ if $x \notin E$. The domain of an extended-valued function $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ is dom$f = \{x \in \mathbb{R}^n : f(x) \in \mathbb{R}\}$. The definitions and notations related to the Jordan algebra are described in Appendix A.

# 2 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 \in \mathbb{N}$ be the number of contact points and let $m \in \mathbb{N}$ be the number of degrees of freedom. The mechanical system is described by means of three vectors: the global velocity $v \in \mathbb{R}^m$, the velocity $u \in \mathbb{R}^{nd}$ and the reaction $r \in \mathbb{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

$$ 
\begin{aligned}
Mv + f &= H^\top r, \\
Hv + w &= u, \\
K^* &\ni u + \Phi(u) \perp r \in K,
\end{aligned} 
\quad (1) 
$$

where the data of the problem are $M \in \mathbb{R}^{m \times m}$ a symmetric and positive-definite matrix, $f \in \mathbb{R}^m$, $H \in \mathbb{R}^{nd \times m}$, $w \in \mathbb{R}^{nd}$, $K$ is a cone and $K^* := \{u \in \mathbb{R}^{nd} : u^\top r \ge 0, \text{ for all } r \in K\}$ is its dual cone. The cone is of the form $K = \prod_{i=1}^n K_i$, each $K_i$ being a cone whose definition depends on the model. For a FC problem, we have

$$ K_i = \{r_i = (r_{i,N}, r_{i,T})^\top \in \mathbb{R} \times \mathbb{R}^{d-1} : \mu_i r_{i,N} \ge \|r_{i,T}\|\}. $$
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 $\mu \in \mathbb{R}^n$ is the vector of friction coefficients. The dual cone of $\mathcal{K}$, is given by $\mathcal{K}^* = \prod_{i=1}^n \mathcal{K}_i^*$, where

$$ \mathcal{K}_i^* = \{u_i = (u_{i,N}, u_{i,T})^\top \in \mathbb{R} \times \mathbb{R}^{d-1} : u_{i,N} \ge \mu_i \|u_{i,T}\|\}. $$

The function $\Phi: \mathbb{R}^{nd} \to \mathbb{R}^{nd}$ is defined by $n$ components $\Phi_i(u) = (\mu_i \|u_{i,T}\|, 0^\top)^\top \in \mathbb{R}^d$.

For a RF problem, the cones are of the form

$$ \mathcal{K}_i = \{r_i = (r_{i,N}, r_{i,T}, r_{i,R})^\top \in \mathbb{R} \times \mathbb{R}^{d-1} \times \mathbb{R}^{d-1} : \mu_i r_{i,N} \ge \|r_{i,T}\|, \mu_{R,i} r_{i,N} \ge \|r_{i,R}\|\}. $$

The vector $r_{i,R}$ is the rolling friction reaction moment at contact and $\mu_R \in \mathbb{R}^n$ is the vector of rolling friction coefficients. The dual cone of $\mathcal{K}_i$ is

$$ \mathcal{K}_i^* = \{u_i = (u_{i,N}, u_{i,T}, u_{i,R})^\top \in \mathbb{R} \times \mathbb{R}^{d-1} \times \mathbb{R}^{d-1} : u_{i,N} \ge \mu_i \|u_{i,T}\| + \mu_{R,i} \|u_{i,R}\|\}. $$

The vector $u_{i,R}$ is the relative angular velocity at contact. For this model, the components of the function $\Phi$ are defined by $\Phi_i(u) = (\mu_i \|u_{i,T}\| + \mu_{R,i} \|u_{i,R}\|, 0^\top)^\top \in \mathbb{R}^d$.

By observing that $\Phi(u) = \Phi(u + \Phi(u))$ and by making the change of variable $u \leftarrow u + \Phi(u)$, the problem (1) is reformulated as

$$ \begin{aligned} & Mv + f = H^\top r, \\ & Hv + w + \Phi(u) = u, \\ & \mathcal{K}^* \ni u \perp r \in \mathcal{K}. \end{aligned} \quad (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

$$ Hv + w + \phi(s) = u \quad \text{and} \quad s_i = \mu_i \|u_{i,T}\|, \quad i \in \{1, \dots, n\}, $$

where $\phi: \mathbb{R}^n \to \mathbb{R}^{nd}$ is defined by $n$ components of the form $(s_i, 0^\top)^\top \in \mathbb{R}^d$. The idea is to solve the nonlinear system as a parametric system of the form

$$ \begin{aligned} & Mv + f = H^\top r, \\ & Hv + w + \phi(s) = u, \\ & \mathcal{K}^* \ni u \perp r \in \mathcal{K}, \end{aligned} \quad (3) $$

where $s$ is periodically updated according to

$$ s_i = \mu_i \|u_{i,T}\|, \quad i \in \{1, \dots, n\}. $$

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:

$$ \begin{array}{ll} \min_{v,u} & \frac{1}{2} v^\top M v + f^\top v \\ \text{s.t.} & Hv + w + \phi(s) = u, \\ & u \in \mathcal{K}^*. \end{array} $$

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.

## 3 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_\mu$ whose ith diagonal block is diag$(1, \mu_i, \mu_i)$ for the FC model and diag$(1, \mu_i, \mu_i, \mu_{R,i}, \mu_{R,i})$ for the RF model. The parameter $s$ is assumed to be fixed. Let us make the change of variables:

$$ u \leftarrow P_\mu u, \quad H \leftarrow P_\mu H, \quad w \leftarrow P_\mu(w + \phi(s)) \quad \text{and} \quad r \leftarrow P_\mu^{-1} r. $$
We keep the same names for the variables, but change the notation of the components. For the RF model we define

$$u_i = \begin{pmatrix} u_{i,0} \\ \bar{u}_i \\ \tilde{u}_i \end{pmatrix} = \begin{pmatrix} u_{i,N} \\ \mu_i u_{i,T} \\ \mu_{R,i} u_{i,R} \end{pmatrix} \quad \text{and} \quad r_i = \begin{pmatrix} r_{i,0} \\ \bar{r}_i \\ \tilde{r}_i \end{pmatrix} = \begin{pmatrix} r_{i,N} \\ \frac{1}{\mu_i} r_{i,T} \\ \frac{1}{\mu_{R,i}} r_{i,R} \end{pmatrix}.$$

The same notation is used for the FC model, except that the components $\tilde{u}_i$ and $\tilde{r}_i$ do not appear. The system (3) becomes

$$ \begin{aligned} Mv + f &= H^\top r, \\ Hv + w &= u, \\ \mathcal{F}^* &\ni u \perp r \in \mathcal{F}, \end{aligned} \quad (4) $$

where $\mathcal{F}$ is the product of cones. The definition of $\mathcal{F}$ depends on the model. For the FC model ($d=3$) we have

$$ \mathcal{F} = \prod_{i=1}^{n} \mathcal{L}_i := \{ r_i = (r_{i0}, \bar{r}_i) \in \mathbb{R}^d : r_{i0} \ge \|\bar{r}_i\| \}. \quad (5) $$

This cone is self-dual (i.e., $\mathcal{L}_i^* = \mathcal{L}_i$), therefore

$$ \mathcal{F}^* = \prod_{i=1}^{n} \mathcal{L}_i^* := \{ u_i = (u_{i0}, \bar{u}_i) \in \mathbb{R}^d : u_{i0} \ge \|\bar{u}_i\| \}. \quad (6) $$

For the RF model ($d=5$), the cone of constraints is

$$ \mathcal{F} = \prod_{i=1}^{n} \mathcal{R}_i := \{ r_i = (r_{i0}, \bar{r}_i, \tilde{r}_i) \in \mathbb{R}^d : r_{i0} \ge \max\{\|\bar{r}_i\|, \|\tilde{r}_i\|\} \}. \quad (7) $$

This cone is not self-dual. The dual cone is then

$$ \mathcal{F}^* = \prod_{i=1}^{n} \mathcal{R}_i^* := \{ u_i = (u_{i0}, \bar{u}_i, \tilde{u}_i) \in \mathbb{R}^d : u_{i0} \ge \|\bar{u}_i\| + \|\tilde{u}_i\| \}. \quad (8) $$

The system (4) can be viewed as the first order optimality conditions of the following second order cone optimization problem:

$$ \begin{aligned} & \min_{v,u} \quad \frac{1}{2} v^\top M v + f^\top v \\ & \text{s.t.} \quad \begin{array}{l} Hv + w = u, \\ u \in \mathcal{F}^*. \end{array} \end{aligned} \quad (9) $$

Let $r \in \mathbb{R}^{nd}$ be the dual variable associated to the linear constraint of (9). The Lagrangian function associated to problem (9) is defined by

$$ L(v, u, r) = \frac{1}{2} v^\top M v + f^\top v + r^\top (u - Hv - w). $$

Since $L$ is separable in $v$ and $u$, convex in $v$ and linear in $u$, and since $\mathcal{F}^{**} = \mathcal{F}$, the dual function is readily obtained and defined by

$$ \inf_{(v,u) \in \mathbb{R}^m \times \mathcal{F}^*} L(v, u, r) = \begin{cases} -\frac{1}{2} r^\top W r - q^\top r - \frac{1}{2} f^\top M^{-1} f & \text{if } r \in \mathcal{F}, \\ -\infty & \text{otherwise,} \end{cases} $$

where

$$ W = HM^{-1}H^{\top} \quad \text{and} \quad q = w - HM^{-1}f. \quad (10) $$

The dual problem is then

$$ \begin{aligned} & \min_{r} \quad \frac{1}{2} r^\top W r + q^\top r \\ & \text{s.t.} \quad r \in \mathcal{F}. \end{aligned} \quad (11) $$
We keep the same names for the variables, but change the notation of the components. For the RF model we define

$$u_i = \begin{pmatrix} u_{i,0} \\ \bar{u}_i \\ \tilde{u}_i \end{pmatrix} = \begin{pmatrix} u_{i,N} \\ \mu_i u_{i,T} \\ \mu_{R,i} u_{i,R} \end{pmatrix} \quad \text{and} \quad r_i = \begin{pmatrix} r_{i,0} \\ \bar{r}_i \\ \tilde{r}_i \end{pmatrix} = \begin{pmatrix} r_{i,N} \\ \frac{1}{\mu_i} r_{i,T} \\ \frac{1}{\mu_{R,i}} r_{i,R} \end{pmatrix}.$$

The same notation is used for the FC model, except that the components $\tilde{u}_i$ and $\tilde{r}_i$ do not appear. The system (3) becomes

$$ \begin{aligned} Mv + f &= H^\top r, \\ Hv + w &= u, \\ \mathcal{F}^* &\ni u \perp r \in \mathcal{F}, \end{aligned} \quad (4) $$

where $\mathcal{F}$ is the product of cones. The definition of $\mathcal{F}$ depends on the model. For the FC model ($d=3$) we have

$$ \mathcal{F} = \prod_{i=1}^{n} \mathcal{L}_i := \{r_i = (r_{i0}, \bar{r}_i) \in \mathbb{R}^d : r_{i0} \ge \|\bar{r}_i\|\}. \quad (5) $$

This cone is self-dual (i.e., $\mathcal{L}_i^* = \mathcal{L}_i$), therefore

$$ \mathcal{F}^* = \prod_{i=1}^{n} \mathcal{L}_i^* := \{u_i = (u_{i0}, \bar{u}_i) \in \mathbb{R}^d : u_{i0} \ge \|\bar{u}_i\|\}. \quad (6) $$

For the RF model ($d=5$), the cone of constraints is

$$ \mathcal{F} = \prod_{i=1}^{n} \mathcal{R}_i := \{r_i = (r_{i0}, \bar{r}_i, \tilde{r}_i) \in \mathbb{R}^d : r_{i0} \ge \max\{\|\bar{r}_i\|, \|\tilde{r}_i\|\}\}. \quad (7) $$

This cone is not self-dual. The dual cone is then

$$ \mathcal{F}^* = \prod_{i=1}^{n} \mathcal{R}_i^* := \{u_i = (u_{i0}, \bar{u}_i, \tilde{u}_i) \in \mathbb{R}^d : u_{i0} \ge \|\bar{u}_i\| + \|\tilde{u}_i\|\}. \quad (8) $$

The system (4) can be viewed as the first order optimality conditions of the following second order cone optimization problem:

$$ \begin{aligned} & \min_{v,u} \quad \frac{1}{2} v^\top M v + f^\top v \\ & \text{s.t.} \quad \begin{array}{l} Hv + w = u, \\ u \in \mathcal{F}^*. \end{array} \end{aligned} \quad (9) $$

Let $r \in \mathbb{R}^{nd}$ be the dual variable associated to the linear constraint of (9). The Lagrangian function associated to problem (9) is defined by

$$ L(v, u, r) = \frac{1}{2} v^\top M v + f^\top v + r^\top (u - Hv - w). $$

Since $L$ is separable in $v$ and $u$, convex in $v$ and linear in $u$, and since $\mathcal{F}^{**} = \mathcal{F}$, the dual function is readily obtained and defined by

$$ \inf_{(v,u) \in \mathbb{R}^m \times \mathcal{F}^*} L(v, u, r) = \begin{cases} -\frac{1}{2} r^\top W r - q^\top r - \frac{1}{2} f^\top M^{-1} f & \text{if } r \in \mathcal{F}, \\ -\infty & \text{otherwise,} \end{cases} $$

where

$$ W = HM^{-1}H^{\top} \quad \text{and} \quad q = w - HM^{-1}f. \quad (10) $$

The dual problem is then

$$ \begin{aligned} & \min_{r} \quad \frac{1}{2} r^\top W r + q^\top r \\ & \text{s.t.} \quad r \in \mathcal{F}. \end{aligned} \quad (11) $$
Let us denote by $\hat{R}$ the optimal set of this problem. It is a convex set, characterized by the first order optimality conditions:

$$ \mathcal{F}^* \ni W r + q \perp r \in \mathcal{F}. \qquad (12) $$

Note that these conditions are equivalent to (4), thanks to the definitions (10) and the equalities

$$ v = M^{-1}(H^{\top}r - f) \quad \text{and} \quad u = W r + q. \qquad (13) $$

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

$$ \frac{1}{2} v^{\top} M v + f^{\top} v - \left( -\frac{1}{2} r^{\top} W r - q^{\top} r - \frac{1}{2} f^{\top} M^{-1} f \right) = u^{\top} r \geq 0, $$

the inequality arising from the fact that $r \in \mathcal{F}$ and $u \in \mathcal{F}^*$. The strong duality and the existence of an optimal solution for (11), we need a constraints qualification assumption, known as the Slater hypothesis.

**Assumption 3.1.** There exists $v \in \mathbb{R}^m$ such that $Hv + w \in \text{int}(\mathcal{F}^*)$.

With respect to the reduced problem (11), the Slater hypothesis can be equivalently formulated as follows: for all nonzero $d \in \mathcal{F}$, if $Wd = 0$ then $q^{\top}d > 0$.

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 \in \mathbb{R}^{m \times m}$ is symmetric and positive-definite, $f \in \mathbb{R}^m$, $H \in \mathbb{R}^{nd \times m}$, $w \in \mathbb{R}^n$ and $\mathcal{F}$ is the product of second order cones of the form (5) or (7).

(i) Any solution of (4), provides a primal-dual optimal solution of (9)-(11).

(ii) If the problem (9) is feasible, then it has a unique optimal solution $(\hat{v}, \hat{u})$.

(iii) Assumption 3.1 is satisfied if and only if the optimal set of (11) is non-empty and compact.

# 4 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., $\mathcal{F} = \mathcal{L}^n$, where $\mathcal{L}$ is the Lorentz cone of dimension $d$. We recall that this cone is self-dual, i.e., $\mathcal{F}^* = \mathcal{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 \circ 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

$$ \begin{aligned} W r + q &= u, \\ r \circ u &= 0, \\ (r, u) &\in \mathcal{L}^{2n}. \end{aligned} \qquad (14) $$

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

$$ \begin{aligned} W r + q &= u, \\ u \circ r &= 2\mu e, \\ (u, r) &\in \text{int}(\mathcal{L}^{2n}), \end{aligned} \qquad (15) $$

where $e$ is the unit vector associated with the Jordan product and $\mu > 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 $\mu$ is called *barrier parameter*, because (15) can be interpreted as the optimality condition of the barrier problem

$$ \min_{r} \varphi_{\mu}(r) := \frac{1}{2} r^{\top} W r + q^{\top} r - \mu \sum_{k=1}^{n} \log \det r_k. \qquad (16) $$
Let us denote by $\hat{R}$ the optimal set of this problem. It is a convex set, characterized by the first order optimality conditions:

$$ \mathcal{F}^* \ni Wr + q \perp r \in \mathcal{F}. \qquad (12) $$

Note that these conditions are equivalent to (4), thanks to the definitions (10) and the equalities

$$ v = M^{-1}(H^\top r - f) \quad \text{and} \quad u = Wr + q. \qquad (13) $$

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

$$ \frac{1}{2} v^\top M v + f^\top v - \left(-\frac{1}{2} r^\top W r - q^\top r - \frac{1}{2} f^\top M^{-1} f\right) = u^\top r \geq 0, $$

the inequality arising from the fact that $r \in \mathcal{F}$ and $u \in \mathcal{F}^*$. The strong duality and the existence of an optimal solution for (11), we need a constraints qualification assumption, known as the Slater hypothesis.

**Assumption 3.1.** There exists $v \in \mathbb{R}^m$ such that $Hv + w \in \text{int}(\mathcal{F}^*)$.

With respect to the reduced problem (11), the Slater hypothesis can be equivalently formulated as follows: for all nonzero $d \in \mathcal{F}$, if $Wd = 0$ then $q^\top d > 0$.

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 \in \mathbb{R}^{m \times m}$ is symmetric and positive-definite, $f \in \mathbb{R}^m$, $H \in \mathbb{R}^{nd \times m}$, $w \in \mathbb{R}^n$ and $\mathcal{F}$ is the product of second order cones of the form (5) or (7).

(i) Any solution of (4), provides a primal-dual optimal solution of (9)-(11).

(ii) If the problem (9) is feasible, then it has a unique optimal solution $(\hat{v}, \hat{u})$.

(iii) Assumption 3.1 is satisfied if and only if the optimal set of (11) is non-empty and compact.

# 4 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., $\mathcal{F} = \mathcal{L}^n$, where $\mathcal{L}$ is the Lorentz cone of dimension $d$. We recall that this cone is self-dual, i.e., $\mathcal{F}^* = \mathcal{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 \circ 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

$$ \begin{aligned} W r + q &= u, \\ r \circ u &= 0, \\ (r, u) &\in \mathcal{L}^{2n}. \end{aligned} \qquad (14) $$

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

$$ \begin{aligned} W r + q &= u, \\ u \circ r &= 2\mu e, \\ (u, r) &\in \text{int}(\mathcal{L}^{2n}), \end{aligned} \qquad (15) $$

where $e$ is the unit vector associated with the Jordan product and $\mu > 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 $\mu$ is called *barrier parameter*, because (15) can be interpreted as the optimality condition of the barrier problem

$$ \min_r \varphi_\mu(r) := \frac{1}{2} r^\top W r + q^\top r - \mu \sum_{k=1}^n \log \det r_k. \qquad (16) $$
The first order optimality condition of (16) is $Wr + q - 2\mu r^{-1} = 0$. By introducing the variable $u = 2\mu 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 $\mu$. The proof is given in Appendix C

**Proposition 4.1.** Under Assumption 3.1, for all $\mu > 0$, the perturbed KKT system (15) has a unique solution $(r(\mu), u(\mu)) \in \text{int}(\mathcal{L}^{2n})$. For all $\bar{\mu} > 0$, the set $\{(r(\mu), u(\mu)) : 0 \le \mu \le \bar{\mu}\}$ is bounded.

From Proposition 3.2-(ii), the optimal solution of (9) is unique, therefore $\lim_{\mu \to 0} u(\mu) = \hat{u}$. For the curve $r(\cdot)$, it can be shown that $\lim_{\mu \to 0} r(\mu) = \hat{r}$, where $\hat{r} \in \text{ri}(\hat{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, \dots, n\}$ of the Lorentz cones is partitioned into six sets (Mohammad-Nezhad and Terlaky, 2021):

$$B = \{i : \exists r \in \hat{R}, r_i \in \text{int}(\mathcal{L})\}, N = \{i : \hat{u}_i \in \text{int}(\mathcal{L})\},$$

$$R = \{i : \hat{u}_i \in \text{bd}(\mathcal{L}) \setminus \{0\} \text{ and } \exists r \in \hat{R}, r_i \in \text{bd}(\mathcal{L}) \setminus \{0\}\},$$

$$T_0 = \{i \in T : \forall r \in \hat{R}, \hat{u}_i = r_i = 0\},$$

$$T_1 = \{i \in T : \hat{u}_i = 0 \text{ and } \exists r \in \hat{R}, r_i \in \text{bd}(\mathcal{L}) \setminus \{0\}\},$$

$$T_2 = \{i \in T : \hat{u}_i \in \text{bd}(\mathcal{L}) \setminus \{0\} \text{ and } \forall r \in \hat{R}, r_i = 0\}.$$

An optimal solution $r \in \hat{R}$ is maximally complementary if and only if for all $i \in B$, $r_i \in \text{int}(\mathcal{L})$ and for all $i \in \mathbb{R} \cup T_1$, $r_i \in \text{bd}(\mathcal{L}) \setminus \{0\}$.

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.

**Proposition 4.2.** Under Assumption 3.1, the central path $(r(\cdot), u(\cdot))$ converges to $(\hat{r}, \hat{u})$, where $\hat{r}$ is a maximally complementary optimal solution of (11).

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.

**Assumption 4.3.** There exists $r \in \hat{R}$ such that $\hat{u} + r \in \text{int}(\mathcal{L}^n)$.

This assumption is equivalent to $T_i = \emptyset$, for $i = 1, 2, 3$. In that case, $(B, N, \mathbb{R})$ is a partition of $\{1, \dots, n\}$. Under this assumption, we can define the analytic center of the optimal set $\hat{R}$. If $\hat{R}$ is a singleton, then the analytic center of $\hat{R}$ is this point. Otherwise, it is the unique optimal solution of the problem

$$\min_{r \in \text{ri}(\hat{R})} \psi(r) := - \sum_{i \in B} \log \det r_i - \sum_{i \in \mathbb{R}} \log r_{i,0}. \quad (17)$$

The next proposition shows that the analytic center is well defined.

**Proposition 4.4.** Under Assumptions 3.1 and 4.3, if $\hat{R}$ is not reduced to a singleton, then the problem (17) has a unique optimal solution $\hat{r} \in \text{ri}(\hat{R})$, characterized by the following property: for all $r \in \text{ri}(\hat{R})$,

$$\sum_{i \in B} r_i^\top \hat{r}_i^{-1} + \frac{1}{2} \sum_{i \in \mathbb{R}} \frac{r_{i,0}}{\hat{r}_{i,0}} \le |B| + \frac{1}{2} |\mathbb{R}|. \quad (18)$$

*Proof.* Let $r \in \text{ri}(\hat{R})$. For all $i \in B$, $\det r_i > 0$, and for all $i \in \mathbb{R}$, $r_{i,0} = \|\bar{r}_i\| > 0$. Since $\text{ri}(\hat{R}) \ne \emptyset$, one has $\text{dom}\psi \cap \text{ri}(\hat{R}) \ne \emptyset$. To show that problem (17) has at least one solution, it suffices to show that the function $\psi + \delta_{\text{ri}(\hat{R})}$ is coercive. This last property is a direct consequence of the compactness of the set $\hat{R}$.

The uniqueness of the minimum comes from the strict convexity of $\psi$ on $\text{ri}(\hat{R})$. Indeed, for a conic component $i \in B$, we have $\nabla_r^2(-\log \det r)_{r=r_i} = 2Q_{r_i^{-1}}$, which is positive definite for all $r_i \in \text{int}(\mathcal{L})$, see Appendix A. For
The first order optimality condition of (16) is $Wr + q - 2\mu r^{-1} = 0$. By introducing the variable $u = 2\mu 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 $\mu$. The proof is given in Appendix C

**Proposition 4.1.** Under Assumption 3.1, for all $\mu > 0$, the perturbed KKT system (15) has a unique solution $(r(\mu), u(\mu)) \in \text{int}(\mathcal{L}^{2n})$. For all $\bar{\mu} > 0$, the set $\{(r(\mu), u(\mu)) : 0 \le \mu \le \bar{\mu}\}$ is bounded.

From Proposition 3.2-(ii), the optimal solution of (9) is unique, therefore $\lim_{\mu \to 0} u(\mu) = \hat{u}$. For the curve $r(\cdot)$, it can be shown that $\lim_{\mu \to 0} r(\mu) = \hat{r}$, where $\hat{r} \in \text{ri}(\hat{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, \dots, n\}$ of the Lorentz cones is partitioned into six sets (Mohammad-Nezhad and Terlaky, 2021):

$$B = \{i : \exists r \in \hat{R}, r_i \in \text{int}(\mathcal{L})\}, N = \{i : \hat{u}_i \in \text{int}(\mathcal{L})\},$$

$$R = \{i : \hat{u}_i \in \text{bd}(\mathcal{L}) \setminus \{0\} \text{ and } \exists r \in \hat{R}, r_i \in \text{bd}(\mathcal{L}) \setminus \{0\}\},$$

$$T_0 = \{i \in T : \forall r \in \hat{R}, \hat{u}_i = r_i = 0\},$$

$$T_1 = \{i \in T : \hat{u}_i = 0 \text{ and } \exists r \in \hat{R}, r_i \in \text{bd}(\mathcal{L}) \setminus \{0\}\},$$

$$T_2 = \{i \in T : \hat{u}_i \in \text{bd}(\mathcal{L}) \setminus \{0\} \text{ and } \forall r \in \hat{R}, r_i = 0\}.$$

An optimal solution $r \in \hat{R}$ is maximally complementary if and only if for all $i \in B$, $r_i \in \text{int}(\mathcal{L})$ and for all $i \in \mathbb{R} \cup T_1$, $r_i \in \text{bd}(\mathcal{L}) \setminus \{0\}$.

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.

**Proposition 4.2.** Under Assumption 3.1, the central path $(r(\cdot), u(\cdot))$ converges to $(\hat{r}, \hat{u})$, where $\hat{r}$ is a maximally complementary optimal solution of (11).

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.

**Assumption 4.3.** There exists $r \in \hat{R}$ such that $\hat{u} + r \in \text{int}(\mathcal{L}^n)$.

This assumption is equivalent to $T_i = \emptyset$, for $i = 1, 2, 3$. In that case, $(B, N, \mathbb{R})$ is a partition of $\{1, \dots, n\}$. Under this assumption, we can define the analytic center of the optimal set $\hat{R}$. If $\hat{R}$ is a singleton, then the analytic center of $\hat{R}$ is this point. Otherwise, it is the unique optimal solution of the problem

$$\min_{r \in \text{ri}(\hat{R})} \psi(r) := - \sum_{i \in B} \log \det r_i - \sum_{i \in \mathbb{R}} \log r_{i,0}. \qquad (17)$$

The next proposition shows that the analytic center is well defined.

**Proposition 4.4.** Under Assumptions 3.1 and 4.3, if $\hat{R}$ is not reduced to a singleton, then the problem (17) has a unique optimal solution $\hat{r} \in \text{ri}(\hat{R})$, characterized by the following property: for all $r \in \text{ri}(\hat{R})$,

$$\sum_{i \in B} r_i^\top \hat{r}_i^{-1} + \frac{1}{2} \sum_{i \in \mathbb{R}} \frac{r_{i,0}}{\hat{r}_{i,0}} \le |B| + \frac{1}{2} |\mathbb{R}|. \qquad (18)$$

*Proof.* Let $r \in \text{ri}(\hat{R})$. For all $i \in B$, $\det r_i > 0$, and for all $i \in \mathbb{R}$, $r_{i,0} = \|\bar{r}_i\| > 0$. Since $\text{ri}(\hat{R}) \ne \emptyset$, one has $\text{dom}\psi \cap \text{ri}(\hat{R}) \ne \emptyset$. To show that problem (17) has at least one solution, it suffices to show that the function $\psi + \delta_{\text{ri}(\hat{R})}$ is coercive. This last property is a direct consequence of the compactness of the set $\hat{R}$.

The uniqueness of the minimum comes from the strict convexity of $\psi$ on $\text{ri}(\hat{R})$. Indeed, for a conic component $i \in B$, we have $\nabla_r^2(-\log \det r)_{r=r_i} = 2Q_{r_i^{-1}}$, which is positive definite for all $r_i \in \text{int}(\mathcal{L})$, see Appendix A. For
all conic components $i \in \mathbb{R}$, there exists $h_i \in \mathbb{R}^{d-1}$, with $\|h_i\| = 1$, such that for all $r \in \mathrm{ri}(\hat{R})$, we have $r_i = r_{i,0}(1, h_i^\top)^\top$. Indeed, suppose that for $i \in \mathbb{R}$, there exist $r$ and $r'$ in $\mathrm{ri}(\hat{R})$, such that $r_i$ and $r'_i$ are not collinear. By the triangle inequality, it is easy to see that $\frac{1}{2}(r_i + r'_i) \in \mathrm{int}(\mathcal{L})$, and thus $i \in \mathbb{B}$, which would contradict $i \in \mathbb{R}$. Finally, for all $r, r' \in \mathrm{ri}(\hat{R})$, if $r \neq r'$ then there exists $i \in \mathbb{B}$ such that $r_i \neq r'_i$ or there exists $i \in \mathbb{R}$ such that $r_{i,0} \neq r'_{i,0}$. In both cases, we have $\psi(r) \neq \psi(r')$, which implies that $\psi$ is strictly convex on $\mathrm{ri}(\hat{R})$.

Since $\psi$ is convex on $\mathrm{ri}(\hat{R})$, $\hat{r}$ is optimal if and only if $-\nabla\psi(\hat{r})$ is in the normal cone to $\hat{R}$ at $\hat{r}$, i.e., for all $r \in \mathrm{ri}(\hat{R})$, $\nabla\psi(\hat{r})^\top(r - \hat{r}) \ge 0$. For $r \in \mathrm{ri}(\hat{R})$, we have

$$ 
\begin{aligned}
\nabla\psi(\hat{r})^\top(r - \hat{r}) &= -2 \sum_{i \in \mathbb{B}} (\hat{r}_i^{-1})^\top (r_i - \hat{r}_i) - \sum_{i \in \mathbb{R}} \frac{r_{i,0} - \hat{r}_{i,0}}{\hat{r}_{i,0}} \\
&= -2 \sum_{i \in \mathbb{B}} r_i^\top \hat{r}_i^{-1} + 2|\mathbb{B}| - \sum_{i \in \mathbb{R}} \frac{r_{i,0}}{\hat{r}_{i,0}} + |\mathbb{R}|,
\end{aligned}
$$

from which we deduce that (18) is satisfied. $\square$

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(\cdot)$ converges to the analytic center of $\hat{R}$.

*Proof.* Suppose that $\hat{R}$ is not reduced to a singleton, otherwise the result is a direct consequence of Proposition 4.2. Let $r \in \mathrm{ri}(\hat{R})$. As in the proof of Proposition 4.2, for all $\mu > 0$, (44) is satisfied. By using $r(\mu) \circ u(\mu) = 2\mu e$ and the definition of the partition $(\mathbb{B}, \mathbb{N}, \mathbb{R})$, we have

$$ \sum_{i \in \mathbb{B}} r_i^\top r_i^{-1}(\mu) + \sum_{i \in \mathbb{N}} \hat{u}_i^\top u_i^{-1}(\mu) + \sum_{i \in \mathbb{R}} (r_i^\top r_i^{-1}(\mu) + \hat{u}_i^\top u_i^{-1}(\mu)) \le n. \quad (19) $$

For any index $i$, by using the Cauchy-Schwarz inequality and the fact that $r_{i,0} \ge \|\bar{r}_i\|$, we have

$$ 
\begin{aligned}
r_i^\top r_i^{-1}(\mu) &= \frac{r_i^\top R r_i(\mu)}{\det r_i(\mu)} \\
&= \frac{r_{i,0} r_{i,0}(\mu) - \bar{r}_i^\top \bar{r}_i(\mu)}{r_{i,0}^2(\mu) - \|\bar{r}_i(\mu)\|^2} \\
&\ge \frac{r_{i,0} r_{i,0}(\mu) - \|\bar{r}_i\| \|\bar{r}_i(\mu)\|}{(r_{i,0}(\mu) - \|\bar{r}_i(\mu)\|)(r_{i,0}(\mu) + \|\bar{r}_i(\mu)\|)} \\
&\ge \frac{r_{i,0}}{r_{i,0}(\mu) + \|\bar{r}_i(\mu)\|}.
\end{aligned}
\quad (20) $$

In the same manner, for all $i \in \{1, \dots, n\}$ we have

$$ \hat{u}_i^\top u_i^{-1}(\mu) \ge \frac{\hat{u}_{i,0}}{u_{i,0}(\mu) + \|\bar{u}_i(\mu)\|}. \quad (21) $$

From (19), (20) and (21), we deduce that

$$ \sum_{i \in \mathbb{B}} r_i^\top r_i^{-1}(\mu) + \sum_{i \in \mathbb{N}} \hat{u}_i^\top u_i^{-1}(\mu) + \sum_{i \in \mathbb{R}} \left( \frac{r_{i,0}}{r_{i,0}(\mu) + \|\bar{r}_i(\mu)\|} + \frac{\hat{u}_{i,0}}{u_{i,0}(\mu) + \|\bar{u}_i(\mu)\|} \right) \le n. $$

When $\mu$ tends to zero, the first sum tends to $\sum_{i \in \mathbb{B}} r_i^\top \hat{r}_i^{-1}$ and the second one to $|\mathbb{N}|$. By using the fact that $\hat{r}_{i,0} = \|\bar{r}_i\|$ and $\hat{u}_{i,0} = \|\bar{u}_i\|$ for $i \in \mathbb{R}$, each term of the third sum tends to $\frac{r_{i,0}}{2\hat{r}_{i,0}} + \frac{1}{2}$. By taking the limit as $\mu \to 0$ and using $n = |\mathbb{B}| + |\mathbb{N}| + |\mathbb{R}|$, we obtain

$$ \sum_{i \in \mathbb{B}} r_i^\top \hat{r}_i^{-1} + \frac{1}{2} \sum_{i \in \mathbb{R}} \frac{r_{i,0}}{\hat{r}_{i,0}} \le |\mathbb{B}| + \frac{1}{2} |\mathbb{R}|, $$

which shows by Proposition 4.4 that $\hat{r}$ is the optimal solution of (17), the analytic center of $\hat{R}$. $\square$
all conic components $i \in \mathbb{R}$, there exists $h_i \in \mathbb{R}^{d-1}$, with $\|h_i\| = 1$, such that for all $r \in \text{ri}(\hat{R})$, we have $r_i = r_{i,0}(1, h_i^\top)^\top$. Indeed, suppose that for $i \in \mathbb{R}$, there exist $r$ and $r'$ in $\text{ri}(\hat{R})$, such that $r_i$ and $r'_i$ are not collinear. By the triangle inequality, it is easy to see that $\frac{1}{2}(r_i + r'_i) \in \text{int}(\mathcal{L})$, and thus $i \in \mathbb{B}$, which would contradict $i \in \mathbb{R}$. Finally, for all $r, r' \in \text{ri}(\hat{R})$, if $r \neq r'$ then there exists $i \in \mathbb{B}$ such that $r_i \neq r'_i$ or there exists $i \in \mathbb{R}$ such that $r_{i,0} \neq r'_{i,0}$. In both cases, we have $\psi(r) \neq \psi(r')$, which implies that $\psi$ is strictly convex on $\text{ri}(\hat{R})$.

Since $\psi$ is convex on $\text{ri}(\hat{R})$, $\hat{r}$ is optimal if and only if $-\nabla\psi(\hat{r})$ is in the normal cone to $\hat{R}$ at $\hat{r}$, i.e., for all $r \in \text{ri}(\hat{R})$, $\nabla\psi(\hat{r})^\top(r - \hat{r}) \ge 0$. For $r \in \text{ri}(\hat{R})$, we have

$$ 
\begin{aligned}
\nabla\psi(\hat{r})^\top(r - \hat{r}) &= -2 \sum_{i \in \mathbb{B}} (\hat{r}_i^{-1})^\top (r_i - \hat{r}_i) - \sum_{i \in \mathbb{R}} \frac{r_{i,0} - \hat{r}_{i,0}}{\hat{r}_{i,0}} \\
&= -2 \sum_{i \in \mathbb{B}} r_i^\top \hat{r}_i^{-1} + 2|\mathbb{B}| - \sum_{i \in \mathbb{R}} \frac{r_{i,0}}{\hat{r}_{i,0}} + |\mathbb{R}|,
\end{aligned}
$$

from which we deduce that (18) is satisfied. $\square$

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(\cdot)$ converges to the analytic center of $\hat{R}$.

*Proof.* Suppose that $\hat{R}$ is not reduced to a singleton, otherwise the result is a direct consequence of Proposition 4.2. Let $r \in \text{ri}(\hat{R})$. As in the proof of Proposition 4.2, for all $\mu > 0$, (44) is satisfied. By using $r(\mu) \circ u(\mu) = 2\mu e$ and the definition of the partition $(\mathbb{B}, \mathbb{N}, \mathbb{R})$, we have

$$ \sum_{i \in \mathbb{B}} r_i^\top r_i^{-1}(\mu) + \sum_{i \in \mathbb{N}} \hat{u}_i^\top u_i^{-1}(\mu) + \sum_{i \in \mathbb{R}} (r_i^\top r_i^{-1}(\mu) + \hat{u}_i^\top u_i^{-1}(\mu)) \le n. \quad (19) $$

For any index $i$, by using the Cauchy-Schwarz inequality and the fact that $r_{i,0} \ge \|\bar{r}_i\|$, we have

$$ 
\begin{aligned}
r_i^\top r_i^{-1}(\mu) &= \frac{r_i^\top R r_i(\mu)}{\det r_i(\mu)} \\
&= \frac{r_{i,0} r_{i,0}(\mu) - \bar{r}_i^\top \bar{r}_i(\mu)}{r_{i,0}^2(\mu) - \|\bar{r}_i(\mu)\|^2} \\
&\ge \frac{r_{i,0} r_{i,0}(\mu) - \|\bar{r}_i\| \|\bar{r}_i(\mu)\|}{(r_{i,0}(\mu) - \|\bar{r}_i(\mu)\|)(r_{i,0}(\mu) + \|\bar{r}_i(\mu)\|)} \\
&\ge \frac{r_{i,0}}{r_{i,0}(\mu) + \|\bar{r}_i(\mu)\|}.
\end{aligned}
\quad (20) $$

In the same manner, for all $i \in \{1, \dots, n\}$ we have

$$ \hat{u}_i^\top u_i^{-1}(\mu) \ge \frac{\hat{u}_{i,0}}{u_{i,0}(\mu) + \|\bar{u}_i(\mu)\|}. \quad (21) $$

From (19), (20) and (21), we deduce that

$$ \sum_{i \in \mathbb{B}} r_i^\top r_i^{-1}(\mu) + \sum_{i \in \mathbb{N}} \hat{u}_i^\top u_i^{-1}(\mu) + \sum_{i \in \mathbb{R}} \left( \frac{r_{i,0}}{r_{i,0}(\mu) + \|\bar{r}_i(\mu)\|} + \frac{\hat{u}_{i,0}}{u_{i,0}(\mu) + \|\bar{u}_i(\mu)\|} \right) \le n. $$

When $\mu$ tends to zero, the first sum tends to $\sum_{i \in \mathbb{B}} r_i^\top \hat{r}_i^{-1}$ and the second one to $|\mathbb{N}|$. By using the fact that $\hat{r}_{i,0} = \|\bar{r}_i\|$ and $\hat{u}_{i,0} = \|\bar{u}_i\|$ for $i \in \mathbb{R}$, each term of the third sum tends to $\frac{r_{i,0}}{2\hat{r}_{i,0}} + \frac{1}{2}$. By taking the limit as $\mu \to 0$ and using $n = |\mathbb{B}| + |\mathbb{N}| + |\mathbb{R}|$, we obtain

$$ \sum_{i \in \mathbb{B}} r_i^\top \hat{r}_i^{-1} + \frac{1}{2} \sum_{i \in \mathbb{R}} \frac{r_{i,0}}{\hat{r}_{i,0}} \le |\mathbb{B}| + \frac{1}{2} |\mathbb{R}|, $$

which shows by Proposition 4.4 that $\hat{r}$ is the optimal solution of (17), the analytic center of $\hat{R}$. $\square$
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.

# 5 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 \times 2$ block linear system into a $3 \times 3$ block system. We also show, by numerical experiments, that it is better to reduce again this system back to a $2 \times 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) \in \text{int}(\mathcal{L}^{2n})$. They only differ on the right-hand side and are of the form:

$$ \begin{pmatrix} W & -I \\ U & R \end{pmatrix} \begin{pmatrix} d^r \\ d^u \end{pmatrix} = \begin{pmatrix} -Wr - q + u \\ -u \circ r - [d_a^u \circ d_a^r - 2\sigma\mu e] \end{pmatrix} \quad (22) $$

where $U = \text{Arw}(u)$ and $R = \text{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_a^u, d_a^r)$ 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

$$ u \circ d_a^r + r \circ d_a^u = -u \circ r. $$

The affine scaling direction is a Newton step on the original optimality system (14). The barrier parameter is set to $\mu = \frac{u^\top r}{n}$. The centralization parameter $\sigma \in (0, 1]$ is fixed by comparing the current value of $\mu$ with its expected reduction obtained along the affine step. The second direction is a linear combination of the affine
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.

# 5 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 \times 2$ block linear system into a $3 \times 3$ block system. We also show, by numerical experiments, that it is better to reduce again this system back to a $2 \times 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) \in \text{int}(\mathcal{L}^{2n})$. They only differ on the right-hand side and are of the form:

$$ \begin{pmatrix} W & -I \\ U & R \end{pmatrix} \begin{pmatrix} d^r \\ d^u \end{pmatrix} = \begin{pmatrix} -Wr - q + u \\ -u \circ r - [d_a^u \circ d_a^r - 2\sigma\mu e] \end{pmatrix} \quad (22) $$

where $U = \text{Arw}(u)$ and $R = \text{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_a^u, d_a^r)$ 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

$$ u \circ d_a^r + r \circ d_a^u = -u \circ r. $$

The affine scaling direction is a Newton step on the original optimality system (14). The barrier parameter is set to $\mu = \frac{u^\top r}{n}$. The centralization parameter $\sigma \in (0, 1]$ is fixed by comparing the current value of $\mu$ 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 \circ d^r + r \circ d^u = -u \circ r - d_a^u \circ d_a^r + 2\sigma \frac{u^\top r}{n} e.$$

The next iterate is set to $(u^+, r^+) = (u, r) + \alpha(d^u, d^r)$, where $\alpha$ is the largest value in $(0, 1]$ such that

$$(u^+, r^+) \in (1 - \tau)(u, r) + \mathcal{L}^{2n},$$

for some value $\tau \in (0, 1]$ (typically $\tau = 0.995$). Contrary to common practice, different primal and dual stelengths 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 $Wr+q-u = 0$ and that different stelengths $\alpha \neq \alpha'$ are taken with $d^u \neq 0$. We then have

$$W(r + \alpha d^r) + q - (u + \alpha' d^u) = (\alpha - \alpha')d^u.$$

If $d^u \neq 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): If $U = \text{Arw}([1, 0.7, 0.7]^{\top})$, $R = \text{Arw}([1, 0.8, 0.5]^{\top})$ and $W = \text{diag}([0.3, 1, 0]^{\top})$, then $\det(W + R^{-1}U) = 0$.

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 \in \text{int} \mathcal{K}$ and let $Q_p$ be the associated quadratic representation (see Appendix A). Since $p$ is in the interior of the cone $\mathcal{K}$, the matrix $Q_p$ is positive definite. From (Alizadeh and Goldfarb, 2003, Theorem 9), we have $Q_p(\mathcal{L}^n) = \mathcal{L}^n$ and $Q_p(\text{int} \mathcal{L}^n) = \text{int} \mathcal{L}^n$. Let us consider the change of variables

$$\tilde{r} = Q_{p-1}r \quad \text{and} \quad \hat{u} = Q_p u.$$

The problem (11) becomes

$$\begin{array}{ll} \displaystyle \min_{\tilde{r}} & \frac{1}{2} \tilde{r}^{\top} Q_p W Q_p \tilde{r} + (Q_p q)^{\top} \tilde{r} \\ \text{s.t.} & \tilde{r} \in \mathcal{L}^n. \end{array} \qquad (23)$$

The corresponding perturbed KKT system is

$$W r + q = u \quad \text{and} \quad \hat{u} \circ \tilde{r} = 2\mu e,$$

with $(\hat{u}, \tilde{r}) \in \text{int} \mathcal{L}^{2n}$. The linearization of these equations leads the following linear system:

$$\begin{pmatrix} W & -I \\ \hat{U} Q_{p-1} & \tilde{R} Q_p \end{pmatrix} \begin{pmatrix} d^r \\ d^u \end{pmatrix} = \begin{pmatrix} -W r - q + u \\ -\hat{u} \circ \tilde{r} - [\hat{d}_a^u \circ \tilde{d}_a^r - 2\sigma\mu e] \end{pmatrix}, \qquad (24)$$

where $\hat{U} = \text{Arw}(\hat{u})$ and $\tilde{R} = \text{Arw}(\tilde{r})$. The choice of the vector $p \in \text{int} \mathcal{K}$ is made so that $\hat{u}$ and $\tilde{r}$ commute, which implies that the matrix of the linear system (24) is nonsingular. Indeed, this matrix is nonsingular if and only if $\det(W + (\tilde{R}Q_p)^{-1}\hat{U}Q_{p-1}) \neq 0$. Since $\hat{U}$ and $\tilde{R}$ are positive definite, $\hat{u}$ and $\tilde{r}$ commute, and $Q_{p-1} = Q_p^{-1}$, we have

$$(\tilde{R}Q_p)^{-1}\hat{U}Q_{p-1} = Q_{p-1}\tilde{R}^{-1/2}\hat{U}^{1/2}\hat{U}^{1/2}\tilde{R}^{-1/2}Q_{p-1},$$

which is symmetric and positive definite.
scaling direction and of a corrector step, to keep the iterates near the central path. It then satisfies the following linear equation

$$u \circ d^r + r \circ d^u = -u \circ r - d_a^u \circ d_a^r + 2\sigma \frac{u^\top r}{n} e.$$

The next iterate is set to $(u^+, r^+) = (u, r) + \alpha(d^u, d^r)$, where $\alpha$ is the largest value in $(0, 1]$ such that

$$(u^+, r^+) \in (1 - \tau)(u, r) + \mathcal{L}^{2n},$$

for some value $\tau \in (0, 1]$ (typically $\tau = 0.995$). Contrary to common practice, different primal and dual stelengths 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 $Wr+q-u = 0$ and that different stelengths $\alpha \neq \alpha'$ are taken with $d^u \neq 0$. We then have

$$W(r + \alpha d^r) + q - (u + \alpha' d^u) = (\alpha - \alpha')d^u.$$

If $d^u \neq 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): If $U = \text{Arw}([1, 0.7, 0.7]^{\top})$, $R = \text{Arw}([1, 0.8, 0.5]^{\top})$ and $W = \text{diag}([0.3, 1, 0]^{\top})$, then $\det(W + R^{-1}U) = 0$.

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 \in \text{int} \mathcal{K}$ and let $Q_p$ be the associated quadratic representation (see Appendix A). Since $p$ is in the interior of the cone $\mathcal{K}$, the matrix $Q_p$ is positive definite. From (Alizadeh and Goldfarb, 2003, Theorem 9), we have $Q_p(\mathcal{L}^n) = \mathcal{L}^n$ and $Q_p(\text{int} \mathcal{L}^n) = \text{int} \mathcal{L}^n$. Let us consider the change of variables

$$\tilde{r} = Q_{p-1}r \quad \text{and} \quad \hat{u} = Q_p u.$$

The problem (11) becomes

$$\begin{array}{ll} \displaystyle \min_{\tilde{r}} & \frac{1}{2} \tilde{r}^{\top} Q_p W Q_p \tilde{r} + (Q_p q)^{\top} \tilde{r} \\ \text{s.t.} & \tilde{r} \in \mathcal{L}^n. \end{array} \qquad (23)$$

The corresponding perturbed KKT system is

$$W r + q = u \quad \text{and} \quad \hat{u} \circ \tilde{r} = 2\mu e,$$

with $(\hat{u}, \tilde{r}) \in \text{int} \mathcal{L}^{2n}$. The linearization of these equations leads the following linear system:

$$\left( \begin{array}{cc} W & -I \\ \hat{U} Q_{p-1} & \tilde{R} Q_p \end{array} \right) \left( \begin{array}{c} d^r \\ d^u \end{array} \right) = \left( \begin{array}{c} -W r - q + u \\ -\hat{u} \circ \tilde{r} - [\hat{d}_a^u \circ \tilde{d}_a^r - 2\sigma\mu e] \end{array} \right), \qquad (24)$$

where $\hat{U} = \text{Arw}(\hat{u})$ and $\tilde{R} = \text{Arw}(\tilde{r})$. The choice of the vector $p \in \text{int} \mathcal{K}$ is made so that $\hat{u}$ and $\tilde{r}$ commute, which implies that the matrix of the linear system (24) is nonsingular. Indeed, this matrix is nonsingular if and only if $\det(W + (\tilde{R}Q_p)^{-1}\hat{U}Q_{p-1}) \neq 0$. Since $\hat{U}$ and $\tilde{R}$ are positive definite, $\hat{u}$ and $\tilde{r}$ commute, and $Q_{p-1} = Q_p^{-1}$, we have

$$(\tilde{R}Q_p)^{-1}\hat{U}Q_{p-1} = Q_{p-1}\tilde{R}^{-1/2}\hat{U}^{1/2}\hat{U}^{1/2}\tilde{R}^{-1/2}Q_{p-1},$$

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):

$$p = (Q_{u^{1/2}}(Q_{u^{1/2}r})^{-1/2})^{-1/2} = (Q_{r^{-1/2}}(Q_{r^{1/2}}u)^{1/2})^{-1/2}. \quad (25)$$

The main property of the NT direction is that $\hat{u} = \tilde{r}$, which implies that

$$\hat{U}^{-1}\tilde{R} = I.$$

The symmetrization of the system (24) is done by left-multiplying the last row by $-Q_p\hat{U}^{-1}$, leading to a symmetric matrix of the form

$$\begin{pmatrix} W & -I \\ -I & -Q_{p^2} \end{pmatrix}.$$

To take advantage of the sparsity of the matrices $M$ and $H$, the system (22) is considered in the following equivalent augmented form:

$$\begin{pmatrix} M & -H^{\top} & 0 \\ -H & 0 & I \\ 0 & U & R \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ d^u \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ Hv + w - u \\ -u \circ r - [d_a^u \circ d_a^r - 2\sigma\mu e] \end{pmatrix}. \quad (26)$$

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

$$\begin{pmatrix} M & -H^{\top} & 0 \\ -H & 0 & I \\ 0 & I & Q_{p^2} \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ d^u \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ Hv + w - u \\ -r - [Q_{p^2}(\hat{u}^{-1} \circ (\hat{d}_a^u \circ \check{d}_a^r)) - 2\sigma\mu u^{-1}] \end{pmatrix}. \quad (27)$$

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

$$\begin{pmatrix} M & -H^{\top} \\ -H & -Q_{p^{-2}} \end{pmatrix} \begin{pmatrix} d^v \\ d^r \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ Hv + w + [Q_{p^{-1}}(\hat{u}^{-1} \circ (\hat{d}_a^u \circ \check{d}_a^r)) - 2\sigma\mu r^{-1}] \end{pmatrix}. \quad (28)$$

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 \in \{1, \dots, n\}$, three situations can occur (Cai and Toh, 2006):

* $u_i^* \in \text{int}(\mathcal{L})$ and $r_i^* = 0$, then all the eigenvalues of $Q_{p_i^2}$ are of order $\mu := u^\top r$;

* $u_i^* = 0$ and $r_i^* \in \text{int}(\mathcal{L})$, then all the eigenvalues of $Q_{p_i^2}$ are of order $1/\mu$;

* $u_i^* \in \text{bd}(\mathcal{L})$, $r_i^* \in \text{bd}(\mathcal{L})$ and $(u_i^*, r_i^*) \neq (0,0)$, then the largest eigenvalue of $Q_{p_i^2}$ is of order $1/\mu$ and the smallest is of order $\mu$.

To overcome the difficulties due to ill-conditioning, we propose to solve the linear systems (27) and (28) under the following equivalent form:

$$\begin{pmatrix} M & -\hat{H}^{\top} & 0 \\ -\hat{H} & 0 & I \\ 0 & I & I \end{pmatrix} \begin{pmatrix} d^v \\ \check{d}^r \\ \hat{d}^u \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ \hat{H}v + \hat{w} - \hat{u} \\ -\check{r} - [\hat{u}^{-1} \circ (\hat{d}_a^u \circ \check{d}_a^r) - 2\sigma\mu\hat{u}^{-1}] \end{pmatrix} \quad (29)$$

and

$$\begin{pmatrix} M & -\hat{H}^{\top} \\ -\hat{H} & -I \end{pmatrix} \begin{pmatrix} d^v \\ \check{d}^r \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ \hat{H}v + \hat{w} + [\hat{u}^{-1} \circ (\hat{d}_a^u \circ \check{d}_a^r) - 2\sigma\mu\hat{u}^{-1}] \end{pmatrix}, \quad (30)$$

where $\hat{H} = Q_p H$, $\hat{w} = Q_p w$, $\hat{d}^u = Q_p d^u$ and $\check{d}^r = Q_{p^{-1}} d^r$.
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):

$$p = (Q_{u^{1/2}}(Q_{u^{1/2}}r)^{-1/2})^{-1/2} = (Q_{r^{-1/2}}(Q_{r^{1/2}}u)^{1/2})^{-1/2}. \quad (25)$$

The main property of the NT direction is that $\hat{u} = \tilde{r}$, which implies that

$$\hat{U}^{-1} \tilde{R} = I.$$

The symmetrization of the system (24) is done by left-multiplying the last row by $-Q_p \hat{U}^{-1}$, leading to a symmetric matrix of the form

$$\begin{pmatrix} W & -I \\ -I & -Q_{p^2} \end{pmatrix}.$$

To take advantage of the sparsity of the matrices $M$ and $H$, the system (22) is considered in the following equivalent augmented form:

$$\begin{pmatrix} M & -H^{\top} & 0 \\ -H & 0 & I \\ 0 & U & R \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ d^u \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ Hv + w - u \\ -u \circ r - [d_a^u \circ d_a^r - 2\sigma\mu e] \end{pmatrix}. \quad (26)$$

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

$$\begin{pmatrix} M & -H^{\top} & 0 \\ -H & 0 & I \\ 0 & I & Q_{p^2} \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ d^u \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ Hv + w - u \\ -r - [Q_{p^2}(\hat{u}^{-1} \circ (\hat{d}_a^u \circ \hat{d}_a^r)) - 2\sigma\mu u^{-1}] \end{pmatrix}. \quad (27)$$

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

$$\begin{pmatrix} M & -H^{\top} \\ -H & -Q_{p^{-2}} \end{pmatrix} \begin{pmatrix} d^v \\ d^r \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ Hv + w + [Q_{p^{-1}}(\hat{u}^{-1} \circ (\hat{d}_a^u \circ \hat{d}_a^r)) - 2\sigma\mu r^{-1}] \end{pmatrix}. \quad (28)$$

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 \in \{1, \dots, n\}$, three situations can occur (Cai and Toh, 2006):

* $u_i^* \in \text{int}(\mathcal{L})$ and $r_i^* = 0$, then all the eigenvalues of $Q_{p_i^2}$ are of order $\mu := u^\top r$;

* $u_i^* = 0$ and $r_i^* \in \text{int}(\mathcal{L})$, then all the eigenvalues of $Q_{p_i^2}$ are of order $1/\mu$;

* $u_i^* \in \text{bd}(\mathcal{L})$, $r_i^* \in \text{bd}(\mathcal{L})$ and $(u_i^*, r_i^*) \neq (0,0)$, then the largest eigenvalue of $Q_{p_i^2}$ is of order $1/\mu$ and the smallest is of order $\mu$.

To overcome the difficulties due to ill-conditioning, we propose to solve the linear systems (27) and (28) under the following equivalent form:

$$\begin{pmatrix} M & -\hat{H}^{\top} & 0 \\ -\hat{H} & 0 & I \\ 0 & I & I \end{pmatrix} \begin{pmatrix} d^v \\ \check{d}^r \\ \hat{d}^u \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ \hat{H}v + \hat{w} - \hat{u} \\ -\check{r} - [\hat{u}^{-1} \circ (\hat{d}_a^u \circ \hat{d}_a^r) - 2\sigma\mu\hat{u}^{-1}] \end{pmatrix} \quad (29)$$

and

$$\begin{pmatrix} M & -\hat{H}^{\top} \\ -\hat{H} & -I \end{pmatrix} \begin{pmatrix} d^v \\ \check{d}^r \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ \hat{H}v + \hat{w} + [\hat{u}^{-1} \circ (\hat{d}_a^u \circ \hat{d}_a^r) - 2\sigma\mu\hat{u}^{-1}] \end{pmatrix}, \quad (30)$$

where $\hat{H} = Q_p H$, $\hat{w} = Q_p w$, $\hat{d}^u = Q_p d^u$ and $\hat{d}^r = Q_{p^{-1}} d^r$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Condition number ($\kappa$) 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).</figcaption></figure>

In our numerical experiments, we also consider the reduced system

$$ (\hat{H}M^{-1}\hat{H}^{\top} + I)\check{d}^r = -\hat{H}M^{-1}(f + H^{\top}r) - \hat{w} - [\hat{u}^{-1} \circ (\hat{d}_a^u \circ \check{d}_a^r) - 2\sigma\mu\hat{u}^{-1}], \quad (31) $$

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: $n = 1, m = 3, M = I, w = 0, f = [3, 3, 3, 1, -1, -3, 1, -1, -3]^{\top}, H = [D, 0, -D]$ where $D = \text{diag}(1, 0.1, 0.1)$. 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)<fnref n="2" />. There are $n = 69$ contact points and $m = 450$ degrees of freedom, $H \in \mathbb{R}^{207 \times 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 $\sqrt{\epsilon}$, where $\epsilon$ 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) \in \mathcal{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^\top y)x - (\det x)R_d y$. Therefore, the product of a vector by a matrix $Q_p$, where $p$ is the NT-vector (25), 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.

<footnote n="2">²https://frictionalcontactlibrary.github.io</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Condition number ($\kappa$) 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).</figcaption></figure>

In our numerical experiments, we also consider the reduced system

$$ (\hat{H}M^{-1}\hat{H}^{\top} + I)\check{d}^r = -\hat{H}M^{-1}(f + H^{\top}r) - \hat{w} - [\hat{u}^{-1} \circ (\hat{d}_a^u \circ \check{d}_a^r) - 2\sigma\mu\hat{u}^{-1}], \quad (31) $$

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: $n = 1, m = 3, M = I, w = 0, f = [3, 3, 3, 1, -1, -3, 1, -1, -3]^{\top}, H = [D, 0, -D]$ where $D = \text{diag}(1, 0.1, 0.1)$. 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)<fnref n="2" />. There are $n = 69$ contact points and $m = 450$ degrees of freedom, $H \in \mathbb{R}^{207 \times 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 $\sqrt{\epsilon}$, where $\epsilon$ 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) \in \mathcal{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^\top y)x - (\det x)R_d y$. Therefore, the product of a vector by a matrix $Q_p$, where $p$ is the NT-vector (25), 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.

<footnote n="2">²https://frictionalcontactlibrary.github.io</footnote>
<table><caption>Table 1: Sizes of problems of the FCLIB collection</caption><thead><tr><th>Family</th><th># Problems</th><th>n</th><th>m</th></tr></thead><tbody><tr><td>BoxStacks</td><td>28</td><td>[31, 557]</td><td>450</td></tr><tr><td>Capsules</td><td>200</td><td>[15, 314]</td><td>600</td></tr><tr><td>Chute</td><td>182</td><td>[3, 3224]</td><td>[1002, 12672]</td></tr><tr><td>KaplasTower</td><td>240</td><td>[48, 899]</td><td>792</td></tr><tr><td>PrimitveSoup</td><td>200</td><td>[37, 3123]</td><td>6000</td></tr><tr><td>Spheres</td><td>200</td><td>[344, 4290]</td><td>12000</td></tr><tr><td>Spheres1mm</td><td>41</td><td>[715, 4213]</td><td>12000</td></tr></tbody></table>

**Algorithm 1** One iteration of Mehrotra primal-dual algorithm for solving a FC problem

**Parameters:** $\eta_1 \in (0, 1)$, $\eta_2 \ge 1$, $\eta_3 \ge 1$, $\tau_1 \in (0, 1)$, $\tau_2 \in (0, 1 - \tau_1)$, tol > 0

1: **if** the stopping criterion (32) is satisfied **then** return $(v, u, r)$ as primal-dual solution of (9);
2: Set $\mu \leftarrow u^\top r/n$;
3: Compute $(d_a^v, d_a^r, d_a^u)$ solution of (26) (resp. (27)-(31)) without the square bracketed term;
4: Find the greatest $\alpha_a \in (0, 1]$ such that $(u, r) + \alpha_a(d_a^u, d_a^r) \in \mathcal{L}^{2n}$;
5: Set $\mu_a \leftarrow (u + \alpha_a d_a^u)^\top (r + \alpha_a d_a^r)/n$;
6: **if** $\mu > \eta_1$ **then** set $e \leftarrow \max\{1, \eta_2 \alpha_a^2\}$ **else** set $e \leftarrow \eta_3$;
7: Set $\sigma \leftarrow \min\{1, (\mu_a/\mu)^e\}$;
8: Compute $(d_a^v, d_a^r, d_a^u)$ solution of (26) (resp. (27)-(31));
9: Set $\tau \leftarrow \tau_1 + \alpha_a \tau_2$;
10: Find the greatest $\alpha \in (0, 1]$ such that $(u, r) + \alpha(d_a^u, d_a^r) \in (1 - \tau)(u, r) + \mathcal{L}^{2n}$;
11: Set $(v, u, r) \leftarrow (v, u, r) + \alpha(d_a^v, d_a^u, d_a^r)$ and goto 1.

Algorithm 1 details one iteration of the implemented algorithm based on that of SDPT3 (Tütüncü et al., 2003). The stopping criterion used in GFC3D is

$$ \text{res} := \max \left\{ \frac{\|Hv + w - u\|}{\max\{\|Hv\|, \|w\|, \|u\|\}}, \frac{\|Mv + f - H^\top r\|}{\max\{\|Mv\|, \|f\|, \|H^\top r\|\}}, |u^\top r| \right\} \le \text{tol}, \quad (32) $$

The values of the parameters are fixed to $\eta_1 = 10^{-10}$, $\eta_2 = 3$, $\eta_3 = 1$, $\tau_1 = 0.9$, $\tau_2 = 0.09$ and tol = $10^{-10}$. The starting point is set as follows: for all $i \in \{1, \dots, n\}$, $u_i = r_i = (0.1, 0.01, 0.01)^\top$ and $v = M^{-1}(H^\top r + f)$. The experiments were done on 1091 problems of the FCLIB collection (Acary et al., 2014). There are seven families of problems, whose dimensions are described in Table 1, where $n$ is the number of three-dimensional cones and $m$ is six times the number of bodies.

The experiments were done on a MacBook Pro 2.3GHz equipped with a processor Intel Core i9 8 cœurs and 32Gb of memory, running OSX 10.15.7.

The linear system (26) is solved by means of a LU factorization, the symmetric ones with a $LDL^\top$ 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_{p2}$ 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
<table><caption>Table 1: Sizes of problems of the FCLIB collection</caption><thead><tr><th>Family</th><th># Problems</th><th>n</th><th>m</th></tr></thead><tbody><tr><td>BoxStacks</td><td>28</td><td>[31, 557]</td><td>450</td></tr><tr><td>Capsules</td><td>200</td><td>[15, 314]</td><td>600</td></tr><tr><td>Chute</td><td>182</td><td>[3, 3224]</td><td>[1002, 12672]</td></tr><tr><td>KaplasTower</td><td>240</td><td>[48, 899]</td><td>792</td></tr><tr><td>PrimitveSoup</td><td>200</td><td>[37, 3123]</td><td>6000</td></tr><tr><td>Spheres</td><td>200</td><td>[344, 4290]</td><td>12000</td></tr><tr><td>Spheres1mm</td><td>41</td><td>[715, 4213]</td><td>12000</td></tr></tbody></table>

**Algorithm 1** One iteration of Mehrotra primal-dual algorithm for solving a FC problem

**Parameters:** $\eta_1 \in (0, 1)$, $\eta_2 \ge 1$, $\eta_3 \ge 1$, $\tau_1 \in (0, 1)$, $\tau_2 \in (0, 1 - \tau_1)$, tol > 0

1: **if** the stopping criterion (32) is satisfied **then** return $(v, u, r)$ as primal-dual solution of (9);
2: Set $\mu \leftarrow u^\top r/n$;
3: Compute $(d_a^v, d_a^r, d_a^u)$ solution of (26) (resp. (27)-(31)) without the square bracketed term;
4: Find the greatest $\alpha_a \in (0, 1]$ such that $(u, r) + \alpha_a(d_a^u, d_a^r) \in \mathcal{L}^{2n}$;
5: Set $\mu_a \leftarrow (u + \alpha_a d_a^u)^\top (r + \alpha_a d_a^r)/n$;
6: **if** $\mu > \eta_1$ **then** set $e \leftarrow \max\{1, \eta_2 \alpha_a^2\}$ **else** set $e \leftarrow \eta_3$;
7: Set $\sigma \leftarrow \min\{1, (\mu_a/\mu)^e\}$;
8: Compute $(d_a^v, d_a^r, d_a^u)$ solution of (26) (resp. (27)-(31));
9: Set $\tau \leftarrow \tau_1 + \alpha_a \tau_2$;
10: Find the greatest $\alpha \in (0, 1]$ such that $(u, r) + \alpha(d_a^u, d_a^r) \in (1 - \tau)(u, r) + \mathcal{L}^{2n}$;
11: Set $(v, u, r) \leftarrow (v, u, r) + \alpha(d_a^v, d_a^u, d_a^r)$ and goto 1.

Algorithm 1 details one iteration of the implemented algorithm based on that of SDPT3 (Tütüncü et al., 2003). The stopping criterion used in GFC3D is

$$ \text{res} := \max \left\{ \frac{\|Hv + w - u\|}{\max\{\|Hv\|, \|w\|, \|u\|\}}, \frac{\|Mv + f - H^\top r\|}{\max\{\|Mv\|, \|f\|, \|H^\top r\|\}}, |u^\top r| \right\} \le \text{tol}, \quad (32) $$

The values of the parameters are fixed to $\eta_1 = 10^{-10}$, $\eta_2 = 3$, $\eta_3 = 1$, $\tau_1 = 0.9$, $\tau_2 = 0.09$ and tol = $10^{-10}$. The starting point is set as follows: for all $i \in \{1, \dots, n\}$, $u_i = r_i = (0.1, 0.01, 0.01)^\top$ and $v = M^{-1}(H^\top r + f)$. The experiments were done on 1091 problems of the FCLIB collection (Acary et al., 2014). There are seven families of problems, whose dimensions are described in Table 1, where $n$ is the number of three-dimensional cones and $m$ is six times the number of bodies.

The experiments were done on a MacBook Pro 2.3GHz equipped with a processor Intel Core i9 8 cœurs and 32Gb of memory, running OSX 10.15.7.

The linear system (26) is solved by means of a LU factorization, the symmetric ones with a $LDL^\top$ 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_{p2}$ 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.

<table><caption>Table 2: Number of successes and failures when solving the FCLIB problems with tol = $10^{-10}$.</caption><thead><tr><th>Linear system</th><th># Success</th><th># Failure 1</th><th># Failure 2</th></tr></thead><tbody><tr><td>(26)</td><td>1069</td><td>22</td><td>0</td></tr><tr><td>(27)</td><td>947</td><td>0</td><td>144</td></tr><tr><td>(28)</td><td>852</td><td>5</td><td>234</td></tr><tr><td>(29)</td><td>308</td><td>1</td><td>782</td></tr><tr><td>(29) + refin</td><td>1085</td><td>0</td><td>6</td></tr><tr><td>(30)</td><td>1091</td><td>0</td><td>0</td></tr><tr><td>(31)</td><td>1083</td><td>1</td><td>7</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Performance profiles for seven different linear system choices. For $\tau \ge 0$, $\rho_s(\tau)$ is the fraction of problems for which the performance of a given version of the algorithm is within a factor $2^\tau$ of the best one.</figcaption></figure>

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 `maxres` 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`.
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.

<table><caption>Table 2: Number of successes and failures when solving the FCLIB problems with tol = $10^{-10}$.</caption><thead><tr><th>Linear system</th><th># Success</th><th># Failure 1</th><th># Failure 2</th></tr></thead><tbody><tr><td>(26)</td><td>1069</td><td>22</td><td>0</td></tr><tr><td>(27)</td><td>947</td><td>0</td><td>144</td></tr><tr><td>(28)</td><td>852</td><td>5</td><td>234</td></tr><tr><td>(29)</td><td>308</td><td>1</td><td>782</td></tr><tr><td>(29) + refin</td><td>1085</td><td>0</td><td>6</td></tr><tr><td>(30)</td><td>1091</td><td>0</td><td>0</td></tr><tr><td>(31)</td><td>1083</td><td>1</td><td>7</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Performance profiles for seven different linear system choices. For $\tau \ge 0$, $\rho_s(\tau)$ is the fraction of problems for which the performance of a given version of the algorithm is within a factor $2^\tau$ of the best one.</figcaption></figure>

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 maxres 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><caption>Table 3: Performance comparison double versus long double when solving the FCLIB problems with system (30) and tol = 10⁻¹¹.</caption><thead><tr><th>Data type</th><th># Success</th><th># Failure 1</th><th># Failure 2</th><th>max res</th><th>cpu (sec)</th></tr></thead><tbody><tr><td>double</td><td>1047</td><td>5</td><td>39</td><td>9e-9</td><td>616</td></tr><tr><td>long double</td><td>1069</td><td>5</td><td>17</td><td>6e-11</td><td>675</td></tr></tbody></table>

# 6 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 $\mathcal{R}_i$ is no more self-dual and therefore not symmetric. There is no Jordan product such that $\mathcal{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 \ge b + c$ if and only if there exist two real numbers $t \ge b$ and $t' \ge c$ such that $a = t + t'$. By setting $u_{i,0} = \bar{t}_i + \tilde{t}_i$ for all $i \in \{1, \dots, n\}$, the primal-dual pair of problems (9)-(11) can be rewritten as

$$ 
\begin{array}{ll}
\min_{v,z} & \frac{1}{2} v^\top M v + f^\top v \\
\text{s.t.} & Hv + w = Jz, \\
& z \in \mathcal{L}^{2n},
\end{array}
\qquad (33) $$

and

$$ 
\begin{array}{ll}
\min_{r} & \frac{1}{2} r^{\top} W r + q^{\top} r \\
\text{s.t.} & J^{\top} r \in \mathcal{L}^{2n},
\end{array}
\qquad (34) $$

where $z = (z_1^\top, \dots, z_n^\top)^\top$ and $z_i = (\bar{t}_i, \tilde{u}_i^\top, \tilde{t}_i, \tilde{u}_i^\top)^\top \in \mathbb{R}^{d+1}$, $J \in \mathbb{R}^{nd \times n(d+1)}$ is a block diagonal matrix with $n$ blocks of the form

$$ J = \begin{pmatrix} 1 & 1 \\ I & \\ & I \end{pmatrix} \in \mathbb{R}^{d \times (d+1)}. $$

For all $i \in \{1, \dots, n\}$, we have $jz_i = (\bar{t}_i + \tilde{t}_i, \bar{u}_i^\top, \tilde{u}_i^\top)^\top$ and $j^\top r_i = (r_{i0}, \bar{r}_i^\top, r_{i0}, \tilde{r}_i^\top)^\top$.

The perturbed KKT system associated with the problem (34) is then

$$ 
\begin{array}{l}
Wr + q = Jz, \\
z \circ J^\top r = 2\mu e, \\
(z, J^\top r) \in \text{int}(\mathcal{L}^{4n}).
\end{array}
\qquad (35) $$

This optimality system is associated with the barrier problem

$$ \min_{r} \varphi_{\mu}(r) := \frac{1}{2} r^{\top} W r + q^{\top} r - \mu \sum_{i=1}^{n} \log(\det(r_{i,0}, \bar{r}_{i}) \det(r_{i,0}, \tilde{r}_{i})). \qquad (36) $$

As for Proposition 4.1 and by coercivity of the barrier function, under Assumption 3.1, for all $\mu > 0$ the system (35) has a unique solution such that $(z(\mu), J^\top r(\mu)) \in \text{int}(\mathcal{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(\cdot)$ 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.

**Assumption 6.1.** There exists a solution $(z, r) \in \mathbb{R}^{n(d+1)} \times \mathbb{R}^{nd}$ of (35) with $\mu = 0$, such that $z + J^\top r \in \text{int}(\mathcal{L}^{2n})$.

The strict complementarity hypothesis can be equivalently reformulated relatively to the original problems (9) and (11), thanks to the orthogonality condition (12). Let $(z, r)$ be a solution of (35) with $\mu = 0$. Then $(v, u = Jz)$ is the unique optimal solution of (9). Assumption 6.1 is satisfied if and only if for all $i \in \{1, \dots, n\}$, one of the following three assertions is satisfied:
<table><caption>Table 3: Performance comparison double versus long double when solving the FCLIB problems with system (30) and tol = 10⁻¹¹.</caption><thead><tr><th>Data type</th><th># Success</th><th># Failure 1</th><th># Failure 2</th><th>max res</th><th>cpu (sec)</th></tr></thead><tbody><tr><td>double</td><td>1047</td><td>5</td><td>39</td><td>9e-9</td><td>616</td></tr><tr><td>long double</td><td>1069</td><td>5</td><td>17</td><td>6e-11</td><td>675</td></tr></tbody></table>

# 6 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 $\mathcal{R}_i$ is no more self-dual and therefore not symmetric. There is no Jordan product such that $\mathcal{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 \ge b + c$ if and only if there exist two real numbers $t \ge b$ and $t' \ge c$ such that $a = t + t'$. By setting $u_{i,0} = \bar{t}_i + \tilde{t}_i$ for all $i \in \{1, \dots, n\}$, the primal-dual pair of problems (9)-(11) can be rewritten as

$$ 
\begin{array}{ll}
\min_{v,z} & \frac{1}{2} v^\top M v + f^\top v \\
\text{s.t.} & Hv + w = Jz, \\
& z \in \mathcal{L}^{2n},
\end{array}
\qquad (33) $$

and

$$ 
\begin{array}{ll}
\min_{r} & \frac{1}{2} r^{\top} W r + q^{\top} r \\
\text{s.t.} & J^{\top} r \in \mathcal{L}^{2n},
\end{array}
\qquad (34) $$

where $z = (z_1^\top, \dots, z_n^\top)^\top$ and $z_i = (\bar{t}_i, \tilde{u}_i^\top, \tilde{t}_i, \tilde{u}_i^\top)^\top \in \mathbb{R}^{d+1}$, $J \in \mathbb{R}^{nd \times n(d+1)}$ is a block diagonal matrix with $n$ blocks of the form

$$ J = \begin{pmatrix} 1 & 1 \\ I & \\ & I \end{pmatrix} \in \mathbb{R}^{d \times (d+1)}. $$

For all $i \in \{1, \dots, n\}$, we have $jz_i = (\bar{t}_i + \tilde{t}_i, \bar{u}_i^\top, \tilde{u}_i^\top)^\top$ and $j^\top r_i = (r_{i0}, \bar{r}_i^\top, r_{i0}, \tilde{r}_i^\top)^\top$.

The perturbed KKT system associated with the problem (34) is then

$$ 
\begin{array}{l}
Wr + q = Jz, \\
z \circ J^\top r = 2\mu e, \\
(z, J^\top r) \in \text{int}(\mathcal{L}^{4n}).
\end{array}
\qquad (35) $$

This optimality system is associated with the barrier problem

$$ \min_{r} \varphi_{\mu}(r) := \frac{1}{2} r^{\top} W r + q^{\top} r - \mu \sum_{i=1}^{n} \log(\det(r_{i,0}, \bar{r}_{i}) \det(r_{i,0}, \tilde{r}_{i})). \qquad (36) $$

As for Proposition 4.1 and by coercivity of the barrier function, under Assumption 3.1, for all $\mu > 0$ the system (35) has a unique solution such that $(z(\mu), J^\top r(\mu)) \in \text{int}(\mathcal{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(\cdot)$ 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.

**Assumption 6.1.** There exists a solution $(z, r) \in \mathbb{R}^{n(d+1)} \times \mathbb{R}^{nd}$ of (35) with $\mu = 0$, such that $z + J^\top r \in \text{int}(\mathcal{L}^{2n})$.

The strict complementarity hypothesis can be equivalently reformulated relatively to the original problems (9) and (11), thanks to the orthogonality condition (12). Let $(z, r)$ be a solution of (35) with $\mu = 0$. Then $(v, u = Jz)$ is the unique optimal solution of (9). Assumption 6.1 is satisfied if and only if for all $i \in \{1, \dots, n\}$, one of the following three assertions is satisfied:
* $r_{i0} > \max\{\|\tilde{r}_i\|, \|\tilde{r}_i\|\}$ and $u_i = 0$,

* $u_{i0} > \|\tilde{u}_i\| + \|\tilde{u}_i\|$ and $r_i = 0$,

* $r_{i0} = \max\{\|\tilde{r}_i\|, \|\tilde{r}_i\|\} > 0$ and $u_{i0} = \|\tilde{u}_i\| + \|\tilde{u}_i\| > 0$.

As in Section 4, under the strict complementarity assumption, the index set $\{1, \dots, n\}$ can be partitioned into two partitions $(\tilde{B}, \tilde{N}, \tilde{R})$ and $(\hat{B}, \hat{N}, \hat{R})$, whose definition is a direct extension of the previous one to the current framework. Let $\tilde{Z}$ and $\tilde{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:

$$ \tilde{B} = \{i : \exists r \in \tilde{R}, (r_{i,0}, \tilde{r}_i) \in \text{int}(\mathcal{L})\}, \quad \tilde{N} = \{i : \exists z \in \tilde{Z}, (\tilde{t}_i, \tilde{u}_i) \in \text{int}(\mathcal{L})\}, $$

$$ \tilde{R} = \{i : \exists (z, r) \in \tilde{Z} \times \tilde{R}, ((\tilde{t}_i, \tilde{u}_i), (r_{i,0}, \tilde{r}_i)) \in \text{bd}(\mathcal{L}^2) \setminus \{0\}\}. $$

We can then define the analytic center of the optimal set $\tilde{R}$ as follows. If $\tilde{R}$ is reduced to a singleton, then it is this point, otherwise it is the minimum of the problem

$$ \min_{r \in \text{ri}(\tilde{R})} - \sum_{i \in \tilde{B}} \log \det(r_{i,0}, \tilde{r}_i) - \sum_{i \in \tilde{B}} \log \det(r_{i,0}, \tilde{r}_i) - \sum_{i \in \tilde{R}, i \in \tilde{R}} \log r_{i,0}. \quad (37) $$

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.

<caption>**Algorithm 2** One iteration of Mehrotra primal-dual algorithm for solving a RF problem</caption>

**Parameters:** $\eta_1 \in (0, 1)$, $\eta_2 \ge 1$, $\eta_3 \ge 1$, $\tau_1 \in (0, 1)$, $\tau_2 \in (0, 1 - \tau_1)$, $\text{tol} > 0$
1: **if** the stopping criterion (32) with $u = Jz$ is satisfied **then** return $(v, z, r)$ as primal-dual solution of (33);
2: Set $\mu \leftarrow r^\top Jz/\text{tol}$;
3: Compute $(d_a^v, d_a^r, d_a^z)$ solution of (38) (resp. (39)-(40)) without the square bracketed term;
4: Find the greatest $\alpha_a \in (0, 1]$ such that $(z, J^\top r) + \alpha_a(d_a^z, J^\top d_a^r) \in \mathcal{L}^{4n}$;
5: Set $\mu_a \leftarrow (r + \alpha_a d_a^r)^\top J(z + \alpha_a d_a^z)/n$;
6: **if** $\mu > \eta_1$ **then** set $e \leftarrow \max\{1, \eta_2 \alpha_a^2\}$ **else** set $e \leftarrow \eta_3$;
7: Set $\sigma \leftarrow \min\{1, (\mu_a/\mu)^e\}$;
8: Compute $(d_a^v, d_a^r, d_a^z)$ solution of (38) (resp. (39)-(40));
9: Set $\tau \leftarrow \tau_1 + \alpha_a \tau_2$;
10: Find the greatest $\alpha \in (0, 1]$ such that $(z, J^\top r) + \alpha(d_a^z, J^\top d_a^r) \in (1 - \tau)(u, r) + \mathcal{L}^{4n}$;
11: Set $(v, z, r) \leftarrow (v, z, r) + \alpha(d_a^v, d_a^z, d_a^r)$ and goto 1.

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

$$ \begin{pmatrix} M & -H^\top & 0 \\ -H & 0 & J \\ 0 & ZJ^\top & R \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ d^z \end{pmatrix} = \begin{pmatrix} -Mv - f + H^\top r \\ Hv + w - Jz \\ -z \circ J^\top r - [d_a^z \circ (J^\top d_a^r) - 2\sigma\mu e] \end{pmatrix}, \quad (38) $$

where $Z = \text{Arw}(z)$, $R = \text{Arw}(J^\top r)$ and $\mu = \frac{r^\top Jz}{n}$. As in Algorithm 1, the affine scaling direction $(d_a^v, d_a^r, d_a^z)$ is the solution of (38) without the square bracketed term in the right-hand side, while the full step $(d_a^v, d_a^r, d_a^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^\top r$. The change of variables is done by setting

$$ \hat{z} := Q_p z \quad \text{and} \quad \hat{y} := Q_{p-1} J^\top r. $$

Recall that $\hat{z} = \hat{y}$, which allows to symmetrize the linear system under the form

$$ \begin{pmatrix} M & -H^\top & 0 \\ -H & 0 & J \\ 0 & J^\top & Q_{p2} \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ d^z \end{pmatrix} = \begin{pmatrix} -Mv - f + H^\top r \\ Hv + w - Jz \\ -J^\top r - [Q_p(\hat{z}^{-1} \circ (\hat{d}_a^z \circ \hat{d}_a^y)) - 2\sigma\mu z^{-1}] \end{pmatrix}. \quad (39) $$
• $r_{i0} > \max\{\|\tilde{r}_i\|, \|\tilde{r}_i\|\}$ and $u_i = 0$,

• $u_{i0} > \|\tilde{u}_i\| + \|\tilde{u}_i\|$ and $r_i = 0$,

• $r_{i0} = \max\{\|\tilde{r}_i\|, \|\tilde{r}_i\|\} > 0$ and $u_{i0} = \|\tilde{u}_i\| + \|\tilde{u}_i\| > 0$.

As in Section 4, under the strict complementarity assumption, the index set $\{1, \dots, n\}$ can be partitioned into two partitions $(\mathbb{B}, \mathbb{N}, \mathbb{R})$ and $(\hat{\mathbb{B}}, \hat{\mathbb{N}}, \hat{\mathbb{R}})$, whose definition is a direct extension of the previous one to the current framework. Let $\hat{Z}$ and $\hat{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:

$$ \bar{\mathbb{B}} = \{i : \exists r \in \hat{R}, (r_{i,0}, \bar{r}_i) \in \text{int}(\mathcal{L})\}, \quad \bar{\mathbb{N}} = \{i : \exists z \in \hat{Z}, (\bar{t}_i, \bar{u}_i) \in \text{int}(\mathcal{L})\}, $$

$$ \bar{\mathbb{R}} = \{i : \exists (z, r) \in \hat{Z} \times \hat{R}, ((\bar{t}_i, \bar{u}_i), (r_{i,0}, \bar{r}_i)) \in \text{bd}(\mathcal{L}^2) \setminus \{0\}\}. $$

We can then define the analytic center of the optimal set $\hat{R}$ as follows. If $\hat{R}$ is reduced to a singleton, then it is this point, otherwise it is the minimum of the problem

$$ \min_{r \in \text{ri}(\hat{R})} - \sum_{i \in \bar{\mathbb{B}}} \log \det(r_{i,0}, \bar{r}_i) - \sum_{i \in \bar{\mathbb{B}}} \log \det(r_{i,0}, \tilde{r}_i) - \sum_{i \in \hat{R}, i \in \bar{\mathbb{R}}} \log r_{i,0}. \quad (37) $$

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.

<caption>**Algorithm 2** One iteration of Mehrotra primal-dual algorithm for solving a RF problem</caption>

**Parameters:** $\eta_1 \in (0, 1)$, $\eta_2 \ge 1$, $\eta_3 \ge 1$, $\tau_1 \in (0, 1)$, $\tau_2 \in (0, 1 - \tau_1)$, $\text{tol} > 0$
1: **if** the stopping criterion (32) with $u = Jz$ is satisfied **then** return $(v, z, r)$ as primal-dual solution of (33);
2: Set $\mu \leftarrow r^\top J z / n$;
3: Compute $(d_a^v, d_a^r, d_a^z)$ solution of (38) (resp. (39)-(40)) without the square bracketed term;
4: Find the greatest $\alpha_a \in (0, 1]$ such that $(z, J^\top r) + \alpha_a (d_a^z, J^\top d_a^r) \in \mathcal{L}^{4n}$;
5: Set $\mu_a \leftarrow (r + \alpha_a d_a^r)^\top J (z + \alpha_a d_a^z) / n$;
6: **if** $\mu > \eta_1$ **then** set $e \leftarrow \max\{1, \eta_2 \alpha_a^2\}$ **else** set $e \leftarrow \eta_3$;
7: Set $\sigma \leftarrow \min\{1, (\mu_a / \mu)^e\}$;
8: Compute $(d_a^v, d_a^r, d_a^z)$ solution of (38) (resp. (39)-(40));
9: Set $\tau \leftarrow \tau_1 + \alpha_a \tau_2$;
10: Find the greatest $\alpha \in (0, 1]$ such that $(z, J^\top r) + \alpha (d_a^z, J^\top d_a^r) \in (1 - \tau)(u, r) + \mathcal{L}^{4n}$;
11: Set $(v, z, r) \leftarrow (v, z, r) + \alpha(d_a^v, d_a^z, d_a^r)$ and goto 1.

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

$$ \begin{pmatrix} M & -H^\top & 0 \\ -H & 0 & J \\ 0 & ZJ^\top & R \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ d^z \end{pmatrix} = \begin{pmatrix} -Mv - f + H^\top r \\ Hv + w - Jz \\ -z \circ J^\top r - [d_a^z \circ (J^\top d_a^r) - 2\sigma\mu e] \end{pmatrix}, \quad (38) $$

where $Z = \text{Arw}(z)$, $R = \text{Arw}(J^\top r)$ and $\mu = \frac{r^\top J z}{n}$. As in Algorithm 1, the affine scaling direction $(d_a^v, d_a^r, d_a^z)$ is the solution of (38) without the square bracketed term in the right-hand side, while the full step $(d_a^v, d_a^r, d_a^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^\top r$. The change of variables is done by setting

$$ \hat{z} := Q_p z \quad \text{and} \quad \hat{y} := Q_{p-1} J^\top r. $$

Recall that $\hat{z} = \hat{y}$, which allows to symmetrize the linear system under the form

$$ \begin{pmatrix} M & -H^\top & 0 \\ -H & 0 & J \\ 0 & J^\top & Q_{p2} \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ d^z \end{pmatrix} = \begin{pmatrix} -Mv - f + H^\top r \\ Hv + w - Jz \\ -J^\top r - [Q_p(\hat{z}^{-1} \circ (\hat{d}_a^z \circ \hat{d}_a^y)) - 2\sigma\mu z^{-1}] \end{pmatrix}. \quad (39) $$
<table><caption>Table 4: Sizes of rolling friction problems of the FCLIB collection</caption><thead><tr><th>Family</th><th># problems</th><th>n</th><th>m</th></tr></thead><tbody><tr><td>Chute</td><td>155</td><td>[4, 1372]</td><td>[768, 6528]</td></tr><tr><td>PrimitiveSoup</td><td>171</td><td>[37, 2269]</td><td>6000</td></tr><tr><td>SpherePile</td><td>200</td><td>[2, 542]</td><td>[24, 1500]</td></tr></tbody></table>

Because of the ill-conditioning of the matrix $Q_{p^2}$, an equivalent form of (39) has been considered:

$$ \begin{pmatrix} M & -H^{\top} & 0 \\ -H & 0 & JQ_{p-1} \\ 0 & Q_{p-1}J^{\top} & I \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ \hat{d}^z \end{pmatrix} = \begin{pmatrix} -Mv - f + H^{\top}r \\ Hv + w - Jz \\ -\tilde{y} - [\hat{z}^{-1} \circ (\hat{d}_a^z \circ \tilde{d}_a^y - 2\sigma\mu e)] \end{pmatrix}, \quad (40) $$

where $\hat{d}^z = 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).

<figure><img src="image_4.png" /><figcaption>Figure 3: Condition number ($\kappa$) 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.</figcaption></figure>

We also tried several reductions to a $2 \times 2$ form as in (28) or (30), leading to matrices of the form

$$ \begin{pmatrix} M & -H^{\top} \\ -H & -JQ_{p-2}J^{\top} \end{pmatrix} \quad \text{or} \quad \begin{pmatrix} M & \hat{H}^{\top} \\ -\hat{H} & -I \end{pmatrix}, \quad (41) $$

where $\hat{H} = P^{-1}H$ and $PP^{\top} = JQ_{p-2}J^{\top}$. 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 \times 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.

# 7 Conclusion

At the beginning of this work, the solution of the relaxed convex form of the friction contact models seemed very simple. However, after our initial experience with existing software and our first implementation of Mehrotra's
<table><caption>Table 4: Sizes of rolling friction problems of the FCLIB collection</caption><thead><tr><th>Family</th><th># problems</th><th>n</th><th>m</th></tr></thead><tbody><tr><td>Chute</td><td>155</td><td>[4, 1372]</td><td>[768, 6528]</td></tr><tr><td>PrimitiveSoup</td><td>171</td><td>[37, 2269]</td><td>6000</td></tr><tr><td>SpherePile</td><td>200</td><td>[2, 542]</td><td>[24, 1500]</td></tr></tbody></table>

Because of the ill-conditioning of the matrix $Q_{p^2}$, an equivalent form of (39) has been considered:

$$ \begin{pmatrix} M & -H^\top & 0 \\ -H & 0 & JQ_{p-1} \\ 0 & Q_{p-1}J^\top & I \end{pmatrix} \begin{pmatrix} d^v \\ d^r \\ \hat{d}^z \end{pmatrix} = \begin{pmatrix} -Mv - f + H^\top r \\ Hv + w - Jz \\ -\tilde{y} - [\hat{z}^{-1} \circ (\hat{d}_a^z \circ \tilde{d}_a^y - 2\sigma\mu e)] \end{pmatrix}, \quad (40) $$

where $\hat{d}^z = 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).

<figure><img src="image_4.png" /><figcaption>Figure 3: Condition number ($\kappa$) 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.</figcaption></figure>

We also tried several reductions to a $2 \times 2$ form as in (28) or (30), leading to matrices of the form

$$ \begin{pmatrix} M & -H^\top \\ -H & -JQ_{p-2}J^\top \end{pmatrix} \quad \text{or} \quad \begin{pmatrix} M & \hat{H}^\top \\ -\hat{H} & -I \end{pmatrix}, \quad (41) $$

where $\hat{H} = P^{-1}H$ and $PP^\top = JQ_{p-2}J^\top$. 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 \times 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.

# 7 Conclusion

At the beginning of this work, the solution of the relaxed convex form of the friction contact models seemed very simple. However, after our initial experience with existing software and our first implementation of Mehrotra's
<table><caption>Table 5: Number of successes and failures when solving the RF problems with tol = 10⁻¹⁰</caption><thead><tr><th>linear system</th><th># success</th><th># failure 1</th><th># failure 2</th><th>max res</th></tr></thead><tbody><tr><td>(38)</td><td>391</td><td>135</td><td>0</td><td>3e-3</td></tr><tr><td>(39)</td><td>451</td><td>0</td><td>75</td><td>4e-9</td></tr><tr><td>(40)</td><td>508</td><td>0</td><td>18</td><td>5e-9</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Performance profiles for the three linear systems for the RF model</figcaption></figure>

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 $\mu > 0$:

$$ 
\begin{aligned}
Mv + f &= H^\top r, \\
Hv + w + \phi(s) &= u, \\
s_i &= \|\bar{u}_i\|, \quad i \in \{1, \dots, n\}, \\
u \circ r &= 2\mu e.
\end{aligned}
 $$

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.
<table><caption>Table 5: Number of successes and failures when solving the RF problems with tol = 10⁻¹⁰</caption><thead><tr><th>linear system</th><th># success</th><th># failure 1</th><th># failure 2</th><th>max res</th></tr></thead><tbody><tr><td>(38)</td><td>391</td><td>135</td><td>0</td><td>3e-3</td></tr><tr><td>(39)</td><td>451</td><td>0</td><td>75</td><td>4e-9</td></tr><tr><td>(40)</td><td>508</td><td>0</td><td>18</td><td>5e-9</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Performance profiles for the three linear systems for the RF model</figcaption></figure>

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 $\mu > 0$:

$$ 
\begin{aligned}
Mv + f &= H^\top r, \\
Hv + w + \phi(s) &= u, \\
s_i &= \|\bar{u}_i\|, \quad i \in \{1, \dots, n\}, \\
u \circ r &= 2\mu e.
\end{aligned}
 $$

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.
References

Acary, V. et al. (2011). "A formulation of the linear discrete Coulomb friction problem via convex optimization". In: *ZAMM - Journal of Applied Mathematics and Mechanics* / *Zeitschrift für Angewandte Mathematik und Mechanik* 91.2, pp. 155-175. ISSN: 0044-2267. [DOI], [HAL].

Acary, V. and F. Bourrier (2021). "Coulomb friction with rolling resistance as a cone complementarity problem". In: *European Journal of Mechanics - A/Solids* 85, p. 104046. ISSN: 0997-7538. [DOI], [OA].

Acary, V., M. Brémond, and O. Huber (2018). "On solving contact problems with Coulomb friction: formulations and numerical comparisons". In: *Advanced Topics in Nonsmooth Dynamics*. Springer International Publishing, pp. 375-457. ISBN: 978-3-319-75971-5; 978-3-319-75972-2. [DOI], [HAL].

Acary, V. and B. Brogliato (2008). *Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics*. Springer Science & Business Media.

Acary, V. and F. Cadoux (2013). "Applications of an existence result for the Coulomb friction problem". In: *Lecture Notes in Applied and Computational Mechanics*. Springer Berlin Heidelberg, pp. 45-66. [DOI], [OA].

Acary, V. et al. (2014). *FCLIB: a collection of discrete 3D Frictional Contact problems*. Technical Report RT-0444. INRIA, p. 34.

Alart, P. and A. Curnier (1991). "A mixed formulation for frictional contact problems prone to Newton like solution method". In: *Computer Methods in Applied Mechanics and Engineering* 92.3, pp. 353-375. [DOI], [HAL].

Alizadeh, F. and D. Goldfarb (2003). "Second-order cone programming". In: *Mathematical Programming* 95.1. ISMP 2000, Part 3 (Atlanta, GA), pp. 3-51. ISSN: 0025-5610. [DOI].

Amestoy, P. R. et al. (2001). "A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling". In: *SIAM Journal on Matrix Analysis and Applications* 23.1, pp. 15-41. [DOI], [OA].

Andersen, E., C. Roos, and T. Terlaky (2003). "On implementing a primal-dual interior-point method for conic quadratic optimization". In: *Mathematical Programming* 95.2. Computational semidefinite and second order cone programming: the state of the art, pp. 249-277. ISSN: 0025-5610, 1436-4646. [DOI].

Andersen, E. D. and K. D. Andersen (2000). "The Mosek interior point optimizer for linear programming: an implementation of the homogeneous algorithm". In: *High Performance Optimization*. Vol. 33. Appl. Optim. Springer US, pp. 197-232. ISBN: 0-7923-6013-3. [DOI].

Auslender, A. and M. Teboulle (2003). *Asymptotic cones and functions in optimization and variational inequalities*. Springer Monographs in Mathematics. Springer-Verlag, pp. xii+249. ISBN: 0-387-95520-8. [DOI], [OA].

Benson, H. Y. and R. J. Vanderbei (2003). "Solving problems with semidefinite and related constraints using interior-point methods for nonlinear programming". In: *Mathematical Programming* 95.2. Computational semidefinite and second order cone programming: the state of the art, pp. 279-302. ISSN: 0025-5610. [DOI].

Bicchi, A. and V. Kumar (2000). "Robotic grasping and contact: A review". In: *Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No.00CH37065)*. Vol. 1. IEEE. IEEE, pp. 348-353. [DOI].

Brogliato, B. (1999). *Nonsmooth mechanics*. Vol. 3. Springer.

Buss, M., L. Faybusovich, and J. Moore (1997). "Recursive algorithms for real-time grasping force optimization". In: *Proceedings of International Conference on Robotics and Automation*. Vol. 1. IEEE. IEEE, pp. 682-687. [DOI].

Buss, M., H. Hashimoto, and J. Moore (1996). "Dextrous hand grasping force optimization". In: *IEEE Transactions on Robotics and Automation* 12.3, pp. 406-418. [DOI].

Byrd, R. H., J. Nocedal, and R. A. Waltz (2006). "KNITRO: An integrated package for nonlinear optimization". In: *Large-Scale Nonlinear Optimization*. Vol. 83. Nonconvex Optim. Appl. Springer US, pp. 35-59. [DOI].

Cai, Z. and K.-C. Toh (2006). "Solving second order cone programming via a reduced augmented system approach". In: *SIAM Journal on Optimization* 17.3, pp. 711-737. ISSN: 1052-6234. [DOI].

Cambou, B., M. Jean, and F. Radjaï (2013). *Micromechanics of granular materials*. Wiley. [DOI], [OA].

De Klerk, E., C. Roos, and T. Terlaky (1998). "Infeasible-start semidefinite programming algorithms via self-dual embeddings". In: *Topics in semidefinite and interior-point methods (Toronto, ON, 1996)*. Vol. 18. Fields Inst. Commun. American Mathematical Society, pp. 215-236. [DOI].

Dolan, E. D. and J. J. Moré (2002). "Benchmarking optimization software with performance profiles". In: *Mathematical Programming* 91.2, pp. 201-213. ISSN: 0025-5610. [DOI], [ARXIV].
References

Acary, V. et al. (2011). "A formulation of the linear discrete Coulomb friction problem via convex optimization". In: *ZAMM - Journal of Applied Mathematics and Mechanics* / *Zeitschrift für Angewandte Mathematik und Mechanik* 91.2, pp. 155-175. ISSN: 0044-2267. [DOI], [HAL].

Acary, V. and F. Bourrier (2021). "Coulomb friction with rolling resistance as a cone complementarity problem". In: *European Journal of Mechanics - A/Solids* 85, p. 104046. ISSN: 0997-7538. [DOI], [OA].

Acary, V., M. Brémond, and O. Huber (2018). "On solving contact problems with Coulomb friction: formulations and numerical comparisons". In: *Advanced Topics in Nonsmooth Dynamics*. Springer International Publishing, pp. 375-457. ISBN: 978-3-319-75971-5; 978-3-319-75972-2. [DOI], [HAL].

Acary, V. and B. Brogliato (2008). *Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics*. Springer Science & Business Media.

Acary, V. and F. Cadoux (2013). "Applications of an existence result for the Coulomb friction problem". In: *Lecture Notes in Applied and Computational Mechanics*. Springer Berlin Heidelberg, pp. 45-66. [DOI], [OA].

Acary, V. et al. (2014). *FCLIB: a collection of discrete 3D Frictional Contact problems*. Technical Report RT-0444. INRIA, p. 34.

Alart, P. and A. Curnier (1991). "A mixed formulation for frictional contact problems prone to Newton like solution method". In: *Computer Methods in Applied Mechanics and Engineering* 92.3, pp. 353-375. [DOI], [HAL].

Alizadeh, F. and D. Goldfarb (2003). "Second-order cone programming". In: *Mathematical Programming* 95.1. ISMP 2000, Part 3 (Atlanta, GA), pp. 3-51. ISSN: 0025-5610. [DOI].

Amestoy, P. R. et al. (2001). "A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling". In: *SIAM Journal on Matrix Analysis and Applications* 23.1, pp. 15-41. [DOI], [OA].

Andersen, E., C. Roos, and T. Terlaky (2003). "On implementing a primal-dual interior-point method for conic quadratic optimization". In: *Mathematical Programming* 95.2. Computational semidefinite and second order cone programming: the state of the art, pp. 249-277. ISSN: 0025-5610, 1436-4646. [DOI].

Andersen, E. D. and K. D. Andersen (2000). "The Mosek interior point optimizer for linear programming: an implementation of the homogeneous algorithm". In: *High Performance Optimization*. Vol. 33. Appl. Optim. Springer US, pp. 197-232. ISBN: 0-7923-6013-3. [DOI].

Auslender, A. and M. Teboulle (2003). *Asymptotic cones and functions in optimization and variational inequalities*. Springer Monographs in Mathematics. Springer-Verlag, pp. xii+249. ISBN: 0-387-95520-8. [DOI], [OA].

Benson, H. Y. and R. J. Vanderbei (2003). "Solving problems with semidefinite and related constraints using interior-point methods for nonlinear programming". In: *Mathematical Programming* 95.2. Computational semidefinite and second order cone programming: the state of the art, pp. 279-302. ISSN: 0025-5610. [DOI].

Bicchi, A. and V. Kumar (2000). "Robotic grasping and contact: A review". In: *Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No.00CH37065)*. Vol. 1. IEEE. IEEE, pp. 348-353. [DOI].

Brogliato, B. (1999). *Nonsmooth mechanics*. Vol. 3. Springer.

Buss, M., L. Faybusovich, and J. Moore (1997). "Recursive algorithms for real-time grasping force optimization". In: *Proceedings of International Conference on Robotics and Automation*. Vol. 1. IEEE. IEEE, pp. 682-687. [DOI].

Buss, M., H. Hashimoto, and J. Moore (1996). "Dextrous hand grasping force optimization". In: *IEEE Transactions on Robotics and Automation* 12.3, pp. 406-418. [DOI].

Byrd, R. H., J. Nocedal, and R. A. Waltz (2006). "KNITRO: An integrated package for nonlinear optimization". In: *Large-Scale Nonlinear Optimization*. Vol. 83. Nonconvex Optim. Appl. Springer US, pp. 35-59. [DOI].

Cai, Z. and K.-C. Toh (2006). "Solving second order cone programming via a reduced augmented system approach". In: *SIAM Journal on Optimization* 17.3, pp. 711-737. ISSN: 1052-6234. [DOI].

Cambou, B., M. Jean, and F. Radjaï (2013). *Micromechanics of granular materials*. Wiley. [DOI], [OA].

De Klerk, E., C. Roos, and T. Terlaky (1998). "Infeasible-start semidefinite programming algorithms via self-dual embeddings". In: *Topics in semidefinite and interior-point methods (Toronto, ON, 1996)*. Vol. 18. Fields Inst. Commun. American Mathematical Society, pp. 215-236. [DOI].

Dolan, E. D. and J. J. Moré (2002). "Benchmarking optimization software with performance profiles". In: *Mathematical Programming* 91.2, pp. 201-213. ISSN: 0025-5610. [DOI], [ARXIV].
Duff, I. S. (2004). "MA57—a code for the solution of sparse symmetric definite and indefinite systems". In: *ACM Transactions on Mathematical Software* 30.2, pp. 118–144. [DOI].

Gay, D. M. (2015). "The AMPL modeling language: an aid to formulating and solving optimization problems". In: *Springer Proceedings in Mathematics & Statistics*. Vol. 134. Springer Proc. Math. Stat. Springer International Publishing, pp. 95–116. [DOI], [OA].

Ghannad, A., D. Orban, and M. A. Saunders (2022). "Linear systems arising in interior methods for convex optimization: a symmetric formulation with bounded condition number". In: *Optimization Methods and Software* 37.4, pp. 1344–1369. [DOI].

Halická, M., E. de Klerk, and C. Roos (2002). "On the convergence of the central path in semidefinite optimization". In: *SIAM Journal on Optimization* 12.4, pp. 1090–1099. ISSN: 1052-6234. [DOI], [OA].

— (2005). "Limiting behavior of the central path in semidefinite optimization". In: *Optimization Methods and Software* 20.1, pp. 99–113. ISSN: 1055-6788. [DOI], [OA].

Han, L., J. Trinkle, and Z. Li (2000). "Grasp analysis as linear matrix inequality problems". In: *IEEE Transactions on Robotics and Automation* 16.6, pp. 663–674. [DOI].

Johnson, K. L. and K. L. Johnson (1987). *Contact mechanics*. Cambridge university press.

Kleinert, J., B. Simeon, and M. Obermayr (2014). "An inexact interior point method for the large-scale simulation of granular material". In: *Computer Methods in Applied Mechanics and Engineering* 278, pp. 567–598. ISSN: 0045-7825. [DOI], [OA].

Klerk, E. de, C. Roos, and T. Terlaky (1997). "Initialization in semidefinite programming via a self-dual skew-symmetric embedding". In: *Operations Research Letters* 20.5, pp. 213–221. ISSN: 0167-6377. [DOI], [OA].

Kučera, R. et al. (2013). "An interior-point algorithm for the minimization arising from 3D contact problems with friction". In: *Optimization Methods and Software* 28.6, pp. 1195–1217. [DOI].

Mangasarian, O. L. (1969). *Nonlinear programming*. McGraw-Hill Book Co., New York-London-Sydney, pp. xiii+220.

Mangoni, D., A. Tasora, and R. Garziera (2018). "A primal-dual predictor-corrector interior point method for non-smooth contact dynamics". In: *Computer Methods in Applied Mechanics and Engineering* 330, pp. 351–367. ISSN: 0045-7825. [DOI].

Mehrotra, S. (1992). "On the implementation of a primal-dual interior point method". In: *SIAM Journal on Optimization* 2.4, pp. 575–601. ISSN: 1052-6234. [DOI].

Melanz, D. et al. (2017). "A comparison of numerical methods for solving multibody dynamics problems with frictional contact modeled via differential variational inequalities". In: *Computer Methods in Applied Mechanics and Engineering* 320, pp. 668–693. ISSN: 0045-7825. [DOI], [OA].

Milnor, J. (1968). *Singular points of complex hypersurfaces*. Annals of Mathematics Studies, No. 61. Princeton University Press, pp. iii+122. [DOI].

Mohammad-Nezhad, A. and T. Terlaky (2021). "On the sensitivity of the optimal partition for parametric second-order conic optimization". In: *Mathematical Programming* 189.1–2, pp. 491–525. ISSN: 0025-5610. [DOI], [ARXIV].

Monteiro, R. D. and F. Zou (1998). "On the existence and convergence of the central path for convex programming and some duality results". In: *Computational Optimization and Applications* 10.1, pp. 51–77. ISSN: 0926-6003. [DOI].

Nesterov, Y. E. and M. J. Todd (1997). "Self-scaled barriers and interior-point methods for convex programming". In: *Mathematics of Operations Research* 22.1, pp. 1–42. ISSN: 0364-765X. [DOI], [OA].

Peng, J., C. Roos, and T. Terlaky (2002). *Self-regularity: a new paradigm for primal-dual interior-point algorithms*. Princeton Series in Applied Mathematics. Princeton University Press, pp. xviii+185. ISBN: 0-691-09192-7; 0-691-09193-5. [DOI].

Pfeiffer, F. and C. Glocker (1996). *Multibody dynamics with unilateral contacts*. Wiley. [DOI].

Popov, V. L. et al. (2010). *Contact mechanics and friction*. Springer.

Radjai, F. and F. Dubois (2011). *Discrete-element modeling of granular materials*. Wiley-Iste.

Rockafellar, R. T. (1970). *Convex analysis*. Princeton Mathematical Series, No. 28. Princeton University Press, pp. xviii+451. [DOI].

Skarpness, B. and V. A. Sposito (1982). "A note on Gordan's theorem over cone domains". In: *International Journal of Mathematics and Mathematical Sciences* 5.4, pp. 809–812. ISSN: 0161-1712. [DOI], [OA].

Tütüncü, R. H., K. C. Toh, and M. J. Todd (2003). "Solving semidefinite-quadratic-linear programs using SDPT3". In: *Mathematical Programming* 95.2. Computational semidefinite and second order cone programming: the state of the art, pp. 189–217. ISSN: 0025-5610. [DOI], [OA].
Duff, I. S. (2004). "MA57—a code for the solution of sparse symmetric definite and indefinite systems". In: *ACM Transactions on Mathematical Software* 30.2, pp. 118–144. [DOI].

Gay, D. M. (2015). "The AMPL modeling language: an aid to formulating and solving optimization problems". In: *Springer Proceedings in Mathematics & Statistics*. Vol. 134. Springer Proc. Math. Stat. Springer International Publishing, pp. 95–116. [DOI], [OA].

Ghannad, A., D. Orban, and M. A. Saunders (2022). "Linear systems arising in interior methods for convex optimization: a symmetric formulation with bounded condition number". In: *Optimization Methods and Software* 37.4, pp. 1344–1369. [DOI].

Halická, M., E. de Klerk, and C. Roos (2002). "On the convergence of the central path in semidefinite optimization". In: *SIAM Journal on Optimization* 12.4, pp. 1090–1099. ISSN: 1052-6234. [DOI], [OA].

— (2005). "Limiting behavior of the central path in semidefinite optimization". In: *Optimization Methods and Software* 20.1, pp. 99–113. ISSN: 1055-6788. [DOI], [OA].

Han, L., J. Trinkle, and Z. Li (2000). "Grasp analysis as linear matrix inequality problems". In: *IEEE Transactions on Robotics and Automation* 16.6, pp. 663–674. [DOI].

Johnson, K. L. and K. L. Johnson (1987). *Contact mechanics*. Cambridge university press.

Kleinert, J., B. Simeon, and M. Obermayr (2014). "An inexact interior point method for the large-scale simulation of granular material". In: *Computer Methods in Applied Mechanics and Engineering* 278, pp. 567–598. ISSN: 0045-7825. [DOI], [OA].

Klerk, E. de, C. Roos, and T. Terlaky (1997). "Initialization in semidefinite programming via a self-dual skew-symmetric embedding". In: *Operations Research Letters* 20.5, pp. 213–221. ISSN: 0167-6377. [DOI], [OA].

Kučera, R. et al. (2013). "An interior-point algorithm for the minimization arising from 3D contact problems with friction". In: *Optimization Methods and Software* 28.6, pp. 1195–1217. [DOI].

Mangasarian, O. L. (1969). *Nonlinear programming*. McGraw-Hill Book Co., New York-London-Sydney, pp. xiii+220.

Mangoni, D., A. Tasora, and R. Garziera (2018). "A primal-dual predictor-corrector interior point method for non-smooth contact dynamics". In: *Computer Methods in Applied Mechanics and Engineering* 330, pp. 351–367. ISSN: 0045-7825. [DOI].

Mehrotra, S. (1992). "On the implementation of a primal-dual interior point method". In: *SIAM Journal on Optimization* 2.4, pp. 575–601. ISSN: 1052-6234. [DOI].

Melanz, D. et al. (2017). "A comparison of numerical methods for solving multibody dynamics problems with frictional contact modeled via differential variational inequalities". In: *Computer Methods in Applied Mechanics and Engineering* 320, pp. 668–693. ISSN: 0045-7825. [DOI], [OA].

Milnor, J. (1968). *Singular points of complex hypersurfaces*. Annals of Mathematics Studies, No. 61. Princeton University Press, pp. iii+122. [DOI].

Mohammad-Nezhad, A. and T. Terlaky (2021). "On the sensitivity of the optimal partition for parametric second-order conic optimization". In: *Mathematical Programming* 189.1–2, pp. 491–525. ISSN: 0025-5610. [DOI], [ARXIV].

Monteiro, R. D. and F. Zou (1998). "On the existence and convergence of the central path for convex programming and some duality results". In: *Computational Optimization and Applications* 10.1, pp. 51–77. ISSN: 0926-6003. [DOI].

Nesterov, Y. E. and M. J. Todd (1997). "Self-scaled barriers and interior-point methods for convex programming". In: *Mathematics of Operations Research* 22.1, pp. 1–42. ISSN: 0364-765X. [DOI], [OA].

Peng, J., C. Roos, and T. Terlaky (2002). *Self-regularity: a new paradigm for primal-dual interior-point algorithms*. Princeton Series in Applied Mathematics. Princeton University Press, pp. xviii+185. ISBN: 0-691-09192-7; 0-691-09193-5. [DOI].

Pfeiffer, F. and C. Glocker (1996). *Multibody dynamics with unilateral contacts*. Wiley. [DOI].

Popov, V. L. et al. (2010). *Contact mechanics and friction*. Springer.

Radjai, F. and F. Dubois (2011). *Discrete-element modeling of granular materials*. Wiley-Iste.

Rockafellar, R. T. (1970). *Convex analysis*. Princeton Mathematical Series, No. 28. Princeton University Press, pp. xviii+451. [DOI].

Skarpness, B. and V. A. Sposito (1982). "A note on Gordan's theorem over cone domains". In: *International Journal of Mathematics and Mathematical Sciences* 5.4, pp. 809–812. ISSN: 0161-1712. [DOI], [OA].

Tütüncü, R. H., K. C. Toh, and M. J. Todd (2003). "Solving semidefinite-quadratic-linear programs using SDPT3". In: *Mathematical Programming* 95.2. Computational semidefinite and second order cone programming: the state of the art, pp. 189–217. ISSN: 0025-5610. [DOI], [OA].
# A Euclidean Jordan algebra

Let us consider the set $\mathcal{K} = \prod_{i=1}^{n} \mathcal{K}_i$ where $\mathcal{K}_i$ is an $n_i$-dimensional Lorentz cone defined by

$$ \mathcal{K}_i = \{ x = (x_{i0}, \bar{x}_i^\top) \in \mathbb{R} \times \mathbb{R}^{n_i-1} : \|\bar{x}_i\| \le x_{i0} \}. $$

Let $N = \sum_{i=1}^{n} n_i$. For $x \in \mathbb{R}^N$, we denote $x = (x_1, \dots, x_n)$, where $x_i = (x_{i0}, \bar{x}_i)$. For two vectors $x$ and $y$ in $\mathbb{R}^N$, the Jordan product is defined by

$$ x \circ y = \begin{pmatrix} x_1 \circ y_1 \\ \vdots \\ x_n \circ y_n \end{pmatrix}, \quad \text{where } x_i \circ y_i = \begin{pmatrix} x_i^\top y_i \\ x_{i0} \bar{y}_i + y_{i0} \bar{x}_i \end{pmatrix}, \quad \text{for } i = 1, \dots, n. $$

Let $x \in \mathbb{R}^N$ and $x^2 = x \circ x$. A fundamental property of the Jordan algebra for interior-point algorithms, is that the Lorentz cone is the cone of squares, that is $\mathcal{K} = \{x^2 : x \in \mathbb{R}^N\}$, see (Alizadeh and Goldfarb, 2003, pp. 18–19).

For matrices $A$ and $B$, we define

$$ A \oplus B := \begin{pmatrix} A & 0 \\ 0 & B \end{pmatrix}. $$

For $x \in \mathbb{R}^N$, the arrow-shaped matrix is defined by

$$ \mathrm{Arw}(x) = \mathrm{Arw}(x_1) \oplus \dots \oplus \mathrm{Arw}(x_n), \quad \text{where } \mathrm{Arw}(x_i) = \begin{pmatrix} x_{i0} & \bar{x}_i^\top \\ \bar{x}_i & x_{i0}I \end{pmatrix} \quad \text{for } i = 1, \dots, n. $$

For $i \in \{1, \dots, n\}$, let $e_i = (1, 0) \in \mathbb{R} \times \mathbb{R}^{n_i-1}$ and let $e = (e_1, \dots, e_n) \in \mathbb{R}^N$ be the neutral element of the Jordan product. For vectors $x$ and $y$ in $\mathbb{R}^N$, one has

$$ x \circ y = \mathrm{Arw}(x)y = \mathrm{Arw}(x)\mathrm{Arw}(y)e $$

For $i \in \{1, \dots, n\}$, let $\det(x_i) = x_{i0}^2 - \|\bar{x}_i\|^2$ be the determinant of $x_i \in \mathbb{R}^{n_i}$. If $x_i$ is nonsingular, i.e., $\det(x_i) \neq 0$, the inverse of $x_i$ is the unique vector of $\mathbb{R}^{n_i}$ such that $x_i \circ x_i^{-1} = e_i$ and is given by $x_i^{-1} = \det(x_i)^{-1} R_{n_i} x_i$, where $R_{n_i}$ is the reflexion matrix defined by

$$ R_{n_i} = \begin{pmatrix} 1 & 0^\top \\ 0 & -I \end{pmatrix} \in \mathbb{R}^{n_i \times n_i}. $$

If for all $i \in \{1, \dots, n\}$ $x_i$ is nonsingular, we have

$$ x^{-1} = \begin{pmatrix} x_i^{-1} \\ \vdots \\ x_n^{-1} \end{pmatrix} = \mathrm{Arw}(x)^{-1}e. $$

For $i \in \{1, \dots, n\}$, the spectral decomposition of a vector $x_i \in \mathbb{R}^{n_i}$ is defined by $x_i = \lambda_i c_i + \lambda'_i c'_i$, where

$$ \lambda_i = x_{i0} + \|\bar{x}_i\|, \quad \lambda'_i = x_{i0} - \|\bar{x}_i\|, \quad c_i = \frac{1}{2} \begin{pmatrix} \frac{1}{\bar{x}_i} \\ \|\bar{x}_i\| \end{pmatrix} \quad \text{and} \quad c'_i = \frac{1}{2} \begin{pmatrix} \frac{1}{\bar{x}_i} \\ -\|\bar{x}_i\| \end{pmatrix}. $$

The scalars $\lambda_i$ and $\lambda'_i$ are eigenvalues of $\mathrm{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 \in \mathcal{K}_i$ (resp. $x_i \in \mathrm{int}(\mathcal{K}_i)$) if and only if $\lambda'_i \ge 0$ (resp. $\lambda'_i > 0$). If $x_i$ is nonsingular, then

$$ x_i^{-1} = \lambda_i^{-1} c_i + \lambda'_i^{-1} c'_i. $$
# A Euclidean Jordan algebra

Let us consider the set $\mathcal{K} = \prod_{i=1}^{n} \mathcal{K}_i$ where $\mathcal{K}_i$ is an $n_i$-dimensional Lorentz cone defined by

$$ \mathcal{K}_i = \{ x = (x_{i0}, \bar{x}_i^\top) \in \mathbb{R} \times \mathbb{R}^{n_i-1} : \|\bar{x}_i\| \le x_{i0} \}. $$

Let $N = \sum_{i=1}^{n} n_i$. For $x \in \mathbb{R}^N$, we denote $x = (x_1, \dots, x_n)$, where $x_i = (x_{i0}, \bar{x}_i)$. For two vectors $x$ and $y$ in $\mathbb{R}^N$, the Jordan product is defined by

$$ x \circ y = \begin{pmatrix} x_1 \circ y_1 \\ \vdots \\ x_n \circ y_n \end{pmatrix}, \quad \text{where } x_i \circ y_i = \begin{pmatrix} x_i^\top y_i \\ x_{i0} \bar{y}_i + y_{i0} \bar{x}_i \end{pmatrix}, \quad \text{for } i = 1, \dots, n. $$

Let $x \in \mathbb{R}^N$ and $x^2 = x \circ x$. A fundamental property of the Jordan algebra for interior-point algorithms, is that the Lorentz cone is the cone of squares, that is $\mathcal{K} = \{x^2 : x \in \mathbb{R}^N\}$, see (Alizadeh and Goldfarb, 2003, pp. 18–19).

For matrices $A$ and $B$, we define

$$ A \oplus B := \begin{pmatrix} A & 0 \\ 0 & B \end{pmatrix}. $$

For $x \in \mathbb{R}^N$, the arrow-shaped matrix is defined by

$$ \mathrm{Arw}(x) = \mathrm{Arw}(x_1) \oplus \dots \oplus \mathrm{Arw}(x_n), \quad \text{where } \mathrm{Arw}(x_i) = \begin{pmatrix} x_{i0} & \bar{x}_i^\top \\ \bar{x}_i & x_{i0}I \end{pmatrix} \quad \text{for } i = 1, \dots, n. $$

For $i \in \{1, \dots, n\}$, let $e_i = (1, 0) \in \mathbb{R} \times \mathbb{R}^{n_i-1}$ and let $e = (e_1, \dots, e_n) \in \mathbb{R}^N$ be the neutral element of the Jordan product. For vectors $x$ and $y$ in $\mathbb{R}^N$, one has

$$ x \circ y = \mathrm{Arw}(x)y = \mathrm{Arw}(x)\mathrm{Arw}(y)e $$

For $i \in \{1, \dots, n\}$, let $\det(x_i) = x_{i0}^2 - \|\bar{x}_i\|^2$ be the determinant of $x_i \in \mathbb{R}^{n_i}$. If $x_i$ is nonsingular, i.e., $\det(x_i) \neq 0$, the inverse of $x_i$ is the unique vector of $\mathbb{R}^{n_i}$ such that $x_i \circ x_i^{-1} = e_i$ and is given by $x_i^{-1} = \det(x_i)^{-1}R_{n_i}x_i$, where $R_{n_i}$ is the reflexion matrix defined by

$$ R_{n_i} = \begin{pmatrix} 1 & 0^\top \\ 0 & -I \end{pmatrix} \in \mathbb{R}^{n_i \times n_i}. $$

If for all $i \in \{1, \dots, n\}$ $x_i$ is nonsingular, we have

$$ x^{-1} = \begin{pmatrix} x_i^{-1} \\ \vdots \\ x_n^{-1} \end{pmatrix} = \mathrm{Arw}(x)^{-1}e. $$

For $i \in \{1, \dots, n\}$, the spectral decomposition of a vector $x_i \in \mathbb{R}^{n_i}$ is defined by $x_i = \lambda_i c_i + \lambda'_i c'_i$, where

$$ \lambda_i = x_{i0} + \|\bar{x}_i\|, \quad \lambda'_i = x_{i0} - \|\bar{x}_i\|, \quad c_i = \frac{1}{2} \begin{pmatrix} \frac{1}{\bar{x}_i} \\ \|\bar{x}_i\| \end{pmatrix} \quad \text{and} \quad c'_i = \frac{1}{2} \begin{pmatrix} \frac{1}{\bar{x}_i} \\ -\|\bar{x}_i\| \end{pmatrix}. $$

The scalars $\lambda_i$ and $\lambda'_i$ are eigenvalues of $\mathrm{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 \in \mathcal{K}_i$ (resp. $x_i \in \mathrm{int}(\mathcal{K}_i)$) if and only if $\lambda'_i \ge 0$ (resp. $\lambda'_i > 0$). If $x_i$ is nonsingular, then

$$ x_i^{-1} = \lambda_i^{-1} c_i + \lambda'_i^{-1} c'_i. $$
If $x_i \in \text{int}(\mathcal{K}_i)$, then

$$x_i^{1/2} = \lambda_i^{1/2} c_i + \lambda_i^{1/2} c_i'$$

More generally, for a continuous function $f$ we can define

$$f(x_i) = f(\lambda_i)c_i + f(\lambda_i')c_i'$$

At last, the quadratic representation of a vector $x_i \in \mathcal{K}_i$ is defined as

$$Q_{x_i} = \text{Arw}^2(x_i) - \text{Arw}(x_i^2) = 2x_i x_i^\top - \det(x_i) R_{n_i}.$$

The scalars $\lambda_i^2$ and $\lambda_i^{1/2}$ are eigenvalues of $Q_{x_i}$ with the eigenvectors $c_i$ and $c_i'$. In particular, two vectors commute if and only if their quadratic representation matrices commute. For $x \in \mathcal{K}$, we define $Q_x = Q_{x_1} \oplus \dots \oplus Q_{x_n}$. For $i \in \{1, \dots, n\}$ and $x \in \text{int}(\mathcal{K}_i)$, we have $\nabla_x(\log \det(x_i)) = 2x_i^{-1}$ and $\nabla_x^2(\log \det(x_i)) = -2Q_{x_i^{-1}}$.

These operators $\text{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 \subset \mathbb{R}^n$ be a closed pointed convex cone, $A \in \mathbb{R}^{m \times n}$ and $B \in \mathbb{R}^{p \times n}$. One and only one of the following statements is true.

(i) There exists a non-zero $x \in K$ such that $Ax = 0$ and $Bx \le 0$.

(ii) There exists $(y, z) \in \mathbb{R}^m \times \mathbb{R}_+^p$ such that $A^\top y + B^\top z \in \text{int}(K^*)$.

The following lemma provides equivalent formulations of Assumption 3.1.

**Lemma B.2.** Let $\mathcal{F}$ be the product of second order cones of the form (5) or (7). Let $M \in \mathbb{R}^{m \times m}$ be symmetric and positive-definite, $f \in \mathbb{R}^m$, $H \in \mathbb{R}^{nd \times m}$, $w \in \mathbb{R}^n$ and $W = HM^{-1}H^\top$. The following four assertions are equivalent.

(i) There exists $v \in \mathbb{R}^m$ such that $Hv + w \in \text{int}(\mathcal{F}^*)$.

(ii) There exists $(v, t) \in \mathbb{R}^m \times \mathbb{R}_+^p$ such that $Hv + tw \in \text{int}(\mathcal{F}^*)$.

(iii) There does not exist a nonzero $d \in \mathcal{F}$, such that $H^\top d = 0$ and $w^\top d \le 0$.

(iv) There does not exist a nonzero $d \in \mathcal{F}$, such that $Wd = 0$ and $q^\top d \le 0$.

*Proof.* The implication (i) ⇒ (ii) is obvious. The equivalence (ii) ⇔ (iii) is a direct consequence of the fact that $\mathcal{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) \in \mathbb{R}^m \times \mathbb{R}_+$ such that $Hv + tw \in \text{int}(\mathcal{F}^*)$. If $t > 0$, set $v' = \frac{1}{t}v$. Since $Hv + tw = t(Hv' + w)$ and $\text{int}(\mathcal{F}^*)$ is a cone, $Hv' + w \in \text{int}(\mathcal{F}^*)$. Suppose now that $t = 0$. There exists $\epsilon > 0$ such that $Hv + B_\epsilon \subset \mathcal{F}^*$, where $B_\epsilon$ is the open ball centered at 0 with radius $\epsilon$. Since $\mathcal{F}^*$ is a cone, for all $t > 0$, $tHv + B_{t\epsilon} \subset \mathcal{F}^*$. Let us choose $t > \frac{1}{\epsilon}\|w\|$. We then have $tHv + w \in \text{int}(\mathcal{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
If $x_i \in \text{int}(\mathcal{K}_i)$, then

$$x_i^{1/2} = \lambda_i^{1/2} c_i + \lambda_i^{1/2} c_i'$$

More generally, for a continuous function $f$ we can define

$$f(x_i) = f(\lambda_i)c_i + f(\lambda_i')c_i'$$

At last, the quadratic representation of a vector $x_i \in \mathcal{K}_i$ is defined as

$$Q_{x_i} = \text{Arw}^2(x_i) - \text{Arw}(x_i^2) = 2x_i x_i^\top - \det(x_i) R_{n_i}.$$

The scalars $\lambda_i^2$ and $\lambda_i^{1/2}$ are eigenvalues of $Q_{x_i}$ with the eigenvectors $c_i$ and $c_i'$. In particular, two vectors commute if and only if their quadratic representation matrices commute. For $x \in \mathcal{K}$, we define $Q_x = Q_{x_1} \oplus \dots \oplus Q_{x_n}$. For $i \in \{1, \dots, n\}$ and $x \in \text{int}(\mathcal{K}_i)$, we have $\nabla_x(\log \det(x_i)) = 2x_i^{-1}$ and $\nabla_x^2(\log \det(x_i)) = -2Q_{x_i^{-1}}$.

These operators $\text{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 \subset \mathbb{R}^n$ be a closed pointed convex cone, $A \in \mathbb{R}^{m \times n}$ and $B \in \mathbb{R}^{p \times n}$. One and only one of the following statements is true.

(i) There exists a non-zero $x \in K$ such that $Ax = 0$ and $Bx \le 0$.

(ii) There exists $(y, z) \in \mathbb{R}^m \times \mathbb{R}_+^p$ such that $A^\top y + B^\top z \in \text{int}(K^*)$.

The following lemma provides equivalent formulations of Assumption 3.1.

**Lemma B.2.** Let $\mathcal{F}$ be the product of second order cones of the form (5) or (7). Let $M \in \mathbb{R}^{m \times m}$ be symmetric and positive-definite, $f \in \mathbb{R}^m$, $H \in \mathbb{R}^{nd \times m}$, $w \in \mathbb{R}^n$ and $W = HM^{-1}H^\top$. The following four assertions are equivalent.

(i) There exists $v \in \mathbb{R}^m$ such that $Hv + w \in \text{int}(\mathcal{F}^*)$.

(ii) There exists $(v, t) \in \mathbb{R}^m \times \mathbb{R}_+^p$ such that $Hv + tw \in \text{int}(\mathcal{F}^*)$.

(iii) There does not exist a nonzero $d \in \mathcal{F}$, such that $H^\top d = 0$ and $w^\top d \le 0$.

(iv) There does not exist a nonzero $d \in \mathcal{F}$, such that $Wd = 0$ and $q^\top d \le 0$.

*Proof.* The implication (i) $\Rightarrow$ (ii) is obvious. The equivalence (ii) $\Leftrightarrow$ (iii) is a direct consequence of the fact that $\mathcal{F}^*$ is a closed pointed convex cone and of Lemma B.1. The equivalence (iii) $\Leftrightarrow$ (iv) follows from the definitions of $W$ and $q$ in (10), and of the positive definiteness of $M$. It remains to prove that (ii) $\Rightarrow$ (i).

Let $(v, t) \in \mathbb{R}^m \times \mathbb{R}_+$ such that $Hv + tw \in \text{int}(\mathcal{F}^*)$. If $t > 0$, set $v' = \frac{1}{t}v$. Since $Hv + tw = t(Hv' + w)$ and $\text{int}(\mathcal{F}^*)$ is a cone, $Hv' + w \in \text{int}(\mathcal{F}^*)$. Suppose now that $t = 0$. There exists $\epsilon > 0$ such that $Hv + B_\epsilon \subset \mathcal{F}^*$, where $B_\epsilon$ is the open ball centered at 0 with radius $\epsilon$. Since $\mathcal{F}^*$ is a cone, for all $t > 0$, $tHv + B_{t\epsilon} \subset \mathcal{F}^*$. Let us choose $t > \frac{1}{\epsilon}\|w\|$. We then have $tHv + w \in \text{int}(\mathcal{F}^*)$, which proves (i). $\square$

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
$g(r) = \frac{1}{2}r^\top W r + q^\top r + \delta_F(r)$, where $\delta_F$ is the indicator function of the set $\mathcal{F}$. The asymptotic function of $g$ is defined by

$$g_{\infty}(s) = \sup_{t>0} \frac{g(ts) - g(0)}{t} = \begin{cases} q^{\top}s & \text{if } s \in \mathcal{F} \text{ and } Ws = 0, \\ \infty & \text{otherwise,} \end{cases}$$

for $s \in \mathbb{R}^{nd}$, see (Auslender and Teboulle, 2003, Proposition 2.5.2). The optimal set $\hat{R}$ coincides with the set of minima of $g$. It is non-empty and compact if and only if $g$ is coercive, that is $g_{\infty}(s) > 0$ for all nonzero $s \in \mathbb{R}^{nd}$ (Auslender and Teboulle, 2003, Proposition 3.1.3). It follows that $\hat{R}$ is nonempty and compact if and only if there is no nonzero $s \in \mathcal{F}$, such that $Ws = 0$ and $q^\top s \le 0$. By Lemma B.2, this is equivalent to the Slater hypothesis. □

## C Proofs of Propositions 4.1 and 4.2

*Proof of Proposition 4.1.* Let $\mu > 0$. Let us show that the barrier function defined in (16) is coercive. As in the proof of Proposition 3.2, we show that for all nonzero $s \in \mathbb{R}^{nd}$, $(\varphi_{\mu})_{\infty}(s) > 0$. Let us write $\varphi_{\mu} = p + \psi_{\mu}$, where $p$ is the quadratic part. We have $(\varphi_{\mu})_{\infty} = p_{\infty} + (\psi_{\mu})_{\infty}$ (Auslender and Teboulle, 2003, Proposition 2.6.1.). We have

$$p_{\infty}(s) = \sup_{t>0} \frac{p(ts) - p(0)}{t} = \begin{cases} q^{\top}s & \text{if } Ws = 0, \\ \infty & \text{otherwise.} \end{cases}$$

To compute the asymptotic derivative of the function $\psi_{\mu}$, let us look at one element of the sum. Let $s \in \mathbb{R}^d$ and let us define $\rho(s) := - \log \det s$. We have

$$\begin{aligned} \rho_{\infty}(s) &= \sup_{t>0} \frac{1}{t}(\rho(e + ts) - \rho(e)) \\ &= \sup_{t>0} \frac{-1}{t} \log(1 + 2te^{\top}s + t^2 \det s) \\ &= \delta_{\mathcal{L}}(s). \end{aligned}$$

We then have $(\psi_{\mu})_{\infty} = \delta_{\mathcal{L}^n}$. It follows that $(\varphi_{\mu})_{\infty} = g_{\infty}$, where $g$ is defined in the proof of Proposition 3.2. Therefore, the same conclusion holds, that is: under Assumption 3.1, $(\varphi_{\mu})_{\infty}(s) > 0$ for all nonzero $s \in \mathcal{L}^n$. This implies that the set of minima of $\varphi_{\mu}$ is nonempty and compact. Since $\varphi_{\mu}$ is strictly convex, this set is reduced to a singleton.

To show that the central path is bounded, let us proceed to a reasoning by contradiction. Suppose that there exists a positive sequence $\{\mu_k\}$, converging to zero and such that $\{(r_k, u_k)\}$ is unbounded, where we denote $r_k := r(\mu_k)$ and $u_k := u(\mu_k)$. Let $t_k := \|(r_k, u_k)\|$. By taking a subsequence, suppose that $\{\frac{1}{t_k}(r_k, u_k)\} \to (r^*, u^*) \ne 0$. By definition of the central path, for all $k \in \mathbb{N}$ we have

$$Wr_k + q = u_k \quad \text{and} \quad r_k \circ u_k = 2\mu_k e. \qquad (42)$$

Dividing the first equation by $t_k$ and the second one by $t_k^2$, then taking the limit $k \to \infty$, we obtain $Wr^* = u^*$ and $r^* \circ u^* = 0$. Multiplying the first equation by $r^*$ on the left, we obtain $r^{*\top}W r^* = r^{*\top}u^* = 0$. Since $W = HM^{-1}H^\top$ and $M$ is positive definite, it follows that $H^\top r^* = 0$. Multiplying the first equation of (42) by $r_k$ we have $r_k^\top W r_k + q^\top r_k = r_k^\top u_k = 2n\mu_k$. Since $W$ is positive semi-definite, we then get $q^\top r_k \le 2n\mu_k$. Dividing by $t_k$ and passing through the limit, we obtain $q^\top r^* \le 0$. Since for all $k \in \mathbb{N}$, $r_k \in \text{int}(\mathcal{L}^n)$, we also have $r^* \in \mathcal{L}^n$. Finally, we have shown that there exists a non-zero vector $r^* \in \mathcal{L}^n$, such that $H^\top r^* = 0$ and $w^\top r^* = q^\top r^* \le 0$. In view of Lemma B.2, this is equivalent to the fact that the Slater hypothesis does not hold, which is in contradiction with the assumption. □

The proof of Proposition 4.2 is done in two steps. The convergence is proved first. It follows the one of (Halická et al., 2002, Theorem A.3). It is based on the use of the curve selection lemma from algebraic geometry, see (Milnor, 1968, Lemma 3.1). In the second part of the proof, it is shown that the limit is maximally complementary. It is inspired by (Klerk et al., 1997).
$g(r) = \frac{1}{2}r^\top W r + q^\top r + \delta_F(r)$, where $\delta_F$ is the indicator function of the set $\mathcal{F}$. The asymptotic function of $g$ is defined by

$$g_{\infty}(s) = \sup_{t>0} \frac{g(ts) - g(0)}{t} = \begin{cases} q^{\top}s & \text{if } s \in \mathcal{F} \text{ and } Ws = 0, \\ \infty & \text{otherwise,} \end{cases}$$

for $s \in \mathbb{R}^{nd}$, see (Auslender and Teboulle, 2003, Proposition 2.5.2). The optimal set $\hat{R}$ coincides with the set of minima of $g$. It is non-empty and compact if and only if $g$ is coercive, that is $g_{\infty}(s) > 0$ for all nonzero $s \in \mathbb{R}^{nd}$ (Auslender and Teboulle, 2003, Proposition 3.1.3). It follows that $\hat{R}$ is nonempty and compact if and only if there is no nonzero $s \in \mathcal{F}$, such that $Ws = 0$ and $q^\top s \le 0$. By Lemma B.2, this is equivalent to the Slater hypothesis. $\square$

## C Proofs of Propositions 4.1 and 4.2

*Proof of Proposition 4.1.* Let $\mu > 0$. Let us show that the barrier function defined in (16) is coercive. As in the proof of Proposition 3.2, we show that for all nonzero $s \in \mathbb{R}^{nd}$, $(\varphi_{\mu})_{\infty}(s) > 0$. Let us write $\varphi_{\mu} = p + \psi_{\mu}$, where $p$ is the quadratic part. We have $(\varphi_{\mu})_{\infty} = p_{\infty} + (\psi_{\mu})_{\infty}$ (Auslender and Teboulle, 2003, Proposition 2.6.1.). We have

$$p_{\infty}(s) = \sup_{t>0} \frac{p(ts) - p(0)}{t} = \begin{cases} q^{\top}s & \text{if } Ws = 0, \\ \infty & \text{otherwise.} \end{cases}$$

To compute the asymptotic derivative of the function $\psi_{\mu}$, let us look at one element of the sum. Let $s \in \mathbb{R}^d$ and let us define $\rho(s) := -\log \det s$. We have

$$\begin{aligned} \rho_{\infty}(s) &= \sup_{t>0} \frac{1}{t}(\rho(e + ts) - \rho(e)) \\ &= \sup_{t>0} \frac{-1}{t} \log(1 + 2te^{\top}s + t^2 \det s) \\ &= \delta_{\mathcal{L}}(s). \end{aligned}$$

We then have $(\psi_{\mu})_{\infty} = \delta_{\mathcal{L}^n}$. It follows that $(\varphi_{\mu})_{\infty} = g_{\infty}$, where $g$ is defined in the proof of Proposition 3.2. Therefore, the same conclusion holds, that is: under Assumption 3.1, $(\varphi_{\mu})_{\infty}(s) > 0$ for all nonzero $s \in \mathcal{L}^n$. This implies that the set of minima of $\varphi_{\mu}$ is nonempty and compact. Since $\varphi_{\mu}$ is strictly convex, this set is reduced to a singleton.

To show that the central path is bounded, let us proceed to a reasoning by contradiction. Suppose that there exists a positive sequence $\{\mu_k\}$, converging to zero and such that $\{(r_k, u_k)\}$ is unbounded, where we denote $r_k := r(\mu_k)$ and $u_k := u(\mu_k)$. Let $t_k := \|(r_k, u_k)\|$. By taking a subsequence, suppose that $\{\frac{1}{t_k}(r_k, u_k)\} \to (r^*, u^*) \ne 0$. By definition of the central path, for all $k \in \mathbb{N}$ we have

$$Wr_k + q = u_k \quad \text{and} \quad r_k \circ u_k = 2\mu_k e. \tag{42}$$

Dividing the first equation by $t_k$ and the second one by $t_k^2$, then taking the limit $k \to \infty$, we obtain $Wr^* = u^*$ and $r^* \circ u^* = 0$. Multiplying the first equation by $r^*$ on the left, we obtain $r^{*\top}W r^* = r^{*\top}u^* = 0$. Since $W = HM^{-1}H^\top$ and $M$ is positive definite, it follows that $H^\top r^* = 0$. Multiplying the first equation of (42) by $r_k$ we have $r_k^\top W r_k + q^\top r_k = r_k^\top u_k = 2n\mu_k$. Since $W$ is positive semi-definite, we then get $q^\top r_k \le 2n\mu_k$. Dividing by $t_k$ and passing through the limit, we obtain $q^\top r^* \le 0$. Since for all $k \in \mathbb{N}$, $r_k \in \text{int}(\mathcal{L}^n)$, we also have $r^* \in \mathcal{L}^n$. Finally, we have shown that there exists a non-zero vector $r^* \in \mathcal{L}^n$, such that $H^\top r^* = 0$ and $w^\top r^* = q^\top r^* \le 0$. In view of Lemma B.2, this is equivalent to the fact that the Slater hypothesis does not hold, which is in contradiction with the assumption. $\square$

The proof of Proposition 4.2 is done in two steps. The convergence is proved first. It follows the one of (Halická et al., 2002, Theorem A.3). It is based on the use of the curve selection lemma from algebraic geometry, see (Milnor, 1968, Lemma 3.1). In the second part of the proof, it is shown that the limit is maximally complementary. It is inspired by (Klerk et al., 1997).
Proof of Proposition 4.2. Let $(\hat{r}, \hat{u})$ be a limit point of the central path. Since $\hat{u}$ is unique, it remains to show that $r(\mu) \to \hat{r}$ as $\mu \to 0$. Let us define the subsets of $\mathbb{R}^{nd+1}$:

$$ \mathcal{V} = \{(r, \mu) : (W(\hat{r} + r) + q) \circ (\hat{r} + r) = 2\mu e\}, $$

$$ \mathcal{U} = \{(r, \mu) : \mu > 0, \det(\hat{r}_i + r_i) > 0 \text{ and } \hat{r}_{i0} + r_{i0} > 0, i = 1, \dots, n\}. $$

The set $\mathcal{V}$ is a real algebraic set and $\mathcal{U}$ is an open set defined by a finite number of polynomial inequalities. Moreover, $0 \in \mathbb{R}^{nd+1}$ is in the closure of $\mathcal{U} \cap \mathcal{V}$. By the curve selection lemma, there exists $\varepsilon > 0$ and a real analytic curve $p: [0, \varepsilon) \to \mathbb{R}^{nd+1}$, with $p(0) = 0$ and for all $t \in (0, \varepsilon)$, $p(t) \in \mathcal{U} \cap \mathcal{V}$. By setting $p(t) = (\rho(t), \nu(t))$, we have $(W(\hat{r} + \rho(t)) + q) \circ (\hat{r} + \rho(t)) = 2\nu(t)e$, for all $t \in (0, \varepsilon)$. Since the central path is uniquely defined by (15), for all $t \in (0, \varepsilon)$ we have

$$ r(\nu(t)) = \hat{r} + \rho(t). \qquad (43) $$

Since $\nu(0) = 0$, $\nu(t) > 0$ for all $t \in (0, \varepsilon)$ and $\nu$ is analytic on $[0, \varepsilon)$, there exists $\varepsilon' > 0$ such that for all $t \in (0, \varepsilon')$, $\nu'(t) > 0$. It follows that $\nu$ is invertible on $(0, \varepsilon')$. Moreover, for all $t \in (0, \nu(\varepsilon'))$, $\nu^{-1}(t) > 0$ and $\lim_{t \to 0} \nu^{-1}(t) = 0$. We then have, for $\mu > 0$ small enough,

$$ r(\mu) = r(\nu(\nu^{-1}(\mu))) = \hat{r} + \rho((\nu^{-1}(\mu))). $$

It follows that $\lim_{\mu \to 0} r(\mu) = \hat{r}$.

Let us show now that $\hat{r}$ is maximally complementary. That is, for all $i \in \mathbb{B}$, $\hat{r}_i \in \det(\mathcal{L})$ and for all $i \in \mathbb{R} \cup \mathbb{T}_1$, $\hat{r}_i \in \text{bd}(\mathcal{L}) \setminus \{0\}$. Let $r \in \hat{\mathbb{R}}$ be a maximally complementary optimal solution. For all $\mu > 0$, by using (14) and (15) and the positive semi-definiteness of $W$, we have

$$ (r(\mu) - r)^{\top}(u(\mu) - \hat{u}) = (r(\mu) - r)^{\top}W(r(\mu) - r) \geq 0. $$

It follows that

$$ \sum_{i=1}^{n} \hat{u}_{i}^{\top} r_{i}(\mu) + \sum_{i=1}^{n} r_{i}^{\top} u_{i}(\mu) \leq 2n\mu. \qquad (44) $$

By dividing by $2\mu$ and by using the partition of the index set, we obtain

$$ \sum_{i \in \mathbb{N} \cup \mathbb{R} \cup \mathbb{T}_2}^{n} \hat{u}_i^\top u_i^{-1}(\mu) + \sum_{i \in \mathbb{B} \cup \mathbb{R} \cup \mathbb{T}_1}^{n} r_i^\top r_i^{-1}(\mu) \leq n. \qquad (45) $$

Since each term in these sums is nonnegative, we have that for all $i \in \mathbb{B} \cup \mathbb{R} \cup \mathbb{T}_1$ and for all $\mu > 0$, $r_i^\top r_i^{-1}(\mu) \leq n$. By using the inequality (20) in the proof of Theorem 4.5, we deduce that for all $i \in \mathbb{B} \cup \mathbb{R} \cup \mathbb{T}_1$ and all $\mu > 0$, we have

$$ 0 < \frac{1}{n} r_{i,0} \leq r_{i,0}(\mu) + \|\bar{r}_i(\mu)\| < 2r_{i,0}(\mu). $$

Taking the limit as $\mu \to 0$, we have $\hat{r}_{i,0} > 0$. It remains to show that for all $i \in \mathbb{B}$, $\det(\hat{r}_i) \neq 0$. For all $i$ and $\mu > 0$, we have $r_i^\top r_i^{-1}(\mu) = \frac{r_i^\top R r_i(\mu)}{\det r_i(\mu)}$. It follows that for all $i \in \mathbb{B}$ and all $\mu > 0$,

$$ 0 \leq r_i^\top R r_i(\mu) \leq n \det(r_i(\mu)). $$

Passing through the limit $\mu \to 0$, we obtain $0 \leq r_i R \hat{r}_i \leq n \det(\hat{r}_i)$. If $\det(\hat{r}_i) = 0$, we then have

$$ r_{i,0} \hat{r}_{i,0} = \bar{r}_i^\top \bar{r}_i \leq \|\bar{r}_i\| \|\bar{r}_i\| < r_{i,0} \|\bar{r}_i\|, $$

from which we deduce that $\hat{r}_i \notin \mathcal{L}$, which is in contradiction with the fact that $\hat{r} \in \mathcal{K}$ and allows to conclude that $\hat{r}$ is maximally complementary. $\square$
Proof of Proposition 4.2. Let $(\hat{r}, \hat{u})$ be a limit point of the central path. Since $\hat{u}$ is unique, it remains to show that $r(\mu) \to \hat{r}$ as $\mu \to 0$. Let us define the subsets of $\mathbb{R}^{nd+1}$:

$$ \mathcal{V} = \{(r, \mu) : (W(\hat{r} + r) + q) \circ (\hat{r} + r) = 2\mu e\}, $$

$$ \mathcal{U} = \{(r, \mu) : \mu > 0, \det(\hat{r}_i + r_i) > 0 \text{ and } \hat{r}_{i0} + r_{i0} > 0, i = 1, \dots, n\}. $$

The set $\mathcal{V}$ is a real algebraic set and $\mathcal{U}$ is an open set defined by a finite number of polynomial inequalities. Moreover, $0 \in \mathbb{R}^{nd+1}$ is in the closure of $\mathcal{U} \cap \mathcal{V}$. By the curve selection lemma, there exists $\varepsilon > 0$ and a real analytic curve $p: [0, \varepsilon) \to \mathbb{R}^{nd+1}$, with $p(0) = 0$ and for all $t \in (0, \varepsilon)$, $p(t) \in \mathcal{U} \cap \mathcal{V}$. By setting $p(t) = (\rho(t), \nu(t))$, we have $(W(\hat{r} + \rho(t)) + q) \circ (\hat{r} + \rho(t)) = 2\nu(t)e$, for all $t \in (0, \varepsilon)$. Since the central path is uniquely defined by (15), for all $t \in (0, \varepsilon)$ we have

$$ r(\nu(t)) = \hat{r} + \rho(t). \qquad (43) $$

Since $\nu(0) = 0$, $\nu(t) > 0$ for all $t \in (0, \varepsilon)$ and $\nu$ is analytic on $[0, \varepsilon)$, there exists $\varepsilon' > 0$ such that for all $t \in (0, \varepsilon')$, $\nu'(t) > 0$. It follows that $\nu$ is invertible on $(0, \varepsilon')$. Moreover, for all $t \in (0, \nu(\varepsilon'))$, $\nu^{-1}(t) > 0$ and $\lim_{t \to 0} \nu^{-1}(t) = 0$. We then have, for $\mu > 0$ small enough,

$$ r(\mu) = r(\nu(\nu^{-1}(\mu))) = \hat{r} + \rho((\nu^{-1}(\mu))). $$

It follows that $\lim_{\mu \to 0} r(\mu) = \hat{r}$.

Let us show now that $\hat{r}$ is maximally complementary. That is, for all $i \in \mathbb{B}$, $\hat{r}_i \in \det(\mathcal{L})$ and for all $i \in \mathbb{R} \cup \mathbb{T}_1$, $\hat{r}_i \in \text{bd}(\mathcal{L}) \setminus \{0\}$. Let $r \in \hat{\mathbb{R}}$ be a maximally complementary optimal solution. For all $\mu > 0$, by using (14) and (15) and the positive semi-definiteness of $W$, we have

$$ (r(\mu) - r)^{\top}(u(\mu) - \hat{u}) = (r(\mu) - r)^{\top}W(r(\mu) - r) \geq 0. $$

It follows that

$$ \sum_{i=1}^{n} \hat{u}_{i}^{\top} r_{i}(\mu) + \sum_{i=1}^{n} r_{i}^{\top} u_{i}(\mu) \leq 2n\mu. \qquad (44) $$

By dividing by $2\mu$ and by using the partition of the index set, we obtain

$$ \sum_{i \in \mathbb{N} \cup \mathbb{R} \cup \mathbb{T}_2}^{n} \hat{u}_i^\top u_i^{-1}(\mu) + \sum_{i \in \mathbb{B} \cup \mathbb{R} \cup \mathbb{T}_1}^{n} r_i^\top r_i^{-1}(\mu) \leq n. \qquad (45) $$

Since each term in these sums is nonnegative, we have that for all $i \in \mathbb{B} \cup \mathbb{R} \cup \mathbb{T}_1$ and for all $\mu > 0$, $r_i^\top r_i^{-1}(\mu) \leq n$. By using the inequality (20) in the proof of Theorem 4.5, we deduce that for all $i \in \mathbb{B} \cup \mathbb{R} \cup \mathbb{T}_1$ and all $\mu > 0$, we have

$$ 0 < \frac{1}{n} r_{i,0} \leq r_{i,0}(\mu) + \|\bar{r}_i(\mu)\| < 2r_{i,0}(\mu). $$

Taking the limit as $\mu \to 0$, we have $\hat{r}_{i,0} > 0$. It remains to show that for all $i \in \mathbb{B}$, $\det(\hat{r}_i) \neq 0$. For all $i$ and $\mu > 0$, we have $r_i^\top r_i^{-1}(\mu) = \frac{r_i^\top R r_i(\mu)}{\det r_i(\mu)}$. It follows that for all $i \in \mathbb{B}$ and all $\mu > 0$,

$$ 0 \leq r_i^\top R r_i(\mu) \leq n \det(r_i(\mu)). $$

Passing through the limit $\mu \to 0$, we obtain $0 \leq r_i R \hat{r}_i \leq n \det(\hat{r}_i)$. If $\det(\hat{r}_i) = 0$, we then have

$$ r_{i,0} \hat{r}_{i,0} = \bar{r}_i^\top \bar{r}_i \leq \|\bar{r}_i\| \|\bar{r}_i\| < r_{i,0} \|\bar{r}_i\|, $$

from which we deduce that $\hat{r}_i \notin \mathcal{L}$, which is in contradiction with the fact that $\hat{r} \in \mathcal{K}$ and allows to conclude that $\hat{r}$ is maximally complementary. $\square$