# The contact problem in Lagrangian systems with redundant frictional bilateral and unilateral constraints and singular mass matrix. The all-sticking contacts problem

Bernard Brogliato, Jozsef Kovecses, Vincent Acary

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

## ► To cite this version:

Bernard Brogliato, Jozsef Kovecses, Vincent Acary. The contact problem in Lagrangian systems with redundant frictional bilateral and unilateral constraints and singular mass matrix. The all-sticking contacts problem. Multibody System Dynamics, 2020, 48 (2), pp.151-192. 10.1007/s11044-019-09712-1. hal-04847393

HAL Id: hal-04847393

https://inria.hal.science/hal-04847393v1

Submitted on 19 Dec 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The contact problem in Lagrangian systems with redundant frictional bilateral and unilateral constraints and singular mass matrix. The all-sticking contacts problem

Bernard Brogliato

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France. bernard.brogliato@inria.fr

Jozsef Kovecses

Mechanical Engineering Department, Macdonald Engineering Building, Montreal, Canada. jozsef.kovecses@mcgill.ca

Vincent Acary

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France. vincent.acary@inria.fr

December 19, 2024

## Abstract

In this article we analyze the following problem: given a mechanical system subject to (possibly redundant) bilateral and unilateral constraints with set-valued Coulomb's friction, provide conditions such that the state which consists of all contacts sticking in both tangential and normal directions, is solvable. The analysis uses complementarity problems, variational inequalities, and linear algebra, hence it provides criteria which are, in principle, numerically tractable. An algorithm and several illustrating examples are proposed.

# 1 Introduction

The contact problem in constrained multibody systems, with or without friction, consists of analysing the system's behaviour under the action of external and contact wrenches. It has considerable practical interest in robotics tasks like biped locomotion stability analysis [47, 22, 16, 9, 48, 49], and dexterous object manipulation, grasping and workpiece fixture design or static equilibrium stability [62, 14, 68, 52, 53, 33, 10, 34, 51], to cite a few articles. The case with both unilateral and bilateral constraints, without reduction of the degrees of freedom, has not been as much studied as the pure bilaterally constrained, or unilaterally constrained problems, despite the fact that it represents the
most frequently met case in practice and industry. See [15, 17, 19, 20, 28, 44, 46, 65] for various
results in this setting, including frictionless and frictional constraints. Continuous-time contact
complementarity problems have been formulated in [1, 31, 30, 40] [29, Chapter 10], for various
cases of systems subjected to unilateral constraints, with and without friction. The analysis of such
problems has been led in [7, 6, 11, 42, 43, 54, 53, 70, 69]. The contact problem is formulated as a
mixed nonlinear complementarity problem in both frictionless and 3D frictional cases in [69], using
suitable slack variables. It is formulated as a quasi variational inequality in [54]. Sticking, sliding,
and mixed sticking/sliding cases are analysed. The crucial property for existence of solutions stems
from Theorem [25, Theorem 3.8.6] (see Theorem 1 in Appendix). In the frictional case, existence
results are stated for small enough coefficients of friction for the sliding modes, the upperbound
being however not provided explicitly. In the sliding/sticking mode case, Proposition 13 in [54]
states feasibility conditions. One important hypothesis in these articles is that Coulomb’s friction
holds at the acceleration level [19, Section 5.3.4]. Theorem 3.1 in [70] states the existence of normal
and tangential accelerations and multipliers, under the condition that the Delassus’ matrix is full
rank. Smallness of friction for sliding contacts guarantees existence, while smallness of friction at
sticking contacts guarantees uniqueness. In the all-sticking contacts case, existence holds under an
inequality involving local velocities and accelerations. Frictional unilateral constraints are considered
in [54, 53, 69, 70]. Lemma 2 and Theorem 2 in [42] are based on the assumption that matrices of
the contact problems are P-matrices ( i.e., with all principal minors positive), but no conditions on
the friction coefficients are given to guarantee this. The notion of *disassemblability* is introduced in
[7] to analyse systems with bilateral and unilateral constraints, with friction. Frictionless systems
are always disassemblable. The contact complementarity problem with friction is solvable if the
frictionless system is disassemblable and friction is small enough. No explicit upperbound on friction
is given. A contact MLCP is constructed in [6]. Sliding frictional contacts are analysed in [38] for
bilateral constraints. Singularities in the contact force are studied numerically on a three-body system
with one sliding joint. Force and form closure, as well as various stability notions and criteria are
proposed in [47, 22, 16, 9, 48, 49, 62, 68, 52, 53, 33, 10, 34, 51] for object grasping and biped robots.
Static equilibrium problems are analysed and algorithms are presented [51, 52, 33, 68, 49, 16, 48, 9],
dynamical effects are taken into account in [47, 62, 34]. Most of these works rely on the careful
study of contact and external wrenches, and yield precise criteria on wrenches and contact points
arrangements so that stability or static equilibrium are guaranteed [10, 48, 49].

The contact model that is chosen in this article is a rigid-body nonsmooth approach with bilateral
constraints, unilateral contacts and complementarity conditions (called hard-finger contact in the
dexterous manipulation literature), and set-valued friction. Compliant contact models can be used
also (linear or nonlinear spring-dashpot models). This may be needed in dexterous manipulation
tasks when flexible soft fingers are designed to manipulate fragile objects. Choosing a compliant
model among the various existing ones may not be easy (especially because of the difficulty related to
dissipation modeling, and how to choose the right dissipation model). Also one is often led to choose
artificially high damping in order to suppress unwanted oscillations, and artificially low equivalent
contact stiffness to avoid the integration of stiff dynamics which leads to long simulation times. In
some applications, the choice between both approaches (rigid or compliant contact) may not be
straightforward. It is nevertheless important to remind that flexible contacts allowing for detachment
("unilateral springs"), can be recast into a complementarity framework quite similar to the one
developed in this article, see [51, 19]. Compliance at contacts can also sometimes remove the non-uniqueness of contact forces (due to constraints redundancy), though this may fail in some problems where Coulomb friction set-valuedness creates the non-uniqueness issue [19, section 5.5.6]. Then some regularization of Coulomb's model have to be used (henceby destroying the sticking modes), but this creates new problems of robustness and high sensitivity of solutions with respect to parameter uncertainties (see [74] for bilaterally constrained systems, where it is shown that regularized models like the LuGre model and its variants are quite inappropriate).

When using the so-called natural coordinates [37], the mass matrix of a multibody system may become constant but singular, and nonlinear bilateral constraints are necessarily added to the dynamics. It is then important to get criteria of well-posedness of the contact problem for singular mass matrices, and redundant constraints. This has been done in [36] when only frictionless bilateral constraints are present (hence giving rise to the classical KKT conditions for the index one formulation of bilateral constraints [19, Section 5.1.1] [39, 60, 63]), and in [20] for both frictionless bilateral and unilateral constraints. When all contacts are sliding, the contact problem has been studied in [15] in various cases (bilateral and/or unilateral contacts, with or without friction), and with a detailed analysis of the index one KKT conditions.

In this article we analyse the contact problem when the multibody system is subjected to both bilateral and unilateral constraints with set-valued Coulomb's friction. We focus in particular on the case of all-sticking contacts, i.e., contacts are activated and are not slipping at the time instant where the problem is analysed. The presented criteria do not allow one to know whether or not the system attains such an all-sticking mode, however it provides conditions that guarantee that such modes are well-posed (with existence and possibly uniqueness of contact forces and acceleration). Moreover the results on the well-posedness of a class of variational inequalities of the second kind presented in [3] are used to formulate conditions for the well-posedness of the all-sticking contact problem when the mass matrix is singular and constraints are redundant.

## 1.1 Generic problem

In a Lagrangian formalism such systems may be written generically as follows,

$$ M(q)\ddot{q} + F(q, \dot{q}, t) = \nabla h_{n,b}(q)\lambda_{n,b} + \nabla h_{n,u}(q)\lambda_{n,u} + H_{t,b}(q)\lambda_{t,b} + H_{t,u}(q)\lambda_{t,u} \quad (1a) $$

$$ \text{Local friction law (Coulomb)} : (v_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i), \quad 1 \le i \le m \quad (1b) $$

$$ \text{Complementarity conditions} : 0 \le h_{n,u}(q) \perp \lambda_{n,u} \ge 0 \quad (1c) $$

$$ \text{Bilateral (holonomic) constraints} : h_{n,b}(q) = 0, \quad (1d) $$

where $q \in \mathbb{R}^n$ is the vector collecting the generalized coordinates $q_i$, $1 \le i \le n$, assumed to be independent when all the constraints are removed, $\dot{q}$ is the vector of generalized velocities, $M(q) = M(q)^\top$ is the symmetric positive definite (or only semi definite) inertia matrix, $F(q, \dot{q}, t)$ collects internal forces (including forces deriving from a potential as well as Coriolis and centripetal forces), and external actions on the system such as disturbances or control. A more general dynamical framework could be chosen with Newton-Euler dynamics, where generalized velocities $\dot{q}$ are replaced by $v = T(q)\dot{q}$ for some transformation matrix $T(q)$, and $v$ may contain quantities like instantaneous angular velocities (so that $T(q)$ involves Olinde Rodrigues matrices).
We consider $m = m_u+m_b$ constraints (or contacts) consisting of $m_u$ unilateral (inequality) constraints $h_{n,u}(q) \in \mathbb{R}^{m_u}$ and $m_b$ bilateral (equality) constraints $h_{n,b}(q) \in \mathbb{R}^{m_b}$. The vectors of Lagrange multipliers $\lambda_{n,u} \in \mathbb{R}^{m_u}$ and $\lambda_{n,b} \in \mathbb{R}^{m_b}$ are associated with the unilateral and bilateral constraints, respectively. We assume that the gap functions $h_{n,u}(q)$ are defined from the local kinematics at the contact points [2, 19]. Thus $\lambda_{n,u}$ corresponds to the normal components of the contact forces. Let $d \in \mathbb{N}$ be the ambient space dimension ($d=2$ in the planar case and $d=3$ in the three dimensional case). The tangential components of the contact forces are, in turn, collected by $\lambda_{t,u} \in \mathbb{R}^{(d-1)m_u}$ (unilateral constraints) and $\lambda_{t,b} \in \mathbb{R}^{(d-1)m_b}$ (bilateral constraints), with $H_{t,b}(q) \in \mathbb{R}^{n \times (d-1)m_b}$ and $H_{t,u}(q) \in \mathbb{R}^{n \times (d-1)m_u}$ linear operators that map local tangent frames at the contact points to the generalized coordinates, obtained from the so-called contact kinematics and local frames [2, Section 3.3] [19, Chapter 4] [57, Chapter 4] [56]. Note that the size of $\lambda_{t,b}$, $\lambda_{t,u}$, $H_{t,b}(q)$ and $H_{t,u}(q)$ depends on the space dimension $d$ considered.

The unilateral constraints and their associated Lagrange multipliers are related through the complementarity condition (1c), which is to be understood componentwise (per contact). It models the fact that for each contact $i$, the normal contact force should not act at a unilateral contact point if the contact is open (i.e., $h_{n,u,i}(q) > 0$), and that $\lambda_{n,u,i} > 0$ if and only if $h_{n,u,i}(q) = 0$. This type of complementarity conditions does not encompass adhesive contact forces [19, Section 5.4.1], which in turn may be useful in some cases [61]. We note anyway that adhesive effects could be taken into account while remaining in a set-valued framework [2, Section 3.9.4.4], however this is outside the scope of this paper.

The tangential counterpart of the $i^{th}$ contact force is, in turn, driven by the Coulomb friction law (1b). Let us denote $\lambda_i \in \mathbb{R}^d$ the $i^{th}$ contact force, expressed in a local frame attached to the contact point, and $v_i$ the local space velocity at contact. The $\lambda_{u,i}$ vector collects the $i^{th}$ scalar component of the vector $\lambda_{n,u}$ (normal part), together with the $i^{th}$ block component (of size $d-1$) of the vector $\lambda_{t,u}$ (tangential part), i.e. $\lambda_{u,i} = (\lambda_{n,u,i}, \lambda_{t,u,i})^\top$. Similarly for $\lambda_{b,i} = (\lambda_{n,b,i}, \lambda_{t,b,i})^\top$. At each contact point $i$, the Coulomb friction law ($v_i, \lambda_i \in \mathcal{C}(\mathbf{n}_i, \mu_i)$ relates the contact force $\lambda_i$ to the local velocity $v_i$ through a nonsmooth multivalued law involving the normal vector $\mathbf{n}_i$ at contact and the coefficient of friction $\mu_i$. For the sake of clarity, the Coulomb friction law is fully described in Section 1.2. Finally we let $\lambda_{n,u} = (\lambda_{n,u,1}, \dots, \lambda_{n,u,m_u})^\top \in \mathbb{R}^{m_u}$, $\lambda_{n,b} = (\lambda_{n,b,1}, \dots, \lambda_{n,b,m_b})^\top \in \mathbb{R}^{m_b}$, $\lambda_{t,u} = (\lambda_{t,u,1}, \dots, \lambda_{t,u,m_u})^\top \in \mathbb{R}^{(d-1)m_u}$, $\lambda_{t,b} = (\lambda_{t,b,1}, \dots, \lambda_{t,b,m_b})^\top \in \mathbb{R}^{(d-1)m_b}$.

**Remark 1** An underlying assumption, is that all joints (bilateral as well as unilateral ones) can be modelled with a finite number of contact points. Thus Coulomb’s law can be applied at each contact point. This may require a modelling step to properly represent, for instance, bilateral constraints that stem from revolute joints in a mechanism, see e.g. [5].

**Mathematical notation:** A (possibly non symmetric) matrix $M \in \mathbb{R}^{n \times n}$ is positive semidefinite ($M \succ 0$) if $x^\top M x \ge 0$ for all $x \in \mathbb{R}^n$. It is positive definite ($M \succ 0$) if $x^\top M x > 0$ for all $x \ne 0$. We denote $M \succ N$ if $M - N \succ 0$. For a matrix $M$, we denote $M_{i,\bullet}$ its $i^{th}$ row, hence $M_{i,\bullet}x$ is the $i^{th}$ component of the vector $Mx$. Similarly $M_{\bullet,j}$ is its $j^{th}$ column. The generalized (Moore-Penrose) inverse of $M$ is denoted $M^\dagger$. The image of $M$ is denoted $\text{im}(M)$ and its kernel $\ker(M)$. The identity matrix with suitable dimension is denoted $I$. Let the $n \times n$ matrix $M \succ 0$ be symmetric, then it defines an inner product as $\langle x, y \rangle = x^\top M y$ for any $x$ and $y \in \mathbb{R}^n$. With this
metric, the orthogonal projection of a vector $x \in \mathbb{R}^n$ on a closed convex set $K \subset \mathbb{R}^n$ is denoted as $\text{proj}_M[K; x] = \text{argmin}_{z \in K} \frac{1}{2}(z-x)^\top M(z-x)$. Let $x \in \mathbb{R}^n$ and $y = Mx + r \in \mathbb{R}^n$; a Linear Complementarity Problem (LCP) is a problem of the form: $x \ge 0, y \ge 0, x^\top y = 0$, written compactly as $0 \le x \perp y \ge 0$. Let $x \in \mathbb{R}^n, z \in \mathbb{R}^m, w = Cx + Dz + b \in \mathbb{R}^m$; a Mixed LCP (MLCP) is a problem of the form: $x \ge 0, w \ge 0, w^\top z = 0, Ax + Bz + a = 0$, for some matrices $A$ and $B$ and a vector $q$. We denote $0_n = (0, \dots, 0)^\top \in \mathbb{R}^n$, and $0^n = (0_n)^\top$. The polyhedral set $P(A, b) \stackrel{\Delta}{=} \{x \in \mathbb{R}^n | Ax \ge b\}$. Let $A \in \mathbb{R}^{n \times n}$ be a real square matrix, its induced norm is $\|A\|_2 = \max_{x \in \mathbb{R}^n, \|x\|_2=1} \|Ax\|_2$, where $\|x\|_2$ is the Euclidean norm on the vector space $\mathbb{R}^n$. This induced matricial norm is sometimes denoted as $\|A\|_{2,2}$. One has $\|A\|_2 = \sigma_{\max}(A) = \sqrt{\lambda_{\max}(AA^\top)}$ [12, Proposition 9.4.7], where $\sigma_{\max}$ is the maximum singular value of $A$, and $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$ are its smallest and the largest eigenvalues, respectively [12, Proposition 9.4.7]. Moreover, if $A$ is invertible, $\sigma_{\min}(A) = \frac{1}{\sigma_{\max}(A^{-1})}$, where $\sigma_{\min}(A)$ is the smallest singular value of $A$ [12, Fact 6.3.21].

## 1.2 Formulation of the Coulomb friction law

The Coulomb friction law is defined locally at each contact point of the system. At contact $i$, let us assume that the geometry of the surface of contact is smooth enough, so that we can define a tangent plane of contact $S$ and the normal (unitary) vector $\mathbf{n}_i$ at $S$. The vector $\mathbf{n}_i$, together with two arbitrary orthonormal vectors $\mathbf{t}_{1,i}$, $\mathbf{t}_{2,i}$ lying onto $S$, define a local basis $\mathcal{B}_i$. We assume that each contact $i$ involves only two parts of the whole mechanical system, and denote by $v_i$ the local (space) relative velocity between the two parts in contact, which can be decomposed in $\mathcal{B}_i$ as a normal part $v_{n,i} = v_i^\top \mathbf{n}_i \in \mathbb{R}$ and a tangential part $v_{t,i} = \{v_i^\top \mathbf{t}_{1,i}, v_i^\top \mathbf{t}_{2,i}\} \in \mathbb{R}^2$ (when the contact $i$ is closed, then $v_{n,i} = 0$). As mentioned earlier, the contact force $\lambda_i$ can be similarly decomposed onto its normal part $\lambda_{n,i} \in \mathbb{R}$ and its tangential part $\lambda_{t,i} \in \mathbb{R}^2$. Let $\mu_i$ be the coefficient of friction at the contact point. It is assumed that the friction coefficient is constant, though nonlinear Stribeck effects could be added in our framework, see [19, section 5.3.1].

Let us assume that the contact $i$ is active (this trivially holds for bilateral contacts). The Coulomb friction law, symbolically denoted as $(v_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i)$, relates the normal and tangential parts of both the local velocity $v_i$ and the local contact force $\lambda_i$ according to the following formulations [19, §5.3] [2, §3.9.1] [41, §5.3.2] [45]:

$$ (v_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i) \iff \left\{ \begin{array}{ll} \text{either} & \|\lambda_{t,i}\| \le \mu_i |\lambda_{n,i}| \text{ and } v_i = 0 \quad \text{(sticking mode)} \\ \text{or} & \|\lambda_{t,i}\| = \mu_i |\lambda_{n,i}| \text{ and } v_{n,i} = 0, v_{t,i} \ne 0 \quad (2) \\ & \text{and } \exists \alpha_i > 0, \lambda_{t,i} = -\alpha_i v_{t,i} \quad \text{(sliding mode).} \end{array} \right. $$

In the sliding mode, one equivalently has $\lambda_{t,i} = -\mu_i |\lambda_{n,i}| \frac{v_{t,i}}{|v_{t,i}|}$. Let us recall that $(v_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i) \Rightarrow \lambda_{t,i} \in \mathcal{D}(\lambda_{n,i}, \mu_i)$, where $\mathcal{D}(\lambda_{n,i}, \mu_i) = |\lambda_{n,i}|\mathcal{D}_{\mu_i}$ is the Coulomb-Moreau's disk [45] [19, Section 5.3.2], and $\mathcal{D}_{\mu_i} = \{z \in \mathbb{R}^{d-1} | z^\top z \le \mu_i^2\}$. Then Coulomb's law (in velocity) in (2) is equivalently formulated as the variational inequality: Find $\lambda_{t,i} \in \mathcal{D}(\lambda_{n,i}, \mu_i)$ such that

$$ v_{t,i}^\top (y - \lambda_{t,i}) \ge 0 \quad \text{for all } y \in \mathcal{D}(\lambda_{n,i}, \mu_i), \quad (3) $$

which shows that Coulomb's friction follows a maximal dissipation "principle". This variational inequality is in turn equivalent to:

$$ v_{t,i} \in -\partial\psi_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(\lambda_{t,i}) \Leftrightarrow \lambda_{t,i} = \text{proj}[\mathcal{D}(\lambda_{n,i}, \mu_i); \lambda_{t,i} - \rho v_{t,i}] \quad \text{for any } \rho > 0. \quad (4) $$
See Appendix C for the definition of the indicator function of a closed convex set. Notice from (4) that $\lambda_{t,i} \in \text{Int}(\mathcal{D}(\lambda_{n,i}, \mu_i))$ implies that $v_{t,i} = 0$. It is noteworthy that other friction models may be cast into such nonsmooth, set-valued framework (Coulomb-Moreau's disk may be replaced by other convex sets, allowing for anisotropic friction [41, Fig. 5.6]), as well as Stribeck effects in slipping phases [19, Remark 5.16], micro displacements occurring during the sticking phase [2, Section 3.9.4.2], and more elaborated models [40, 66]. A crucial feature is that sticking contacts can be correctly handled when the tangential velocity vanishes. This can be done only with a model of friction that is set-valued at zero tangential velocity.

**Restriction to the 2D case** In the planar case, the tangent space $\mathcal{S}$ to the surface of contact is a line and the local contact basis $\mathcal{B}_i$ is defined by two vectors, the normal $\mathbf{n}_i$ and the tangent $\mathbf{t}_i$ at contact. The two cases of the Coulomb law boil down to a single, compact formulation relating the tangential components of the force and the velocity,

$$ \lambda_{t,i} \in -\mu_i |\lambda_{n,i}| \operatorname{sgn}(v_{t,i}), $$

where $\operatorname{sgn}(x)$ is the set-valued signum function: $\operatorname{sgn}(x) = 1$ if $x > 0$, $\operatorname{sgn}(x) = -1$ if $x < 0$, $\operatorname{sgn}(0) = [-1, 1]$.

**Coulomb's law at the acceleration level** So far we have formulated the Coulomb friction law at the velocity level, regardless of the value of the acceleration. When studying the transition between the sticking case and the sliding case, it is however helpful to make a further assumption pertaining to the value of the acceleration $a_i = \dot{v}_i$. In particular, a common assumption consists in imposing the force $\lambda_i$ to lie on the border of the cone in the sticking case ($v_i = 0$) as soon as the acceleration $a_i$ ceases to vanish in the tangential direction [29, 53]. In this case the tangential component $\lambda_{t,i}$ of the force should be parallel and opposed in sign to $a_{t,i}$. This new model is symbolically denoted as $(v_i, a_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i)$ and summarized as [19, §5.3.4]:

$$ (v_i, a_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i) \iff \left\{ \begin{array}{ll} \text{either} & \begin{cases} \text{either } \|\lambda_{t,i}\| \le \mu_i |\lambda_{n,i}| & \text{and } v_i = 0, a_i = 0 \\ \text{or } \|\lambda_{t,i}\| = \mu_i |\lambda_{n,i}| & \text{and } v_i = 0, a_{n,i} = 0, a_{t,i} \ne 0 \\ \text{and } \exists \beta_i > 0, \lambda_{t,i} = -\beta_i a_{t,i} & \text{(sticking modes)} \end{cases} \\ \text{or} & \begin{cases} \|\lambda_{t,i}\| = \mu_i |\lambda_{n,i}| & \text{and } v_{n,i} = 0, v_{t,i} \ne 0 \\ \text{and } \exists \alpha_i > 0, \lambda_{t,i} = -\alpha_i v_{t,i} & \text{(sliding mode).} \end{cases} \end{array} \right. \tag{5} $$

The second sticking mode yields an equivalent formulation $\lambda_{t,i} = -\mu_i |\lambda_{t,i}| \frac{a_{t,i}}{\|a_{t,i}\|}$, while the sliding mode yields the same $\alpha_i$ as in (2). Note that $(v_i, a_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i) \Rightarrow (v_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i)$. Compared to the previous model formulated at the velocity level, a new subcase of the sticking mode has been added here, which corresponds to a stick→slip transition (vanishing velocity and non-vanishing tangential acceleration $a_{t,i}$). We may rewrite the first part of the model in (5), named “sticking modes”, as follows [64, Equations (5) (6)] : If $v_i = 0$ and $a_{n,i} = 0$, find $\lambda_{t,i} \in \mathcal{D}(\lambda_{n,i}, \mu_i)$ such that:

$$ a_{t,i}^{\top}(y - \lambda_{t,i}) \ge 0 \text{ for all } y \in \mathcal{D}(\lambda_{n,i}, \mu_i). \tag{6} $$
This can in turn be formulated equivalently as: If $v_i = 0$ and $a_{n,i} = 0$,

$$ 
\begin{aligned}
a_{t,i} &\in -\partial\psi_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(\lambda_{t,i}) &\Leftrightarrow \lambda_{t,i} &= \operatorname{proj}[\mathcal{D}(\lambda_{n,i}, \mu_i); \lambda_{t,i} - \rho a_{t,i}] \text{ for any } \rho > 0 \\
&\Leftrightarrow \lambda_{t,i} &\in \partial\psi_{\mathcal{D}(\lambda_{n,i}, \mu_i)}^*(a_{t,i}),
\end{aligned}
\quad (7) $$

where $\psi_{\mathcal{D}(\lambda_{n,i}, \mu_i)}^*(\cdot)$ is the conjugate function of $\psi_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(\cdot)$, i.e., the support function of the set $\mathcal{D}(\lambda_{n,i}, \mu_i)$. All the quantities in the above are considered at a given time $t$, and it is of interest to investigate their evolution on the right of $t$. Notice that if $a_{n,i}(t^+) > 0$, then $v_{n,i}(t^+) > 0$, $\lambda_{n,i}(t^+) = 0$ and it follows from (7) that $\lambda_{t,i}(t^+) = 0$ and $a_{t,i}(t^+) \in \mathbb{R}$, the latter being irrelevant in such a case where the contact opens.

## 1.3 Mechanical notation

Since contacts may be frictionless, or in contrast may involve Coulomb's friction, which in turn may be in a sliding or sticking mode, the following conventions shall be adopted:

* Contacts: $m_b$ bilateral contacts, $m_u$ unilateral contacts, $m_b + m_u = m$.

* Bilateral contacts: $1 \le i \le m_b$, i.e., $i \in \mathcal{I}_b$.

- Frictional bilateral contacts: $1 \le i \le m_b^\mu$, i.e., $i \in \mathcal{I}_b^\mu$.

* Sticking frictional bilateral contacts: $1 \le i \le m_b^{\mu, st}$, i.e., $i \in \mathcal{I}_b^{\mu, st}$.

* Sliding frictional bilateral contacts: $m_b^{\mu, st} + 1 \le i \le m_b^\mu$, i.e., $i \in \mathcal{I}_b^{\mu, st}$.

- Frictionless bilateral contacts: $m_b^\mu + 1 \le i \le m_b$, i.e., $i \in \mathcal{I}_b^0$.

* Unilateral contacts: $m_b + 1 \le i \le m$, i.e., $i \in \mathcal{I}_u$.

- Frictional unilateral contacts: $m_b + 1 \le i \le m_u^\mu$, i.e., $i \in \mathcal{I}_u^\mu$.

* Sticking frictional unilateral contacts: $m_b + 1 \le i \le m_u^{\mu, st}$, i.e., $i \in \mathcal{I}_u^{\mu, st}$.

* Sliding frictional unilateral contacts: $m_u^{\mu, st} + 1 \le i \le m_u^\mu$, i.e., $i \in \mathcal{I}_u^{\mu, st}$.

- Frictionless unilateral contacts: $m_u^\mu + 1 \le i \le m$, i.e., $i \in \mathcal{I}_u^0$.

We thus define ten disjoint index sets with $\mathcal{I}_b \cup \mathcal{I}_u = (\mathcal{I}_b^\mu \cup \mathcal{I}_b^0) \cup (\mathcal{I}_u^\mu \cup \mathcal{I}_u^0) = ((\mathcal{I}_b^{\mu, st} \cup \mathcal{I}_b^{\mu, sl}) \cup \mathcal{I}_b^0) \cup ((\mathcal{I}_u^{\mu, st} \cup \mathcal{I}_u^{\mu, sl}) \cup \mathcal{I}_u^0)$, which correspond to the different modes of the system. We may therefore rewrite (1a) as

$$ 
\begin{aligned}
M(q)\ddot{q} + F(q, \dot{q}, t) &= \sum_{i \in \mathcal{I}_b^{\mu, st}} \nabla h_{n,b,i}^{st}(q) \lambda_{n,b,i}^{st} + \sum_{i \in \mathcal{I}_b^{\mu, st}} \nabla h_{n,b,i}^{sl}(q) \lambda_{n,b,i}^{sl} + \sum_{i \in \mathcal{I}_b^0} \nabla h_{n,b,i}^0(q) \lambda_{n,b,i}^0 \\
&\quad + \sum_{i \in \mathcal{I}_b^{\mu, st}} H_{t,b,i}^{st}(q) \lambda_{t,b,i}^{st} + \sum_{i \in \mathcal{I}_b^{\mu, st}} H_{t,b,i}^{sl}(q) \lambda_{t,b,i}^{sl} \\
&\quad + \sum_{i \in \mathcal{I}_u^{\mu, st}} \nabla h_{n,u,i}^{st}(q) \lambda_{n,u,i}^{st} + \sum_{i \in \mathcal{I}_u^{\mu, st}} \nabla h_{n,u,i}^{sl}(q) \lambda_{n,u,i}^{sl} + \sum_{i \in \mathcal{I}_u^0} \nabla h_{n,u,i}^0(q) \lambda_{n,u,i}^0 \\
&\quad + \sum_{i \in \mathcal{I}_u^{\mu, st}} H_{t,u,i}^{st}(q) \lambda_{t,u,i}^{st} + \sum_{i \in \mathcal{I}_u^{\mu, st}} H_{t,u,i}^{sl}(q) \lambda_{t,u,i}^{sl},
\end{aligned}
\quad (8) $$

where $\sum_{i \in \mathcal{I}_b^{\mu, st}} \nabla h_{n,b,i}^{st}(q) \lambda_{n,b,i}^{st} = \nabla h_{n,b}^{st}(q) \lambda_{n,b}^{st}$, and so on for the remaining terms.
## 1.4 Systems with 3D frictional persistent sticking contacts

Let us assume that $\mathcal{I}_b^0 = \mathcal{I}_u^0 = \mathcal{I}_b^{\mu,sl} = \mathcal{I}_u^{\mu,sl} = \emptyset$, which implies $\mathcal{I}_b^{\mu,sl} \neq \emptyset$ and $\mathcal{I}_u^{\mu,sl} \neq \emptyset$ and $\mathcal{I}_b^{\mu,sl} \cup \mathcal{I}_u^{\mu,sl} = \mathcal{I}_u \cup \mathcal{I}_b$: all contacts are supposed to be in the sticking mode in (2). An active contact point $i$ is in a persistent sticking mode on a non-trivial interval $[t_0, t]$ if and only if the non-holonomic constraint<fnref n="1" /> $H_{t,i}(q)^\top \dot{q} = v_{t,i}(t) = 0$ is satisfied for all $t \in [t_0, t]$. We deduce that $a_{t,i} = H_{t,i}(q)^\top \ddot{q} + \frac{d}{dt}(H_{t,i}(q)^\top)\dot{q} = 0$ on $[t_0, t]$, where $\ddot{q}(\cdot)$ is assumed to be right-continuous, i.e. $\ddot{q}(t) = \ddot{q}(t^+)$ (and same for the contact force multipliers). Since it is supposed that all the contacts are frictional and in sticking mode, it follows that $H_{t,b}^{st}(q)^\top \ddot{q} + \frac{d}{dt}(H_{t,b}^{st}(q)^\top)\dot{q} = 0$ and $H_{t,u}^{st}(q)^\top \ddot{q} + \frac{d}{dt}(H_{t,u}^{st}(q)^\top)\dot{q} = 0$ in a right neighborhood of $t$. Our objective is to analyse the conditions under which the system with all persistent sticking contacts is solvable. In other words, such that the following set of equations has at least one solution (the unknowns being the acceleration and the contact force multipliers):

$$ 
\left\{
\begin{array}{l}
(a) \ M(q)\ddot{q} + F(q, \dot{q}, t) = \nabla h_{n,b}^{st}(q)\lambda_{n,b}^{st} + \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} + H_{t,b}^{st}(q)\lambda_{t,b}^{st} + H_{t,u}^{st}(q)\lambda_{t,u}^{st} \\
(b) \ a_{t,b} = H_{t,b}^{st}(q)^\top \ddot{q} + \frac{d}{dt}(H_{t,b}^{st}(q)^\top)\dot{q} = 0 \\
(c) \ a_{t,u} = H_{t,u}^{st}(q)^\top \ddot{q} + \frac{d}{dt}(H_{t,u}^{st}(q)^\top)\dot{q} = 0 \\
(d) \ a_{n,b} = \nabla h_{n,b}^{st}(q)^\top \ddot{q} + \frac{d}{dt}(\nabla h_{n,b}^{st}(q)^\top)\dot{q} = 0 \\
(e) \ 0 \le \lambda_{n,u}^{st} \perp a_{n,u} = \nabla h_{n,u}^{st}(q)^\top \ddot{q} + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^\top)\dot{q} \ge 0,
\end{array}
\right.
\quad (9) $$

where the time argument has been dropped. The set of equations in (9) forms a mixed LCP. Indeed if the system (9) has a solution that respects Coulomb's constraint, the system persists in the all-sticking mode in a right-neighborhood of $t$. Assume that $M(q) \succ 0$. Using the first line to express the acceleration one obtains:

$$ 
\left\{
\begin{array}{l}
(a) \ M(q)\ddot{q} + F(q, \dot{q}, t) = \nabla h_{n,b}^{st}(q)\lambda_{n,b}^{st} + \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} + H_{t,b}^{st}(q)\lambda_{t,b}^{st} + H_{t,u}^{st}(q)\lambda_{t,u}^{st} \\
(b) \ a_{t,b} = H_{t,b}^{st}(q)^\top M(q)^{-1} H_{t,b}^{st}(q)\lambda_{t,b}^{st} + R_{t,b}^{st}(q, \dot{q}, t) = 0 \\
(c) \ a_{t,u} = H_{t,u}^{st}(q)^\top M(q)^{-1} H_{t,u}^{st}(q)\lambda_{t,u}^{st} + R_{t,u}^{st}(q, \dot{q}, t) = 0 \\
(d) \ a_{n,b} = \nabla h_{n,b}^{st}(q)^\top M(q)^{-1} \nabla h_{n,b}^{st}(q)\lambda_{n,b}^{st} + R_{n,b}^{st}(q, \dot{q}, t) = 0 \\
(e) \ 0 \le \lambda_{n,u}^{st} \perp a_{n,u} = \nabla h_{n,u}^{st}(q)^\top \ddot{q} + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^\top)\dot{q} \ge 0,
\end{array}
\right.
\quad (10) $$

which is a mixed LCP with unknowns $\ddot{q}$, $\lambda_{n,b}^{st}$, $\lambda_{n,u}^{st}$, $\lambda_{t,b}^{st}$, $\lambda_{t,u}^{st}$, where:

$$ R_{t,b}^{st}(q, \dot{q}, t) = H_{t,b}^{st}(q)^\top M(q)^{-1} [-F(q, \dot{q}, t) + \nabla h_{n,b}^{st}(q)\lambda_{n,b}^{st} + \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} + H_{t,u}^{st}(q)\lambda_{t,u}^{st}] \\ + \frac{d}{dt}(H_{t,b}^{st}(q)^\top)\dot{q}, \quad (11) $$

<footnote n="1">¹This is sometimes called rolling friction, or friction without sliding. The expression for the velocities comes from the principle of virtual work.</footnote>
$$ R_{t,u}^{st}(q, \dot{q}, t) = H_{t,u}^{st}(q)^{\top} M(q)^{-1} [-F(q, \dot{q}, t) + \nabla h_{n,b}^{st}(q) \lambda_{n,b}^{st} + \nabla h_{n,u}^{st}(q) \lambda_{n,u}^{st} + H_{t,b}^{st}(q) \lambda_{t,b}^{st}] \\ + \frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q}, \quad (12) $$

$$ R_{n,b}^{st}(q, \dot{q}, t) = \nabla h_{n,b}^{st}(q)^{\top} M(q)^{-1} [-F(q, \dot{q}, t) + \nabla h_{n,u}^{st}(q) \lambda_{n,u}^{st} + H_{t,b}^{st}(q) \lambda_{t,b}^{st} + H_{t,u}^{st}(q) \lambda_{t,u}^{st}] \\ + \frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top})\dot{q}. \quad (13) $$

## 2 Bilaterally constrained systems

Let us deal briefly with the bilaterally constrained case. Assume that $m_u = 0$, from (9) (a) (b) (d) one obtains the system:

$$ (a) \underbrace{\begin{pmatrix} M(q) & -\nabla h_{n,b}^{st}(q) & -H_{t,b}^{st}(q) \\ \nabla h_{n,b}^{st}(q)^{\top} & 0 & 0 \\ H_{t,b}^{st}(q)^{\top} & 0 & 0 \end{pmatrix}}_{\triangleq M_b(q) \in \mathbb{R}^{(n+dm) \times (n+dm)}} \underbrace{\begin{pmatrix} \ddot{q} \\ \lambda_{n,b}^{st} \\ \lambda_{t,b}^{st} \end{pmatrix}}_{\triangleq z_b} = \underbrace{\begin{pmatrix} -F(q, \dot{q}, t) \\ -\frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top})\dot{q} \\ -\frac{d}{dt}(H_{t,b}^{st}(q)^{\top})\dot{q} \end{pmatrix}}_{\triangleq F_b(q, \dot{q}, t) \in \mathbb{R}^{n+dm}} \quad (14) $$

$$ (b) \quad \lambda_{t,b,i}^{st} \in \mathcal{D}(\lambda_{n,b,i}^{st}, \mu_i) \text{ for all } 1 \le i \le m. $$

Part (a) possesses the form of a classical KKT system which is ubiquitous in bilaterally constrained systems [63, 74, 27, 76] and whose well-posedness conditions have been deeply analysed in [15, Appendix E] and [36]. We may rewrite (14) (a) equivalently as an optimization problem with equality constraints, where $\lambda_{n,b}^{st}$ and $\lambda_{t,b}^{st}$ are multipliers associated with the constraints (9) (b) (e). The complete problem therefore involves constraints (14) (b) on the classical KKT system multipliers. Notice that (9) and (14) do not rely on the assumption that $M(q) > 0$. The next result relies on $M(q) \ge 0$ solely. It has been tackled in [36] and [15].

**Proposition 1 (i)** The system (14) (a) has a unique solution $z_b$ for any $F_b(q, \dot{q}, t)$ (equivalently $M_b(q)$ is invertible) if and only if $\text{rank}(\nabla h_{n,b}^{st}(q) H_{t,b}^{st}(q)) = dm_b$ ($\Rightarrow$ $dm_b \le n$), and

$$ \ker(M(q)) \cap \ker \begin{pmatrix} \nabla h_{n,b}^{st}(q)^{\top} \\ H_{t,b}^{st}(q)^{\top} \end{pmatrix} = \{0\}. \quad (15) $$

(ii) Consider the system (14) (a). The necessary and sufficient condition for the existence of $\ddot{q}$, $\lambda_{n,b}^{st}$ and $\lambda_{t,b}^{st}$, such that uniqueness of $\ddot{q}$, $\nabla h_{n,b}^{st}(q)\lambda_{n,b}$ and $H_{t,b}^{st}(q)\lambda_{t,b}$, holds for arbitrary $F(q, \dot{q}, t)$, is that

$$ -\frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top})\dot{q} \in \text{im}(\nabla h_{n,b}^{st}(q)^{\top}), \quad -\frac{d}{dt}(H_{t,b}^{st}(q)^{\top})\dot{q} \in \text{im}(H_{t,b}^{st}(q)^{\top}), \\ \ker(M(q)) \cap \ker \begin{pmatrix} \nabla h_{n,b}^{st}(q)^{\top} \\ H_{t,b}^{st}(q)^{\top} \end{pmatrix} = \{0\}. \quad (16) $$

**Proof 1** The proof of (i) can be established using [15, Appendix E, item (iii')], the proof of (ii) follows from [15, Appendix E, item (iv)].
The interest of item (ii) is that it does not restrict the number of constraints and the rank of the acceleration constraints matrices, and thus allows for an arbitrary number of contact points. In case (ii) the matrix $M_b(q)$ is not necessarily full-rank. But since a solution exists that solves $M_b(q)z_b = F_b(q, \dot{q}, t)$, this solution can be expressed as $z_b = M_b(q)^{\dagger}F_b(q, \dot{q}, t) + (I - M_b(q)^{\dagger}M_b(q))y$ for any $y \in \mathbb{R}^{n+dm}$, and $M_b(q)^{\dagger}$ is the generalized inverse of $M_b(q)$ [12, Proposition 6.1.7]. The vector $y$ reflects that there may be more contact multipliers than independent equations to calculate them.

## 2.1 Solvability of (14) (a) and (b) (independent constraints)

The next step is to incorporate (14) (b) in the analysis. For this we proceed with the explicit calculation of the solution to (14) (a), doing the assumptions which allow us to use the matrix inversion lemma [12, Proposition 2.8.7]. Let us assume that $M(q) > 0$, $\nabla h_{n,b}^{st}(q)^{\top}M(q)^{-1}\nabla h_{n,b}^{st}(q) > 0$, and denote

$$ G(q) \triangleq \begin{pmatrix} \nabla h_{n,b}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,b}^{st}(q) & \nabla h_{n,b}^{st}(q)^{\top} M(q)^{-1} H_{t,b}^{st}(q) \\ H_{t,b}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,b}^{st}(q) & H_{t,b}^{st}(q)^{\top} M(q)^{-1} H_{t,b}^{st}(q) \end{pmatrix} = \begin{pmatrix} G_{nn}(q) & G_{nt}(q) \\ G_{nt}(q)^{\top} & G_{tt}(q) \end{pmatrix}, \quad (17) $$

where $G(q) \in \mathbb{R}^{dm \times dm}$, and

$$ S_G(q) \triangleq G_{tt}(q) - G_{nt}(q)^{\top} G_{nn}(q)^{-1} G_{nt}(q) = S_G(q)^{\top} \in \mathbb{R}^{(d-1)m \times (d-1)m} \quad (18) $$

is the Schur complement of $G_{nn}(q)$ in $G(q)$. Then supposing that $G(q) > 0$ ($\implies S_G(q) > 0$ by the Schur complement Theorem), which implies that the gradients $\nabla h_{n,b,i}^{st}$, $1 \le i \le m$ are independent, as well as the operators $H_{n,b,i}^{st}$, $1 \le i \le dm$, the matrix $M_b(q)$ is invertible and the Banachiewicz-Schur form of its inverse is [67]<fnref n="2" />:

$$ M_b(q)^{-1} = \begin{pmatrix} M^{-1} - M^{-1}(\nabla h_{n,b}^{st} H_{t,b}^{st}) G^{-1} \begin{pmatrix} \nabla h_{n,b}^{st,\top} \\ H_{t,b}^{st,\top} \end{pmatrix} M^{-1} & M^{-1}(\nabla h_{n,b}^{st} H_{t,b}^{st}) G^{-1} \\ -G^{-1} \begin{pmatrix} \nabla h_{n,b}^{st,\top} \\ H_{t,b}^{st,\top} \end{pmatrix} M^{-1} & G^{-1} \end{pmatrix} \quad (19) $$

and

$$ G(q)^{-1} = \begin{pmatrix} G_{nn}^{-1} + G_{nn}^{-1} G_{nt} S_G^{-1} G_{nt}^{\top} G_{nn}^{-1} & -G_{nn}^{-1} G_{nt} S_G^{-1} \\ -S_G^{-1} G_{nt}^{\top} G_{nn}^{-1} & S_G^{-1} \end{pmatrix} \quad (20) $$

We infer that (dropping the $q$ argument)

$$ \begin{pmatrix} \lambda_{n,b}^{st} \\ \lambda_{t,b}^{st} \end{pmatrix} = G^{-1} \begin{pmatrix} \nabla h_{n,b}^{st,\top} \\ H_{t,b}^{st,\top} \end{pmatrix} M^{-1} F(q, \dot{q}, t) - G^{-1} \begin{pmatrix} \frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top})\dot{q} \\ \frac{d}{dt}(H_{t,b}^{st}(q)^{\top})\dot{q} \end{pmatrix} \quad (21) $$

which could be obtained using (10) (b) (d), (11) and (13). In the following, we call dynamical sticking the mode which corresponds to sticking contact points (in both tangential and normal directions), in a dynamical setting (i.e., accelerations and velocities may be non zero).

<footnote n="2">²The $q$ argument is dropped because all the terms are functions of $q$.</footnote>
**Proposition 2 (Dynamical sticking, independent constraints)** Assume that $M(q) > 0$ and $G(q) > 0$. Let $(q, \dot{q})$ and $F(q, \dot{q}, t)$ be given and consider the three dimensional case (i.e. $d = 3$). Then the all-sticking problem (14) (a) (b) has a solution if and only if:

$$ \begin{aligned} & \left\| \begin{pmatrix} [S_G^{-1} G_{nt}^\top G_{nn}^{-1}]_{2i-1, \bullet} \\ [S_G^{-1} G_{nt}^\top G_{nn}^{-1}]_{2i, \bullet} \end{pmatrix} \left( \frac{d}{dt}(\nabla h_{n,b}^{st,\top})\dot{q} - \nabla h_{n,b}^{st,\top} M^{-1} F(q, \dot{q}, t) \right) \right. \\ & \left. + \begin{pmatrix} [S_G^{-1}]_{2i-1, \bullet} \\ [S_G^{-1}]_{2i, \bullet} \end{pmatrix} \left( H_{t,b}^{st,\top} M^{-1} F(q, \dot{q}, t) - \frac{d}{dt}(H_{t,b}^{st,\top})\dot{q} \right) \right\| \\ & \le \mu_i \left| (G_{nn}^{-1} + G_{nn}^{-1} G_{nt} S_G^{-1} G_{nt}^\top G_{nn}^{-1})_{i, \bullet} \left( \nabla h_{n,b}^{st,\top} M^{-1} F(q, \dot{q}, t) - \frac{d}{dt}(\nabla h_{nb}^{st,\top})\dot{q} \right) \right. \\ & \left. - (G_{nn}^{-1} G_{nt} S_G^{-1})_{i, \bullet} \left( H_{t,b}^{st,\top} M^{-1} F(q, \dot{q}, t) + \frac{d}{dt}(H_{t,b}^{st,\top})\dot{q} \right) \right| \end{aligned} \quad (22) $$

for all $1 \le i \le m$.

**Proof 2** Directly follows from (21) and (14) (b), noting that the conditions for the invertibility of $M_b(q)$ are satisfied. ■

It may be interesting to formulate the problem another way. To that end let us define $F(q, \dot{q}, t) = C(q, \dot{q})\dot{q} + P(q) + F_{ext}(t)$, where $C(q, \dot{q})\dot{q}$ accounts for centrifugal and Coriolis forces, $P(q)$ gathers forces that derive from a potential, $F_{ext}(t)$ are external forces acting on the system, and:

$$ \begin{aligned} W_t(q) &\stackrel{\Delta}{=} S_G^{-1}[-G_{nt}^\top G_{nn}^{-1} \nabla h_{n,b}^{st,\top} M^{-1} + H_{t,b}^{st,\top} M^{-1}] \\ V_t(q, \dot{q}) &\stackrel{\Delta}{=} S_G^{-1}[G_{nt}^\top G_{nn}^{-1} \frac{d}{dt}(\nabla h_{n,b}^{st,\top})\dot{q} - \frac{d}{dt}(H_{t,b}^{st,\top})\dot{q}] + W_t(q)[C(q, \dot{q})\dot{q} + P(q)] \\ W_n(q) &\stackrel{\Delta}{=} G_{nn}^{-1}[(I + G_{nn}^{-1} G_{nt} S_G^{-1} G_{nt}^\top G_{nn}^{-1})\nabla h_{n,b}^{st,\top} M^{-1} - G_{nt} S_G^{-1} M^{-1}] \\ V_n(q, \dot{q}) &\stackrel{\Delta}{=} G_{nn}^{-1}[(I + G_{nn}^{-1} G_{nt} S_G^{-1} G_{nt}^\top G_{nn}^{-1})\frac{d}{dt}(\nabla h_{n,b}^{st,\top})\dot{q} - G_{nt} S_G^{-1} \frac{d}{dt}(H_{t,b}^{st,\top})\dot{q}] + W_n(q)[C(q, \dot{q})\dot{q} + P(q)] \end{aligned} \quad (23) $$

Then from (20) and (21) $\lambda_{n,b}^{st} = W_n(q)F_{ext}(t) + V_n(q, \dot{q})$ and $\lambda_{t,b}^{st} = W_t(q)F_{ext}(t) + V_t(q, \dot{q})$.

**Proposition 3 (Dynamical sticking, independent constraints)** Let $q, \dot{q}$ and $\mu_i$, $1 \le i \le m$, be given, and the conditions of Propositions 2 hold. Then the all-sticking problem has a solution if and only if there exists a vector $F_{ext}(t)$ such that:

$$ \begin{aligned} W_{t,i,\bullet}(q)F_{ext}(t) + V_{t,i}(q, \dot{q}) &= \text{proj}[\mathcal{D}(W_{n,i,\bullet}(q)F_{ext}(t) + V_{n,i}(q, \dot{q}), \mu_i); W_{t,i,\bullet}(q)F_{ext}(t) + V_{t,i}(q, \dot{q})] \\ &\iff 0 \in \mathcal{N}_{\mathcal{D}(W_{n,i,\bullet}(q)F_{ext}(t)+V_{n,i}(q,\dot{q}),\mu_i)}(W_{t,i,\bullet}(q)F_{ext}(t) + V_{t,i}(q, \dot{q})) \end{aligned} \quad (24) $$

for all $1 \le i \le m$.

**Proof 3** The conditions of Proposition 2 guarantee that $\lambda_{n,b}^{st} = W_n(q)F_{ext}(t) + V_n(q, \dot{q})$ and $\lambda_{t,b}^{st} = W_t(q)F_{ext}(t) + V_t(q, \dot{q})$ solve (14) (a), using (19)-(21). Using (4) we see that $F_{ext}(t)$ which solves (24) ensures that (14) (b) holds. The normal cone formulation is a generalized equation and the equivalence follows from (85). ■
Obviously (24) is equivalent to $W_{t,i,\bullet}(q)F_{\text{ext}}(t) + V_{t,i}(q, \dot{q}) \in \mathcal{D}(W_{n,i,\bullet}(q)F_{\text{ext}}(t) + V_{n,i}(q, \dot{q}), \mu_i)$, and to $||\lambda_{t,b,i}^{st}|| \le \mu_i|\lambda_{n,b,i}|$. The advantage of the formalism (24) is that efficient numerical algorithms for fixed-point calculations may be used to solve it [2, Chapters 10 and 12]. In general there is non uniqueness of solutions, because several (may be an infinity) external forces can be applied on a system to guarantee an all-sticking mode. In the following, the above developments are particularized to the static case, where both the acceleration and the velocity are zero.

**Corollary 1 (Statical all-sticking equilibrium, independent constraints)** Let $q$ and $\mu_i$, $1 \le i \le m$, be given, and the conditions of Proposition 2 hold. There exists a static equilibrium $(q, \dot{q}) = (q^*, 0)$ with all sticking contact points, if and only if there exists a vector $F_{\text{ext}}(t)$ such that:

$$ W_{t,i,\bullet}(q)(F_{\text{ext}}(t) + P(q^*)) = \text{proj}[\mathcal{D}(W_{n,i,\bullet}(q^*)(F_{\text{ext}}(t) + P(q^*)), \mu_i); W_{t,i,\bullet}(q)(F_{\text{ext}}(t) + P(q^*))] \\ \iff 0 \in \mathcal{N}_{\mathcal{D}(W_{n,i,\bullet}(q^*)(F_{\text{ext}}(t)+P(q^*)),\mu_i)}(W_{t,i,\bullet}(q)(F_{\text{ext}}(t) + P(q^*))), \quad (25) $$

for all $1 \le i \le m$.

Notice that (25) is nothing else but (4) with $v_{t,i} = 0$. Let us investigate now the system's behaviour in a right neighborhood of the time $t$ at which all-sticking is assumed to hold.

**Proposition 4** Let the conditions of Proposition 2 hold. (i) Suppose that inequality (22) is satisfied strictly for each $i \in \{1, ..., m\}$. If all the data are state and time piece-wise Lipschitz continuous, then the all-sticking mode persists on $[t, t + \varepsilon)$ for some $\varepsilon > 0$. (ii) Suppose that the inequality in (22) is replaced by an equality for some $i \in \mathcal{J} \subseteq \{1, ..., m\}$, while it is satisfied strictly at the other contact points. If $a_{t,i}(t^+) = 0$ for all $i \in \mathcal{J}$, then the all sticking mode exists at $t$. If $a_{t,i}(t^+) \ne 0$ for some $i \in \bar{\mathcal{J}} \subseteq \mathcal{J}$, then a stick $\to$ slip transition occurs at contacts $j \in \bar{\mathcal{J}}$.

**Proof 4 (i)** The dynamics of the system is given at $t$ by (10) (a) (b) (d) and (11) (13), with $\lambda_{n,u}^{st} = 0$ and $\lambda_{t,u}^{st} = 0$. Under the conditions of Proposition 2, $G(q) \succ 0$ so that this differential-algebraic equation can be reduced to an ordinary differential equation by computing the multipliers from (11) (13) and injecting their values into (10) (a). Due to the piece-wise Lipschitz continuity of its right-hand side, this ODE has a unique local solution in a right-neighborhood of $t$. This solution (position and velocity) is continuously differentiable, thus the multipliers $\lambda_{n,b}^{st}$ and $\lambda_{t,b}^{st}$ are, due to (11) (13), time-continuous functions. Consequently in a right-neighborhood of $t$ the multipliers vary continuously, therefore the contact reactions stay in the interior of Coulomb's cone for a non zero period of time, which ends the proof. (ii) The proof follows from (5).

## 2.2 The general case

It is possible to express the generalized inverse $M_b(q)^\dagger$ to obtain a general criterion, just assuming that $M(q) \succ 0$ and without restrictions on the number of contact points, hence not assuming independent constraints. First notice that after suitable rearrangements, we can rewrite equivalently (14) (a) as:

$$ \underbrace{\begin{pmatrix} M(q) & -\nabla h_{n,b}^{st}(q) & -H_{t,b}^{st}(q) \\ -\nabla h_{n,b}^{st}(q)^{\top} & 0 & 0 \\ -H_{t,b}^{st}(q)^{\top} & 0 & 0 \end{pmatrix}}_{\triangleq \bar{M}_b(q)} \begin{pmatrix} \ddot{q} \\ \lambda_{n,b}^{st} \\ \lambda_{t,b}^{st} \end{pmatrix} = \underbrace{\begin{pmatrix} -F(q, \dot{q}, t) \\ \frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top})\dot{q} \\ \frac{d}{dt}(H_{t,b}^{st}(q)^{\top})\dot{q} \end{pmatrix}}_{\triangleq \bar{F}_b(q, \dot{q}, t)} \quad (26) $$
where $\bar{M}_b(q)$ is symmetric while $M_b(q)$ in (14) is not. In this section the symmetry is needed because Proposition 14 is used, which relies on Lemma 3 which needs symmetry. Following Proposition 14 and Lemma 3, let us denote

$$ \mathbb{R}^{n \times n} \ni E(q) \stackrel{\Delta}{=} M(q) + \nabla h_{n,b}^{st}(q) \nabla h_{n,b}^{st}(q)^{\top} + H_{t,b}^{st}(q) H_{t,b}^{st}(q)^{\top} = E(q)^{\top} \succcurlyeq 0, \quad (27) $$

and $D(q) \stackrel{\Delta}{=} \begin{pmatrix} \nabla h_{n,b}^{st}(q)^{\top} E(q)^{\dagger} \nabla h_{n,b}^{st}(q) & \nabla h_{n,b}^{st}(q)^{\top} E(q)^{\dagger} H_{t,b}^{st}(q) \\ H_{t,b}^{st}(q)^{\top} E(q)^{\dagger} \nabla h_{n,b}^{st}(q) & H_{t,b}^{st}(q)^{\top} E(q)^{\dagger} H_{t,b}^{st}(q) \end{pmatrix} = \begin{pmatrix} D_{nn}(q) & D_{nt}(q) \\ D_{nt}(q)^{\top} & D_{tt}(q) \end{pmatrix} = D(q)^{\top} \in \mathbb{R}^{dm \times dm}, D_{nn}(q) \in \mathbb{R}^{m \times m}, D_{tt}(q) \in \mathbb{R}^{(d-1)m \times (d-1)m}, D_{nt}(q) \in \mathbb{R}^{m \times (d-1)m}. \text{ Notice that } M(q) \succcurlyeq 0 \Rightarrow E(q) \succcurlyeq 0, \text{ however we do not assume } M(q) \text{ full-rank in the next developments.}$$

**Lemma 1** The matrix $D(q)$ is positive semidefinite if and only if:

$$ \begin{gathered} D_{nn}(q) \succcurlyeq 0, \quad D_{nt}(q) = D_{nn}(q) D_{nn}(q)^{\dagger} D_{nt}(q), \quad D_{nt}(q)^{\top} D_{nn}(q)^{\dagger} D_{nt}(q) \preceq D_{tt}(q) \tag{28} \\ \Longleftrightarrow D_{tt}(q) \succcurlyeq 0, \quad D_{nt}(q) = D_{nt}(q) D_{tt}(q) D_{tt}(q)^{\dagger}, \quad D_{nt}(q) D_{tt}(q)^{\dagger} D_{nt}(q)^{\top} \preceq D_{nn}(q). \end{gathered} $$

Let $D(q) \succcurlyeq 0$, and define $Q(q) \stackrel{\Delta}{=} D_{tt}(q) - D_{nt}^{\top} D_{nn}(q)^{\dagger} D_{nt}(q)$. Suppose that $Q(q)$ is nonsingular, then the Moore-Penrose generalized inverse of $D(q)$ is given by:

$$ D(q)^{\dagger} = \begin{pmatrix} D_{nn}^{\dagger} + D_{nn}^{\dagger} D_{nt} Q^{-1} D_{nt}^{\top} D_{nn}^{\dagger} & -D_{nn}^{\dagger} D_{nt} Q^{-1} \\ -Q^{-1} D_{nt} D_{nn}^{\dagger} & Q^{-1} \end{pmatrix} \succcurlyeq 0. \quad (29) $$

**Proof 5** Conditions (28) are necessary and sufficient for positive semidefiniteness of $D(q)$ [12, Proposition 8.2.3]. The second part follows from the results in [59], see Lemma 3 in Appendix A, and from [12, Proposition 6.1.6 xxvii)] which states that $D(q) \succcurlyeq 0 \Leftrightarrow D(q)^{\dagger} \succcurlyeq 0$.

The next step is to calculate the generalized inverse of $\bar{M}_b(q)$, since from Proposition 14, if the system (26) is solvable then all solutions are given by $z_b = \bar{M}_b(q)^{\dagger} \bar{F}_b(q, \dot{q}, t) + (I - \bar{M}_b(q)^{\dagger} \bar{M}_b(q))y, y \in \mathbb{R}^{dm}$. In view of (29) and (79), we should be able to get expressions for $\lambda_{n,b}^{st}$ and $\lambda_{t,b}^{st}$. Using (79) we find:

$$ \bar{M}_b(q)^{\dagger} = \begin{pmatrix} E^{\dagger} - E^{\dagger}(\nabla h_{n,b}^{st} H_{t,b}^{st}) D^{\dagger} \begin{pmatrix} \nabla h_{n,b}^{st, \top} \\ H_{t,b}^{st, \top} \end{pmatrix} E^{\dagger} & -E^{\dagger}(\nabla h_{n,b}^{st} H_{t,b}^{st}) D^{\dagger} \\ -(E^{\dagger}(\nabla h_{n,b}^{st} H_{t,b}^{st}) D^{\dagger})^{\top} & DD^{\dagger} - D^{\dagger} \end{pmatrix} \quad (30) $$

from which we deduce:

$$ \begin{pmatrix} \lambda_{n,b}^{st} \\ \lambda_{t,b}^{st} \end{pmatrix} = (E^{\dagger}(\nabla h_{n,b}^{st} H_{t,b}^{st}) D^{\dagger})^{\top} F(q, \dot{q}, t) + (DD^{\dagger} - D^{\dagger}) \begin{pmatrix} \frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top}) \dot{q} \\ \frac{d}{dt}(H_{t,b}^{st}(q)^{\top}) \dot{q} \end{pmatrix} \quad (31) $$
In view of Lemma 3 and Proposition 14, this solution is a minimum norm solution. After calculations using Lemma 1 we find that:

$$ \lambda_{n,b}^{st} = \left[ \begin{aligned} & [(D_{nn}(q) - I)D_{nn}(q)^{\dagger}(I + D_{nt}(q)Q(q)^{\dagger}D_{nt}(q)^{\top}D_{nn}(q)^{\dagger}) - D_{nt}(q)Q(q)^{\dagger}D_{nt}(q)D_{nn}(q)^{\dagger}] \frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top})\dot{q} \\ & - [(D_{nn}(q) - I)D_{nn}(q)^{\dagger}D_{nt}(q) + D_{nt}(q)] Q(q)^{\dagger} \frac{d}{dt}(H_{t,b}^{st}(q)^{\top})\dot{q} \\ & - D_{nn}(q)^{\dagger}[-(I + D_{nt}(q)Q(q)^{\dagger}D_{nt}(q)^{\top}D_{nn}(q)^{\dagger})\nabla h_{n,b}^{st}(q)^{\top} + D_{nt}(q)Q(q)^{\dagger}H_{t,b}^{st}(q)^{\top}] E(q)^{\dagger}F(q, \dot{q}, t) \end{aligned} \right] $$

$$ \lambda_{t,b}^{st} = \left[ \begin{aligned} & [D_{nt}(q)^{\top}D_{nn}(q)^{\dagger}(I + D_{nt}(q)Q(q)^{\dagger}D_{nt}(q)^{\top}D_{nn}(q)^{\dagger}) - (D_{tt}(q) - I)Q(q)^{\dagger}D_{nt}(q)D_{nn}(q)^{\dagger}] \frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top})\dot{q} \\ & + [-D_{nt}(q)^{\top}D_{nn}(q)^{\dagger}D_{nt}(q)Q(q)^{\dagger} + (D_{tt}(q) - I)Q(q)^{\dagger}] \frac{d}{dt}(H_{t,b}^{st}(q)^{\top})\dot{q} \\ & + \dot{Q}(q)^{\dagger}[D_{nt}(q)D_{nn}(q)^{\dagger}\nabla h_{n,b}^{st}(q)^{\top} - H_{t,b}^{st}(q)^{\top}] E(q)^{\dagger}F(q, \dot{q}, t). \end{aligned} \right] \quad (32) $$

Let us denote $\lambda_{n,b}^{st} = W_n(q)F_{ext}(t) + V_n(q, \dot{q})$ and $\lambda_{t,b}^{st} = W_t(q)F_{ext}(t) + V_t(q, \dot{q})$. The matrices $W_n(q)$, $V_n(q, \dot{q})$, $W_t(q)$ and $V_t(q, \dot{q})$ can be obtained from (29) and (31), and $F(q, \dot{q}, t) = C(q, \dot{q})\dot{q} + P(q) + F_{ext}(t)$ (Coriolis and centrifugal torques, forces that derive from a potential, and external forces and torques).

**Proposition 5 (Dynamical sticking, arbitrary constraints)** Let $M(q) \ge 0$, $m(= m_b)$ and $n$ be arbitrary positive integers, $q$ and $\dot{q}$ be given, and the conditions of Lemma 1 hold. The all-sticking problem (14) (a) (b) has a solution if and only if there exists a vector $F_{ext}(t)$ such that:

$$ \begin{aligned} W_{t,i,\bullet}(q)F_{ext}(t) + V_{t,i}(q, \dot{q}) &= \operatorname{proj}[\mathcal{D}(W_{n,i,\bullet}(q)F_{ext}(t) + V_{n,i}(q, \dot{q}), \mu_i); W_{t,i,\bullet}(q)F_{ext}(t) + V_{t,i}(q, \dot{q})] \\ &\iff 0 \in \mathcal{N}_{\mathcal{D}(W_{n,i,\bullet}(q)F_{ext}(t)+V_{n,i}(q,\dot{q}),\mu_i)}(W_{t,i,\bullet}(q)F_{ext}(t) + V_{t,i}(q, \dot{q})), \end{aligned} \quad (33) $$

for all $1 \le i \le m$, where the matrices $W_n(q)$, $W_t(q)$ and the vectors $V_n(q, \dot{q})$, $V_t(q, \dot{q})$ are obtained from (31).

The proof is based on the use of (4). The static equilibrium conditions are easily deduced setting $\dot{q} = 0$ in (31) and (33). KKT problems as (14) (a) have been studied using Moore-Penrose generalized inverses in [75] with full-rank mass matrix, where problems related to the dependence of numerical solutions to chosen units are analysed. Many other results for (26) may be found in [76, 55, 73, 27, 74], with both uniqueness of contact forces analysis and algorithms for solving the KKT problem.

## 3 Unilaterally constrained systems

Now we assume that $m_b = 0$ and $\mathcal{I}_u = \mathcal{I}_u^{st}$, thus $m_u = m_u^{st} = m$. We first analyse the contact problem disregarding the constraints imposed on the contact force by friction (hence implicitly assuming an infinite coefficient of friction), thus stating necessary conditions for the whole problem to be solvable. Then the Coulomb's cone constraint is added to the problem.
## 3.1 The contact MLCP

We assume that the $m$ unilateral contacts are active, thus $h_{n,u}(q(t)) = 0$ and $\nabla h_{n,u}(q(t))^{\top}\dot{q}(t) = 0$, with continuous velocity. If we assume right-continuous acceleration ($\ddot{q}(t) = \ddot{q}(t^+) = \lim_{s\searrow t} \ddot{q}(s)$ and $a_{n,u}(t) = \lim_{s\searrow t} a_{n,u}(s)$), and if $F(q(t), \dot{q}(t), t)$ is such that $a_{n,u}(t) = 0$, the contact remains active in a right-neighborhood of $t$. However in the general case, it could be that some entries $a_{n,u,i}$ satisfy $a_{n,u,i}(t^+) > 0$: some contacts are deactivated in a right-neighborhood of $t$ while they are sticking at $t^3$. Some continuity conditions on the data will have to be imposed to guarantee that sticking persists. The detection of tangentially sticking contacts which undergo a stick $\to$ slip transition with $a_{t,u}(t^-) = 0$ but $a_{t,u}(t^+) \neq 0$, can be done using the Coulomb's law at the acceleration level. When $M(q) > 0$, we can use the system in (10), which in this case reduces to:

$$ \left\{ \begin{array}{l} (a) \ M(q)\ddot{q} + F(q, \dot{q}, t) = \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} + H_{t,u}^{st}(q)\lambda_{t,u}^{st} \\ (c) \ H_{t,u}^{st}(q)^{\top}M(q)^{-1}H_{t,u}^{st}(q)\lambda_{t,u}^{st} + R_{t,u}^{st}(q, \dot{q}, t, \lambda_{n,u}^{st}) = 0 \\ (e) \ 0 \le \lambda_{n,u}^{st} \perp a_{n,u} = \nabla h_{n,u}^{st}(q)^{\top}\ddot{q} + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \ge 0, \end{array} \right. \quad (34) $$

with $R_{t,u}^{st}(q, \dot{q}, t, \lambda_{n,u}^{st}) = H_{t,u}^{st}(q)^{\top}M(q)^{-1}[-F(q, \dot{q}, t) + \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st}] + \frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q}$. If we deal with 3D friction, $d=3$ and $\lambda_{t,u}^{st}$ has dimension $2m_u = 2m$. In case of planar friction $d=2$ and it has dimension $m_u = m$. The problem in (34) is an MLCP with unknowns $\lambda_{n,u}^{st}$ and $\lambda_{t,u}^{st}$. There are various ways to analyse MLCPs [26]. In its full generality we should assume that $M(q) = M(q)^{\top} \ge 0$ and $H_{t,u}^{st}(q)^{\top}M(q)^{-1}H_{t,u}^{st}(q) \ge 0$ only. The following assumption allows us to transform easily the MLCP into an LCP:

**Assumption 1** Let $M(q) > 0$ and the matrix $\mathbb{R}^{(d-1)m \times (d-1)m} \ni H_{t,u}^{st}(q)^{\top}M(q)^{-1}H_{t,u}^{st}(q) > 0$ ($\Leftrightarrow H_{t,u}^{st}(q)$ has full column rank).

Assumption 1 allows us to eliminate $\lambda_{t,u}^{st}$ and to obtain a constrained Delassus' matrix. Following the same lines of calculation as in [17], one obtains that (34) (e) is rewritten as an augmented contact LCP whose matrix is:

$$ A_c^{st}(q) = \nabla h_{n,u}^{st}(q)^{\top} M_c^{st}(q)^{-1} \nabla h_{n,u}^{st}(q) \quad (35) $$

with

$$ M_c^{st}(q)^{-1} \triangleq M(q)^{-1} - M(q)^{-1} H_{t,u}^{st}(q) (H_{t,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q))^{-1} H_{t,u}^{st}(q)^{\top} M(q)^{-1}. \quad (36) $$

The matrix $M_c^{st}(q)^{-1}$ is symmetric positive semi-definite with rank $n-(d-1)m$ (using [17, Lemma 2], and the fact that Assumption 1 implies that $(d-1)m \le n$), and it is denoted this way to mimic the case without bilateral constraints. However $A_c^{st}(q)$ may have full rank, being the Schur complement of $H_{t,u}^{st}(q)^{\top}M(q)^{-1}H_{t,u}^{st}(q)$ in the matrix $\begin{pmatrix} \nabla h_{n,u}^{st}(q)^{\top} \\ H_{t,u}^{st}(q)^{\top} \end{pmatrix} M(q)^{-1}(\nabla h_{n,u}^{st}(q) H_{t,u}^{st}(q))$ (a necessary and sufficient condition for this is that all vectors $\nabla h_{n,u,i}^{st}(q)$ and $H_{t,u,i}^{st}(q)$ be independent [17] [19, p.252]). The contact LCP obtained from (34) is therefore equal to [17, section 4]

$$ 0 \le \lambda_{n,u}^{st} \perp A_c^{st}(q)\lambda_{n,u}^{st} + H_c^{st}(q, \dot{q}, t) \ge 0, \quad (37) $$

<footnote n="3">³This kind of transition is often assumed to exist in biped robots control, for instance, when one foot detaches from the ground at the end of a step where both feet are in persistent all-sticking mode.</footnote>
where:

$$ H_c^{st}(q, \dot{q}, t) = -\nabla h_{n,u}^{st}(q)^{\top} M_c^{st}(q)^{-1} F(q, \dot{q}, t) \\ -\nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q) (H_{t,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q))^{-1} \frac{d}{dt} (H_{t,u}^{st}(q)^{\top}) \dot{q} + \frac{d}{dt} (\nabla h_{n,u}^{st}(q)^{\top}) \dot{q}. \quad (38) $$

The analytical computation of such expressions, is in general quite cumbersome. We have the following result that is an extension of Proposition 1, and which holds because $A_c^{st}(q)$ is symmetric positive (semi) definite. It characterizes the augmented contact LCP (35) (36) (37) (38).

**Proposition 6 (Sticking unilateral constraints) Let Assumption 1 hold. Then:**

* (i) The contact LCP in (37) has a unique solution for any $H_c^{st}(q, \dot{q}, t)$ if and only if $A_c^{st}(q) > 0$.

* (ii) Let $\lambda_{n,u,1}^{st}$ and $\lambda_{n,u,2}^{st}$ be two solutions of the contact LCP in (37). Then $\nabla h_{n,u}^{st}(q)(\lambda_{n,u,1}^{st} - \lambda_{n,u,2}^{st}) \in \ker(M_c^{st}(q)^{-1})$, and $(\lambda_{n,u,1}^{st,T} - \lambda_{n,u,2}^{st,T})H_c^{st}(q, \dot{q}) = 0$.

* (iii) Let $\mathrm{Im}(\nabla h_{n,u}^{st}(q)) \cap \mathrm{Im}(H_{t,u}^{st}(q)) = \{0\}$ for all $q \in \mathbb{R}^n$. Let $\lambda_{n,u,1}^{st}$ and $\lambda_{n,u,2}^{st}$ be two solutions of the LCP in (37). Then $\nabla h_{n,u}^{st}(q)(\lambda_{n,u,1}^{st} - \lambda_{n,u,2}^{st}) = 0$. If the LCP in (37) is solvable then the acceleration $\ddot{q}$ exists and is unique.

* (iv) The LCP in (37) is solvable if for any $z \in \mathbb{R}^n$ in the set of solutions of the homogeneous LCP $0 \le z \perp A_c^{st}(q)z \ge 0$ one has $z^\top H_c^{st}(q, \dot{q}, t) \ge 0$; if $\frac{d}{dt}(\nabla H_{t,u}^{st}(q)) = 0$ and $\frac{d}{dt}(\nabla h_{n,u}^{st}(q)) = 0$ and $\mathrm{Im}(\nabla h_{n,u}^{st}(q)) \cap \mathrm{Im}(H_{t,u}^{st}(q)) = \{0\}$ for all $q \in \mathbb{R}^n$, the LCP in (37) is solvable.

* (v) Suppose that $-\lambda_{n,u}^{\top}(\nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q)^{\top}) \dot{q} + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top}) \dot{q}) \ge 0$, then the LCP in (37) is solvable.

**Proof 6** The proof follows from [15, Proposition 7], just noting that due to the all-sticking assumption, the system (34) (or the system (9) (a) (c) (e) (12)) is identical to the system [15, Equation (16)], replacing the bilateral constraints of [15, Equation (16)] by the frictional sticking constraints.

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

Let us now analyse various cases of the so-called planar rocking block. Let $(x, y)$ be the coordinates of its gravity center (supposed to be the same as the geometric center), $\theta$ is the orientation angle (see Figure 2), $L$ is the block's width, $l$ is its height. This system has two unilateral constraints at the potential contact/impact points $A_1$ and $A_2$ (provided the base line is assumed to be concave) when $y \le \frac{\sqrt{l^2+L^2}}{2}$: $h_{n,u,1}(q) = y - \frac{l}{2}\cos(\theta) + \frac{L}{2}\sin(\theta) \ge 0$, and $h_{n,u,2}(q) = y - \frac{l}{2}\cos(\theta) - \frac{L}{2}\sin(\theta) \ge 0$. An external force $F_{ext}(\cdot)$ acts at the gravity center. The dynamics of the block with Coulomb's friction is given by [19]:
$$
\left\{
\begin{aligned}
m\ddot{x}(t) &= \lambda_{t,u,1}(t) + \lambda_{t,u,2}(t) - F_{ext}^x \\
m\ddot{y}(t) &= \lambda_{n,u,1}(t) + \lambda_{n,u,2}(t) - F_{ext}^y \\
I_G \ddot{\theta}(t) &= \lambda_{n,u,1}(t) \left(\frac{l}{2}\sin(\theta(t)) + \frac{L}{2}\cos(\theta(t))\right) + \lambda_{n,u,2}(t) \left(\frac{l}{2}\sin(\theta(t)) - \frac{L}{2}\cos(\theta(t))\right) \\
&\quad + \left(\frac{l}{2}\cos(\theta(t)) - \frac{L}{2}\sin(\theta(t))\right) \lambda_{t,u,1} + \left(\frac{l}{2}\cos(\theta(t)) + \frac{L}{2}\sin(\theta(t))\right) \lambda_{t,u,2} \\
&0 \le \lambda_{n,u}(t) \perp h_{n,u}(q(t)) \ge 0 \\
\lambda_{t,u,i}(t) &\in -\mu_i \lambda_{n,u,i}(t) \operatorname{sgn}(v_{t,i}(t)), \quad i = 1, 2,
\end{aligned}
\right.
\tag{39}
$$

where $n = 3$, $m = 2$, $d = 2$, we have used the expression of Coulomb's law in 2D (see section 1.2), $I_G = m\frac{l^2+L^2}{12}$, $\mu_i > 0$ is the friction coefficient at contact $i$, and $v_{t,i}$ is the tangential velocity at the point $i$, i.e., $v_{t,1} = \dot{x} + (\frac{l}{2}\cos(\theta) - \frac{L}{2}\sin(\theta))\dot{\theta}$ at $A_1$ and $v_{t,2} = \dot{x} + (\frac{l}{2}\cos(\theta) + \frac{L}{2}\sin(\theta))\dot{\theta}$ at $A_2$ (from which $v_{t,1} = v_{t,2}$ when $\theta = 0$). With $q = (x, y, \theta)^\top$, one can identify $M(q) = M$, $F_{ext}$, $\nabla h_{n,u}(q) = \begin{pmatrix} 0 & \frac{l}{2} & \sin(\theta) + \frac{L}{2}\cos(\theta) \\ 0 & 1 & \frac{l}{2}\sin(\theta) - \frac{L}{2}\cos(\theta) \end{pmatrix}^\top$ and $H_{t,u}(q) = \begin{pmatrix} 1 & 0 & \frac{l}{2}\cos(\theta) - \frac{L}{2}\sin(\theta) \\ 1 & 0 & \frac{l}{2}\cos(\theta) + \frac{L}{2}\sin(\theta) \end{pmatrix}^\top$, in (1) from (39). The contact LCP is given by $0 \le \lambda_{n,u} \perp D_{nn}(\theta)\lambda_{n,u} + \nabla f(\theta)^\top M^{-1}H_t(\theta)\lambda_{t,u} + B(\theta, \dot{\theta}) \ge 0$, with

$$D_{nn}(\theta) = \begin{pmatrix} \frac{1}{m} + \frac{1}{4I_G}(l \sin(\theta) + L \cos(\theta))^2 & \frac{1}{m} + \frac{1}{4I_G}(l^2 \sin^2(\theta) - L^2 \cos^2(\theta)) \\ \frac{1}{m} + \frac{1}{4I_G}(l^2 \sin^2(\theta) - L^2 \cos^2(\theta)) & \frac{1}{m} + \frac{1}{4I_G}(l \sin(\theta) - L \cos(\theta))^2 \end{pmatrix},$$

and

$$B(\theta, \dot{\theta}) = \begin{pmatrix} -g + \frac{1}{2}\dot{\theta}^2(l \cos(\theta) - L \sin(\theta)) \\ -g + \frac{1}{2}\dot{\theta}^2(l \cos(\theta) + L \sin(\theta)) \end{pmatrix}.$$

**Example 1** Let us consider the system in Figure 1, where both contacts are supposed to stick in both normal and tangential directions (hence $v_{t,1} = v_{t,2} = 0$, $\dot{\theta} = 0$, $\theta = 0$, $\dot{x} = 0$, $\dot{y} = 0$). The configuration is $q^{st} = (x, \frac{l}{2}, 0)^\top$. Here we have $\nabla h_{n,u}^{st}(q^{st}) = \begin{pmatrix} 0 & 1 & \frac{L}{2} \\ 0 & 1 & -\frac{L}{2} \end{pmatrix}^\top \in \mathbb{R}^{3 \times 2}$ which is full column rank, but $H_{t,u}^{st}(q^{st}) = \begin{pmatrix} 1 & 0 & \frac{l}{2} \\ 1 & 0 & \frac{l}{2} \end{pmatrix}^\top \in \mathbb{R}^{3 \times 2}$, hence Assumption 1 is not satisfied and one cannot construct the contact LCP in (35) (36) (37) (38). We thus postpone the analysis later, see Example 4.

**Example 2** Let us consider the system in Figure 2. Only point $A_2$ is in contact and sticks in both directions ($v_{t,2} = 0$, $h_{n,u,2}(q) = y - \frac{l}{2}\cos(\theta) - \frac{L}{2}\sin(\theta) = 0$). Since there is only one contact, $H_{t,u}(q) = H_{t,u,2}(q) = (1, 0, \frac{l}{2}\cos(\theta) + \frac{L}{2}\sin(\theta))^\top$ is a vector and Assumption 1 is satisfied since this vector is not zero. We have $\nabla h_{n,u}(q) = \nabla h_{n,u,2}(q) = (0, 1, \frac{l}{2}\sin(\theta) - \frac{L}{2}\cos(\theta))^\top$. The matrix

$$A_c^{st}(q) = \frac{12}{m(l^2 + L^2)} \left( \frac{l^2 + L^2}{12} + \left( \frac{l}{2} \sin(\theta) - \frac{L}{2} \cos(\theta) \right)^2 - \frac{\left( \frac{l}{2} \sin(\theta) - \frac{L}{2} \cos(\theta) \right)^2 \left( \frac{l}{2} \cos(\theta) + \frac{L}{2} \sin(\theta) \right)^2}{\frac{l^2+L^2}{12} + \left( \frac{l}{2} \cos(\theta) + \frac{L}{2} \sin(\theta) \right)^2} \right) \quad (40)$$

is a positive scalar, since the vectors $H_{t,u,2}(q)$ and $\nabla h_{n,u,2}(q)$ are independent. Therefore the contact LCP in (35) (36) (37) (38) always has a unique solution. Let us remind that obviously, this does not
<figure><img src="image_4.png" /><figcaption>Figure 1: A 2D rigid block in static equilibrium on a rough ground.</figcaption></figure>

mean that the contact is sticking for any $F_{ext}$, as we are treating only a necessary condition. The contact LCP being scalar, it is however easy to envision its two possible outcomes. The first term in the right-hand side of (38) is:

$$ -\nabla h_{n,u}^{st}(q)^{\top} M_c^{st}(q)^{-1} F(q, \dot{q}, t) = \frac{1}{m} F_{ext}^{y} + \frac{12}{l^2+L^2} \left( \left( \frac{l}{2} \sin(\theta) - \frac{L}{2} \cos(\theta) \right) \left( \frac{l}{2} \sin(\theta) + \frac{L}{2} \cos(\theta) \right) \right) \times \left( 1 + 12 \frac{\left( \frac{l}{2} \cos(\theta) + \frac{L}{2} \sin(\theta) \right)^2}{l^2 + L^2} \right)^{-1} F_{ext}^{x} \quad (41) $$

It is inferred from (41), (35) and (38), that whatever the values of $F_{ext}^x$ and of $\dot{q}$, a large enough negative $F_{ext}^y$ exists to guarantee that $\lambda_{n,u}^{st} > 0$. Notice also that if $\theta = \arctan(\frac{L}{l})$ (the block's diagonal $A_2A_3$ is normal to the constraint), then $F_{ext}^x$ plays no role. However, it will play a role when Coulomb's cone is taken into account (for the moment, the analysis implicitly assumes an infinite coefficient of friction).

<figure><img src="image_5.png" /><figcaption>Figure 2: A 2D rigid block rotating around $A_1$.</figcaption></figure>
**Example 3** Let us consider the system in Figure 3. It has four degrees of freedom, with $q = (x, y, \theta, \alpha)^T$. Let us denote the length of the pendulum, ($x_p, y_p$) the coordinates of its gravity center, $x_p = x + \frac{a}{2} \cos(\alpha + \theta) - \frac{1}{2} \sin(\theta)$, $y_p = y + \frac{1}{2} \cos(\alpha + \theta) + \frac{a}{2} \sin(\alpha + \theta)$ (in the figure the system is depicted with $\theta = 0$). The dynamics of this system is detailed in section D. The big discrepancy with the above two systems, is that the inertia matrix is no longer constant diagonal. Even if $\theta = 0$ rad, and $\dot{\theta} = 0$ rad/s, the pendulum's motion plays a role in the sticking conditions. The inertial terms in (90) become when $\dot{\theta} = 0$ rad/s, $\theta = 0$ rad, $\dot{x} = 0$ m/s, $\dot{y} = 0$ m/s, and $y = \frac{1}{2}$ m, equal to:

$$ 
\begin{aligned}
F^x(q, \dot{q}, t) &= -\frac{\bar{m}a}{2} \dot{\alpha}^2 \cos(\alpha) + F_{ext}^x \\
F^y(q, \dot{q}, t) &= -\frac{\bar{m}a}{2} \dot{\alpha}^2 \sin(\alpha) + F_{ext}^y \\
F^\theta(q, \dot{q}, t) &= -\frac{\bar{m}al}{4} \dot{\alpha}^2 \sin(\alpha) \\
F^\alpha(q, \dot{q}, t) &= 0.
\end{aligned}
\qquad (42) $$

From (42) the first term in (38), it appears that not only the external force but also the inertial effects due to the pendulum play a role in the sticking problem. We do not continue the calculations, because inverting the mass matrix in (88) yields quite cumbersome calculations, which should be performed numerically.

<figure><img src="image_6.png" /><figcaption>Figure 3: A 2D rigid block in static equilibrium on a rough ground with an inverted pendulum at its top.</figcaption></figure>

**Remark 2** The LCP matrix $A_c^{st}(q)$ is supposed to be full rank only in (i) of Proposition 6. The other items hold for positive semi definite LCP matrices, that is, redundant constraints are allowed. It is noteworthy that sticking friction under Assumption 1 does not destroy the convexity of the contact complementarity problem, since the contact LCP matrix $A_c^{st}(q)$ is always symmetric and positive semi
definite, and may even be positive definite if some conditions are fullfilled as alluded to above. This is a major discrepancy between sticking and sliding friction, and might explain why sticking modes have always been considered as being more tractable than sliding modes in the literature.

The generalized force $F(q, \dot{q}, t)$ plays no role in item (v) of Proposition 6, because of the specific structure of $H_c^{st}(q, \dot{q}, t)$ in (38). We note anyway that item (v) of Proposition 6 is less general than [70, Theorem 3.1 (B)], or [54, Theorem 5], which both rely on a similar condition stemming from [25, Theorem 3.8.6] (see Theorem 1 in Appendix), but incorporate the Coulomb's cone constraint (a similar comment applies to Proposition 7). Finally we can start from (9) (a) (c) (e) and adapt [15, Proposition 8] straightforwardly to the unilateral all-sticking contact case. The interest is that it relaxes Assumption 1, in the sense that $H_{t,u}^{st}(q)$ needs not be full column rank.

**Proposition 7** Let us consider the contact MLCP with unknowns $\ddot{q}$ and $\lambda_{n,u}^{st}$ in (9) (a) (c) (e). Assume that $M(q) \succ 0$. Suppose that $\frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q} \in \text{im}(H_{t,u}^{st}(q)^{\top})$. Then the solvability of the contact MLCP holds if and only if $\nabla h_{n,u}^{st}(q)[z + y] + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \ge 0$ for some $z \in \mathbb{R}^n$ and $y \in \ker(H_{t,u}^{st}(q)^{\top})$.

As alluded to above, this result is also close to [54, Theorem 5], nevertheless the obtained conclusions differ from those in [54, Theorem 5], because we do not introduce Coulomb's friction at the acceleration level as in (5). Indeed this modeling assumption allows one to transform the persistent-sticking mode MLCP in (34) into an LCP, introducing suitable slack variables. The fundamental assumption (inequality (35) in [54]) is in our notations: $(\lambda_{n,u}^{st,\top}, \lambda_{t,u}^{st,\top}) \frac{d}{dt} \begin{pmatrix} \nabla h_{n,u}^{st}(q)^{\top} \\ H_{t,u}^{st}(q)^{\top} \end{pmatrix} \dot{q} \ge 0$ for all $\lambda_{t,u}^{st} \in \mathcal{D}(\mu, \lambda_{n,u}^{st})$ and $(\lambda_{n,u}^{st}, \lambda_{t,u}^{st}) \in \ker \begin{pmatrix} \nabla h_{n,u}^{st} \\ H_{t,u}^{st} \end{pmatrix}$.

## 3.2 Analysis of acceleration/multipliers inclusions

In the foregoing section, the contact LCP is obtained using Assumption 1. If the full-rank conditions are relaxed, one obtains a contact MLCP that is harder to analyse. Inclusions for the acceleration may be obtained, starting from (9) and relaxing Assumption 1 (hence allowing for an arbitrary number of constraints). Starting from (9) (a) (c) (e) one obtains the MLCP with unknowns $\ddot{q}$, $\lambda_{t,u}^{st}$ and $\lambda_{n,u}^{st}$:

$$ \underbrace{\begin{pmatrix} M(q) & -H_{t,u}^{st}(q) \\ H_{t,u}^{st}(q)^{\top} & 0 \end{pmatrix}}_{\triangleq M_u^{st}(q)} \underbrace{\begin{pmatrix} \ddot{q}(t) \\ \lambda_{t,u}^{st}(t) \end{pmatrix}}_{\triangleq z(t)} + \underbrace{\begin{pmatrix} F(q, \dot{q}, t) \\ \frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q} \end{pmatrix}}_{\triangleq F_u^{st}(q, \dot{q}, t)} = \begin{pmatrix} \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st}(t) \\ 0 \end{pmatrix} \quad (43) $$
$$ 0 \le \lambda_{n,u}^{st}(t) \perp \nabla h_{n,u}^{st}(q)^{\top}\ddot{q} + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \ge 0. $$

Similarly to (14), one assumes that $M(q) \succ 0$ in (43). One sees that $M_u^{st}(q)$ is positive semidefinite, non symmetric, but may have full rank (see [15, Equations (67) and (69)] and [12, Lemma 6.4.20] for
rank calculation). Let us transform the MLCP (43) into a more convenient form. First let us define the set of $\mathbb{R}^n$:

$$ K \triangleq \{w \in \mathbb{R}^n | \nabla h_{n,u}^{st}(q)^\top w + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^\top) \dot{q} \ge 0\}. \quad (44) $$

Obviously for given $q$ and $\dot{q}$, $K$ is a closed polyhedral convex set. We assume that in the sequel of this section, the following holds:

**Assumption 2** Given $q$ and $\dot{q}$, the set $K$ is non empty.

This is a kind of constraint qualification. The complementarity conditions in (43) are equivalent to the inclusion: $\lambda_{n,u}^{st} \in -\mathcal{N}_{\mathbb{R}_+^m}(\nabla h_{n,u}^{st}(q)^\top \ddot{q} + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^\top) \dot{q})$ (see e.g. [19, section B.2]). Using the fact that $\partial \psi_{\mathbb{R}_+^m}(x) = \mathcal{N}_{\mathbb{R}_+^m}(x)$, and the chain rule of Convex Analysis [58, Theorem 23.9], the term $\nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st}$ can thus be rewritten equivalently as $-\mathcal{N}_K(\ddot{q})$, provided that $\mathrm{im}(\nabla h_{n,u}^{st}(q)^\top)$ contains a point of $\mathbb{R}_+^m - \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^\top)\dot{q})$ (which is guaranteed by Assumption 2)$^4$. We can rewrite equivalently the MLCP in (43) as the inclusion $M_u^{st}(q)z(t) + F_u^{st}(q, \dot{q}, t) \in -\begin{pmatrix} \mathcal{N}_K(\ddot{q}) \\ 0 \end{pmatrix}$, where the normal multiplier $\lambda_{n,u}^{st}$ no longer appears explicitly but is replaced by a set-valued function of $\ddot{q}$. Consequently we have transformed the all-sticking contact MLCP in (43) into a problem with unknowns $\ddot{q}$ and $\lambda_{t,u}^{st}$ only. Now we have $\mathcal{N}_{\mathbb{R}^{(d-1)m}}(\lambda_{t,u}^{st}) = \{0\}$. Let us define the closed convex set $\tilde{K} = K \times \mathbb{R}^{(d-1)m} \ni (\ddot{q}, \lambda_{t,u}^{st})$, we thus obtain the inclusion equivalent to the MLCP in (43):

$$ M_u^{st}(q)z(t) + F_u^{st}(q, \dot{q}, t) \in -\mathcal{N}_{\tilde{K}}(z(t)) = -\partial \psi_{\tilde{K}}(z(t)), \quad (45) $$

where we used the fact that given two closed convex sets $K_1 \subseteq \mathbb{R}^{n_1}$ and $K_2 \subseteq \mathbb{R}^{n_2}$, and two vectors $x_i \in K_i$, $i = 1, 2$, one has $\psi_{K_1 \times K_2}(x_1, x_2) = \psi_{K_1}(x_1) + \psi_{K_2}(x_2)$, so that $\partial \psi_{K_1 \times K_2}(x_1, x_2) = \begin{pmatrix} \partial \psi_{K_1}(x_1) \\ \partial \psi_{K_2}(x_2) \end{pmatrix}$. The inclusion in (45) is in turn equivalent to the variational inequality: Find $z(t) \in \tilde{K}$ such that

$$ \langle M_u^{st}(q)z(t) + F_u^{st}(q, \dot{q}, t), y - z(t) \rangle \ge 0, \text{ for all } y \in \tilde{K}. \quad (46) $$

The existence and uniqueness of solutions for this type of problems, has been studied in [3] (the existence of solutions to (46) could also be analysed using results in [50], due to the particular form of the function $\varphi(\cdot)$ in (81)), whose results have been used in [20] for the analysis of an inclusion simpler than (45) (equivalently an MLCP simpler than (43), see [20, Equation (16)]; and the MLCP in [20, Equation (34)] has a structure different than (43) because it involves two different multipliers). To analyse the contact MLCP (43), we will rely on the material in Appendix B. It is noteworthy that the problem we are dealing with is hard enough because $M_u^{st}(q)$ is non symmetric and the indicator function $\psi_{\tilde{K}}(\cdot)$ is not strictly convex. Clearly (86) cannot be used. Let us define $S_{n,u}^{st} = \{q \in \mathbb{R}^n | h_{n,u}^{st}(q) \ge 0\}$, which under some suitable conditions on the gap functions $h_{n,u,i}^{st}: \mathbb{R}^n \to \mathbb{R}$ is a prox-regular set [4]. For a matrix $H \in \mathbb{R}^{n \times m}$ and a set $C \subseteq \mathbb{R}^n$, $H^{-1}(C) = \{x \in \mathbb{R}^m | Hx \in C\}$. The MFCQ definition is given in Appendix C, as well as the definition of normal cones to prox-regular sets.

<footnote>$^4$See Appendix B and C for the definition of the indicator function, its subdifferential and the normal cone to a set $K$.</footnote>
**Proposition 8 (Sticking unilateral constraints)** Let us consider the variational inequality in (46). Suppose that the constraints $h_{n,u,i}^{st}(q)$, $1 \le i \le m$, satisfy the MFCQ. Let

$$ \mathcal{E}_{n,u}^{st}(q) \triangleq \{(x, y) \in \mathbb{R}^n \times \mathbb{R}^{(d-1)m} | x \in \mathcal{T}_{S_{n,u}^{st}}(q) \cap \ker(M(q)) \cap \ker(H_{t,u}^{st}(q)^{\top}), y \in (H_{t,u}^{st}(q))^{-1}(\mathcal{N}_{S_{n,u}^{st}}(q))\} \quad (47) $$

Then:

• (i) If $\mathcal{E}_{n,u}^{st}(q) = \{(0_n, 0_{(d-1)m})\}$, the contact MLCP in (43), is solvable,

• (ii) if $\mathcal{E}_{n,u}^{st}(q) \neq \{(0_n, 0_{(d-1)m})\}$ and there exists $x_0 \in K$ such that $F_u^{st}(q, \dot{q}, t)^{\top}v + x_0^{\top}H_{t,u}^{st}(q)v_m > 0$ for all $v = (v_n^{\top}, v_m^{\top})^{\top} \neq (0^n, 0^{(d-1)m})^{\top}$, $v_n \in \mathbb{R}^n$, $v_m \in \mathbb{R}^{(d-1)m}$, and $v_m \in (H_{t,u}^{st})^{-1}(\mathcal{N}_{S_{n,u}^{st}}(q))$, then the contact MLCP is solvable.

• (iii) If $(\ddot{q}^1, \lambda_{t,u}^{st,1})$ and $(\ddot{q}^2, \lambda_{t,u}^{st,2})$ are two solutions of the contact MLCP (43), then $\ddot{q}^1 - \ddot{q}^2 \in \ker(M(q))$.

**Proof 7** The proof relies on Proposition 16 in Appendix B. First we see that (46) is (81) with $\varphi(\cdot) = \psi_{\tilde{K}}(\cdot)$, $\mathbf{M} = M_u^{st}(q)$, $\mathbf{q} = F_u^{st}(q, \dot{q}, t)$, $u = z(t)$, $v = y$. The first step consists in the calculation of the set $(\mathrm{dom}(\varphi))_{\infty} \cap \ker\{\mathbf{M} + \mathbf{M}^{\top}\} \cap \mathcal{K}(\mathbf{M}, \varphi)$ in Proposition 16 (a). We have $(\mathrm{dom}(\varphi))_{\infty} = \tilde{K}_{\infty} = \{w \in \mathbb{R}^n | \nabla h_{n,u}^{st}(q)^{\top}w \ge 0\} \times \mathbb{R}^{(d-1)m}$, which is under the MFCQ equal to $\mathcal{T}_{S_{n,u}^{st}}(q) \times \mathbb{R}^{(d-1)m}$, and we used Proposition 15 e). Secondly $\ker(M_u^{st}(q) + M_u^{st}(q)^{\top}) = \ker\left(\begin{pmatrix} 2M(q) & 0 \\ 0 & 0 \end{pmatrix}\right) = \{(x, y) \in \mathbb{R}^{n \times (d-1)m} | x \in \ker(M(q)), y \in \mathbb{R}^{(d-1)m}\}$. The third set is $\mathcal{K}(M_u^{st}(q), \psi_{\tilde{K}}) = \{(x, y) \in \mathbb{R}^{n \times (d-1)m} | M_u^{st}(q)\begin{pmatrix} x \\ y \end{pmatrix} = \tilde{K}_{\infty}\}$. We have:

$$ \begin{aligned} (\tilde{K}_{\infty})^* &= \{(v, w) \in \mathbb{R}^{n \times (d-1)m} | v^\top x + w^\top y \ge 0, \text{ for all } x \in \mathcal{T}_{S_{n,u}^{st}}(q) \text{ and for all } y \in \mathbb{R}^{(d-1)m}\} \\ &= \{(v, w) \in \mathbb{R}^{n \times (d-1)m} | v^\top x \ge 0 \text{ for all } x \in \mathcal{T}_{S_{n,u}^{st}}(q), w = 0\} \\ &= (\mathcal{T}_{S_{n,u}^{st}}(q))^* \times \{0_{(d-1)m}\} \\ &= -\mathcal{N}_{S_{n,u}^{st}}(q) \times \{0_{(d-1)m}\}. \end{aligned} \quad (48) $$

It is noteworthy that the normal and tangent cones have to be given a precise meaning, being the set $S_{n,u}^{st}$ non convex in general, see Appendix C. The last equality is obtained because we assume that the set $S_{n,u}^{st}$ is prox-regular, so that the dual cone of its tangent cone is the normal cone, and both cones are understood in the sense of Clarke. Consequently $(x, y) \in \mathcal{K}(M_u^{st}(q), \psi_{\tilde{K}}) \Leftrightarrow (M(q)x - H_{t,u}^{st}(q)y) \in (\tilde{K}_{\infty})^* \Leftrightarrow M(q)x - H_{t,u}^{st}(q)y \in -\mathcal{N}_{S_{n,u}^{st}}(q)$ and $x \in \ker(H_{t,u}^{st}(q)^{\top})$. Therefore $(\mathrm{dom}(\varphi))_{\infty} \cap \ker\{\mathbf{M} + \mathbf{M}^{\top}\} \cap \mathcal{K}(\mathbf{M}, \varphi) = \{(x, y) \in \mathbb{R}^n \times \mathbb{R}^{(d-1)m} | x \in \mathcal{T}_{S_{n,u}^{st}}(q) \cap \ker(M(q)) \cap \ker(H_{t,u}^{st}(q)^{\top}), y \in (H_{t,u}^{st}(q))^{-1}(\mathcal{N}_{S_{n,u}^{st}}(q))\} = \mathcal{E}_{n,u}^{st}(q)$, where we used that $x \in \ker(M(q))$, see the kernel in the second set calculation. Item (i) is proved. Let us examine now item (ii). Let $(x_0, y_0) \in \mathrm{dom}(\psi_{\tilde{K}}) = \tilde{K} \Leftrightarrow x_0 \in K$ and $y_0 \in \mathbb{R}^{(d-1)m}$. The condition in (83) is equal to: $\langle F_u^{st}(q, \dot{q}, t) -$
$M_u^{st}(q) \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}, v \rangle + \psi_{\bar{K}_\infty}(v) > 0$ for all $v = \begin{pmatrix} v_n \\ v_m \end{pmatrix} \in \mathcal{E}_{n,u}^{st}(q)$. Notice that $\psi_{\bar{K}_\infty}(v) = \psi_{\mathcal{T}_{S_{n,u}^{st}}(q)}(v_1)$, using Proposition 15 c) and e). Developing the product we find $F_u^{st}(q, \dot{q}, t)^\top v - x_0^\top M(q)v_n - y_0^\top H_{t,u}^{st}(q)^\top v_n + x_0^\top H_{t,u}^{st}(q)v_m > 0$ for all $v_n \in \mathcal{T}_{S_{n,u}^{st}}(q) \cap \ker(M(q)) \cap \ker(H_{t,u}^{st}(q)^\top)$ and for all $v_2 \in (H_{n,u}^{st})^{-1}(\mathcal{N}_{S_{n,u}^{st}}(q))$. Thus equivalently $F_u^{st}(q, \dot{q}, t)^\top v + x_0^\top H_{t,u}^{st}(q)v_2 > 0$ for all $v_m \in (H_{t,u}^{st})^{-1}(\mathcal{N}_{S_{n,u}^{st}}(q))$. Item (iii) follows directly from Proposition 16 c). ■

The definition of the set $\mathcal{E}_{n,u}^{st}(q)$ shows that the proposed criterion is an extension of the bilateral case in Proposition 1. It is worth noting that solvability means that there exists $\ddot{q}(t)$, $\lambda_{t,u}^{st}(t)$, and $\lambda_{n,u}^{st}(t)$ which solve (43) $(\nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st}(t)$ being a selection of the set $-\mathcal{N}_{\bar{K}}(z(t))$ in (45)). Nothing can be deduced about the uniqueness of the multiplier $\lambda_{t,u}^{st}$ in Proposition 8 (iii), which is in agreement with the well-known fact that contact forces may be non unique in frictional problems [19, Section 5.5.6]. However acceleration non-uniqueness is essentially related to the mass matrix singularity. One sees that if there is no friction $(H_{t,u}^{st}(q) = 0)$ and if the constraints are bilateral so that $S_{n,u}^{st}$ is a codimension $m$ submanifold (the tangent cone thus reduces to the usual tangent space), we recover the criterion in [36]. The next result follows directly from Proposition 8.

**Corollary 2** Let the conditions of Proposition 8 hold. If $M(q) > 0$ and $\mathrm{im}(H_{t,u}^{st}(q)) \cap \mathcal{N}_{S_{n,u}^{st}}(q) = \{0_n\}$, then the contact MLCP in (43) is always solvable with uniqueness for $\ddot{q}(t)$.

**Example 4 (Example 1, continued)** Assume that the set $S_{n,u}^{st}$ is prox-regular and satisfies the MFCQ, the normal cone is generated by the normals at the active contact points and is calculated as:

$$ \mathcal{N}_{S_{n,u}^{st}}(q^{st}) = \{z \in \mathbb{R}^3 | z = \lambda_1 \nabla h_{n,u,1}^{st}(q^{st}) + \lambda_2 \nabla h_{n,u,2}^{st}(q^{st}), \lambda_1 \ge 0, \lambda_2 \ge 0, 0 \le \lambda \perp h_{n,u}(q) \ge 0\} \\ = \{z \in \mathbb{R}^3 | z = (0, \lambda_1 + \lambda_2, \frac{L}{2}(\lambda_1 - \lambda_2))^\top, \lambda_1 \ge 0, \lambda_2 \ge 0, 0 \le \lambda \perp h_{n,u}(q) \ge 0\} \quad (49) $$

The set $\mathrm{im}(H_{n,u}^{st}(q^{st})) = \{z \in \mathbb{R}^3 | z = \beta(1, 0, \frac{L}{2})^\top, \beta \in \mathbb{R}\}$. Thus any vector in $\mathrm{im}(H_{n,u}^{st}(q^{st})) \cap \mathcal{N}_{S_{n,u}^{st}}(q^{st})$ must satisfy $\lambda_1 + \lambda_2 = 0 \Rightarrow \lambda_1 = \lambda_2 = 0$, and $\beta = 0$. Corollary 2 applies.

Besides a full rank inertia matrix, the corollary requires some kind of uncoupling between tangential and normal directions. It is noteworthy that the sign conditions in the outwards normal cone's elements, reflect the unilaterality of the contact. Recall that $K$ is defined in (44).

**Corollary 3** Let the conditions of Proposition 8 hold. Assume that $\{0_n\} \in K$. Then the contact MLCP in (43) is solvable provided that $F_u^{st}(q, \dot{q}, t)^\top v > 0$ for all $v = (v_n^\top, v_m^\top)^\top \ne (0^n, 0^{(d-1)m})^\top$, $v_n \in \mathbb{R}^n$, $v_m \in \mathbb{R}^{(d-1)m}$, and $v_m \in (H_{n,u}^{st})^{-1}(\mathcal{N}_{S_{n,u}^{st}}(q))$.

Notice that $\{0_n\} \in K$ is guaranteed if $\frac{d}{dt}(\nabla h_{n,u}^{st}(q)^\top)\dot{q} = 0$, which occurs in statics ($\dot{q} = 0$) or if $\nabla h_{n,u}^{st}(q)$ is a constant.

**Remark 3** We could choose to work with a slightly transformed equivalent problem involving the symmetric matrix $M_u^{st}(q) = \begin{pmatrix} M(q) & -H_{t,u}^{st}(q) \\ -H_{t,u}^{st}(q)^\top & 0 \end{pmatrix}$. This matrix is never positive semidefinite as this would require that $-H_{t,u}^{st}(q)^\top M(q)^\top H_{t,u}^{st}(q) \ge 0$ [12, Proposition 8.2.3], which is clearly impossible
excepted if $H_{t,u}^{st}(q)$ is the zero matrix (the class of quasi-positive-semi-definite (QPSD) matrices [24], a relaxed form of QPD matrices [71], which enjoy block-diagonal factorizations, may be a future path of investigations). This is the reason why we chose to work with (43) instead, with a matrix with a skew-symmetric part, enabling us to apply Proposition 16. This is a big difference with the bilaterally constrained case (26), which does not rely on any positivity of the matrix $\bar{M}_b(q)$, but just on linear algebra arguments, so that the symmetric form in (26) is usually chosen [75, 55, 27, 76]. But unilaterality yields complementarity problems which need some positivity. However skew-symmetry hampers the use of items b'), d) and e) of Proposition 16 in Appendix B.

## 3.3 Adding the Coulomb's cone constraint

In the above two sections 3.1 and 3.2, we have derived the problems to be studied for the all-sticking mode, and some well-posedness results for the contact problem, assuming that the contact force is not constrained. Consequently the above results can be seen as necessary conditions for the complete contact problem to be well-posed, in the sense that if the condition of Corollary 2 applies, one still has to check that the solution satisfies Coulomb's constraint. However, if there is no solution to the above problem, then the complete problem has no solution. Let us now formulate the complete problem, that includes Coulomb's constraint $\Lambda_{t,i}^{st} = \begin{pmatrix} \lambda_{t,1,i}^{st} \\ \lambda_{t,2,i}^{st} \end{pmatrix} \in \mathcal{D}(\mu_i, \lambda_{n,u,i})$ at all the contact points $i$. Let us assume that the contacts have been ordered such that $\lambda_{t,u} = (\lambda_{t,u,1,1} \ \lambda_{t,u,1,2} \ \lambda_{t,u,2,1} \ \lambda_{t,u,2,2} \dots \lambda_{t,u,m,1} \ \lambda_{t,u,m,2})^\top = (\Lambda_{t,u,1} \ \Lambda_{t,u,2} \dots \Lambda_{t,u,m})^\top$. Coulomb's constraint can therefore be written as $\lambda_{t,u}^{st} \in \mathcal{D}(\mu, \lambda_{n,u}^{st})$ with

$$ \mathcal{D}(\mu, \lambda_{n,u}^{st}) \triangleq \mathcal{D}(\mu_1, \lambda_{n,u,1}^{st}) \times \dots \times \mathcal{D}(\mu_m, \lambda_{n,u,m}^{st}). \qquad (50) $$

Proposition 8 makes minimal assumptions on the system, however it does not provide many informations on the solutions. Therefore suppose that Assumption 1 holds, securing that $M_c^{st}(q)^{-1}$ in (36) is defined. Then using (10) (c) and (12) the complete all-sticking problem boils down to analysing whether or not the problem:

$$ \left\{ \begin{array}{l} -\lambda_{t,u}^{st} = (H_{t,u}^{st}(q))^{\top} M(q)^{-1} H_{t,u}^{st}(q))^{-1} R_{t,u}^{st}(q, \dot{q}, t, \lambda_{n,u}^{st}) \in -\mathcal{D}(\mu, \lambda_{n,u}^{st}) \\ 0 \le \lambda_{n,u}^{st} \perp A_c^{st}(q)\lambda_{n,u}^{st} + H_c^{st}(q, \dot{q}, t) \ge 0, \end{array} \right. \qquad (51) $$

is solvable with unknown $\lambda_{n,u}^{st}$, possibly with uniqueness. It follows from the expression of $R_{t,u}^{st}(q, \dot{q}, t, \lambda_{n,u}^{st})$ in (34) that $\lambda_{n,u}^{st}$ is multiplied by the inertial normal/tangential coupling matrix $H_{t,u}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q)$. Using (51) we can state the following result.

**Proposition 9 (Frictional, sticking unilateral constraints with Coulomb's cone)** Let Assumption 1 hold.

(i) The mode corresponding to sticking unilateral contacts has a solution, if and only if the prob-
lem:

$$ egin{gathered} \left( H_{t,u}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q) \quad \lambda_{n,u}^{st} + \left( H_{t,u}^{st}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) - \frac{d}{dt}(H_{t,u}^{st}(q)^{\top}) \dot{q} \right) H_c^{st}(q, \dot{q}, t) \right) \\ \in - \left( \begin{array}{cc} H_{t,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q) & 0 \\ 0 & I \end{array} \right) \left( \begin{array}{c} \mathcal{D}(\mu, \lambda_{n,u}^{st}) \\ \mathcal{N}_{\mathbb{R}_+^{m_u^{st}}}(\lambda_{n,u}^{st}) \end{array} \right) \end{gathered} \tag{52} $$

is solvable.

(ii) Assume now that $A_c^{st}(q)$ in (35) is positive definite. Then the mode corresponding to sticking unilateral contacts has a solution, if and only if the problem:

$$ \left\{ \begin{array}{l} H_{t,u}^{st}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) - \frac{d}{dt} (H_{t,u}^{st}(q)^{\top}) \dot{q} + H_{t,u}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q) \lambda_{n,u}^{st} \\ \in -H_{t,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q) \mathcal{D}(\mu, \lambda_{n,u}^{st}) \quad (53) \\ \lambda_{n,u}^{st} = \text{proj}_{A_c^{st}(q)}[\mathbb{R}_+^{m_u^{st}}; -A_c^{st}(q)^{-1} H_c^{st}(q, \dot{q}, t)] \end{array} \right. $$

is solvable, where $H_c^{st}(q, \dot{q}, t)$ is in (38).

**Proof 8** (i) follows from (51) and the definition of $R_{t,u}^{st}(q, \dot{q}, t, \lambda_{n,u}^{st})$ in (12). Item (ii) follows from (i) and application of (86). ■

Assumption 1 is needed to state (52) because of (38) and (36). Imposing $A_c^{st}(q) \succ 0$ in (53) is even more stringent. The problem (52) is a generalized equation with unknown $\lambda_{n,u}^{st}$, i.e., an equation of the form $0 \in f(\lambda_{n,u}^{st}) + F(\lambda_{n,u}^{st})$, with $f(\cdot)$ single-valued, $F(\cdot)$ multivalued functions, where the external action, $q, \dot{q}$ and $\mu$ are given data. The problem in (53) is a generalized equation where $q, \dot{q}$ and possibly $\mu$ and the external action, are the unknowns, since the multiplier $\lambda_{n,u}^{st}$ is calculated explicitly as a function of $q, \dot{q}$ and external actions.

The solvability of the generalized equations (53) (52) may have several meanings: given $q, \dot{q}$ and $\mu_i$, find an external force $F(q, \dot{q}, t)$; or given $F(q, \dot{q}, t)$ (inertial forces and external forces applied on the system), find the state $q, \dot{q}$ and friction coefficients $\mu_i$; or find external forces, $q, \dot{q}$ and $\mu$. Uniqueness usually fails. The problem tackled here is similar to the problem studied in [54, Section 4]. Instead of assuming that $A_c^{st}(q)$ in (35) is positive definite, we may simply suppose that the LCP is solvable, using Proposition 5 (iv). Then if one solution can be computed, it suffices to check whether or not the first line of (53) is satisfied.

Let us describe a method for solving (53) when the sets $\mathcal{D}_{\mu_i}$ are facetized with an inner approximation as depicted in Figure 4, where we recall that at each contact point $\mathcal{D}(\lambda_{n,i}, \mu_i) = |\lambda_{n,i}|\mathcal{D}_{\mu_i}$ and $\mathcal{D}_{\mu_i} = \{z \in \mathbb{R}^{d-1} | z^\top z \le \mu_i^2\}$. The facetized sets are convex polyhedral sets denoted $\bar{\mathcal{D}}_{\mu_i} = \{w \in \mathbb{R}^{d-1} | C_{\mu_i}w + D_{\mu_i} \ge 0\}$ for $C_{\mu_i} \in \mathbb{R}^{l \times (d-1)}$ and $D_{\mu_i} \in \mathbb{R}^l$. The integer $l$ depends on the number of facets of the approximation. This operation makes sense only if $d = 3$ (3D friction), for otherwise the
<figure><img src="image_7.png" /><figcaption>Figure 4: Inner polyhedral approximation of Coulomb-Moreau’s disk.</figcaption></figure>

Coulomb-Moreau disk is already facetized, being an interval. In order to prepare the next proposition, let us first rewrite more compactly (53). Let us denote

$$ \tilde{F}(q, \dot{q}, t) \triangleq H_{t,u}^{st}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) - \frac{d}{dt} (H_{t,u}^{st}(q)^{\top}) \dot{q} + H_{t,u}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q) \lambda_{n,u}^{st}(q, \dot{q}, t) $$

and

$$ A_{tt}^{st}(q) \triangleq H_{t,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q) \in \mathbb{R}^{(d-1)m \times (d-1)m}. \qquad (54) $$

We further define $\tilde{A}_{tt}^{st}(q, \dot{q}, t) \triangleq A_{tt}^{st}(q)$ blockdiag $\begin{bmatrix} \lambda_{n,u,i}^{st} & 0 \\ 0 & \lambda_{n,u,i}^{st} \end{bmatrix} \in \mathbb{R}^{(d-1)m \times (d-1)m}$, where we recall that from (53) $\lambda_{n,u}^{st} = \lambda_{n,u}^{st}(q, \dot{q}, t)$, being the projection of a known quantity. Thus the inclusion in (53) is rewritten in its facetized form as $\tilde{F}(q, \dot{q}, t) + \tilde{A}_{tt}^{st}(q, \dot{q}, t) z = 0$ where

$$ z = (z_{11} \ z_{12}, z_{21} \ z_{22}, \dots, z_{i1} \ z_{i2}, \dots, z_{m1} \ z_{m2})^{\top} \in \mathbb{R}^{(d-1)m}, \text{ and } z_i = \begin{pmatrix} z_{i1} \\ z_{i2} \end{pmatrix} \in \bar{D}_{\mu_i}. $$

Finally we set $C_{\mu} \triangleq \text{blockdiag}(C_{\mu_i}) \in \mathbb{R}^{ml \times (d-1)m}$, $D_{\mu} \triangleq (D_{\mu_1}^{\top}, \dots, D_{\mu_m}^{\top})^{\top}$, $\bar{D}_{\mu} \triangleq \bar{D}_{\mu_1} \times \dots \times \bar{D}_{\mu_m}$. Hence $\bar{D}_{\mu} = \{z \in \mathbb{R}^{(d-1)m} | C_{\mu}z + D_{\mu} \ge 0\} = P(C_{\mu}, D_{\mu})$ and the facetized inclusion becomes

$$ \tilde{F}(q, \dot{q}, t) = -\tilde{A}_{tt}^{st}(q, \dot{q}, t) z, \quad z \in \bar{D}_{\mu}. \qquad (55) $$

It is noteworthy that both $\bar{D}_{\mu}$ and $D_{\mu}$ are non-empty convex bounded closed (hence compact) sets.

**Proposition 10** Let $q, \dot{q}$ and external actions be given. Assume that $z^*$ is a solution of (55), then $z^*$ also solves the inclusion

$$ \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t) z^* + \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{F}(q, \dot{q}, t) \in -\mathcal{N}_{\bar{D}_{\mu}}(z^*). \qquad (56) $$
**Proof 9** Assume that $z^*$ solves (55). Then $z^*$ is a solution of the quadratic program:

$$ \left\{ \begin{array}{l} \min_z \frac{1}{2}(\tilde{F}(q, \dot{q}, t) + \tilde{A}_{tt}^{st}(q, \dot{q}, t) z)^{\top}(\tilde{F}(q, \dot{q}, t) + \tilde{A}_{tt}^{st}(q, \dot{q}, t) z) \\ \text{subject to: } z \in \bar{D}_{\mu}. \end{array} \right. \qquad (57) $$

From the expression of $\bar{D}_{\mu}$, the (necessary and sufficient) KKT conditions for the problem (57) are given by:

$$ \left\{ \begin{array}{l} \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t) z^{*} + \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{F}(q, \dot{q}, t) - C_{\mu}^{\top} \lambda = 0 \\ 0 \le \lambda \perp C_{\mu} z^{*} + D_{\mu} \ge 0, \end{array} \right. \qquad (58) $$

which is an MLCP. We can rewrite equivalently (58) as $\tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t) z^{*} + \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{F}(q, \dot{q}, t) \in -C_{\mu}^{\top} \partial \psi_{\mathbb{R}^{m+}}(C_{\mu} z^{*} + D_{\mu}).$ Since the indicator function of a polyhedron is a polyhedral function [35], applying the chain rule [21, Proposition A.2], we obtain equivalently (56). ■

The interest of (56) is that the matrix $\tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t) \succeq 0$, and it is symmetric. It is noteworthy that $A_{tt}^{st}(q)$ in (54) is $\succ 0$ in case (ii), because $A_{tt}^{st}(q) \succ 0$ if and only if all constraints are independent, using the Schur complement Theorem. Consequently the rank of $\tilde{A}_{tt}^{st}(q, \dot{q}, t)$, and in turn of $\tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t)$, depends on the rank of blockdiag $\begin{bmatrix} \lambda_{n,u,i}^{st} & 0 \\ 0 & \lambda_{n,u,i}^{st} \end{bmatrix}$ only. We can apply Proposition 16 to analyse the well-posedness of the following variational inequality, which is equivalent to (56):

$$ \begin{array}{l} \text{Find } z^* \in \bar{D}_{\mu} \text{ such that:} \\ \langle \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t) z^* + \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{F}(q, \dot{q}, t), y - z^* \rangle \ge 0, \text{ for all } y \in \bar{D}_{\mu}. \end{array} \qquad (59) $$

The variable $z^*$ is a selection of the facetized Coulomb-Moreau's disk, which guarantees that the active contacts are tangentially sticking.

**Corollary 4 (i)** The variational inequality in (59) is always solvable; (ii) If $z_1^*$ and $z_2^*$ are two solutions then $z_1^* - z_2^* \in \ker(\text{blockdiag} \left[ \begin{pmatrix} \lambda_{n,u,i}^{st} & 0 \\ 0 & \lambda_{n,u,i}^{st} \end{pmatrix} \right])$; (iii) Moreover $\langle \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{F}(q, \dot{q}, t), z_1^* - z_2^* \rangle = 0$.

**Proof 10 (i)** Let us first calculate the set $(\text{dom}(\varphi))_{\infty} \cap \ker\{\mathbf{M}\} \cap \mathcal{K}(\mathbf{M}, \varphi)$ in Proposition 16. We have $\mathbf{M} = \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t)$, $\varphi = \psi_{\bar{D}_{\mu}}$, hence $\text{dom}(\varphi) = \bar{D}_{\mu}$. Using Proposition 15 e) f) we infer that $(\text{dom}(\varphi))_{\infty} = (\bar{D}_{\mu})_{\infty} = \{z \in \mathbb{R}^{(d-1)m} | C_{\mu}z \ge 0\} = P(C_{\mu}, 0) = \{0\}$, because $\bar{D}_{\mu}$ is a bounded closed convex set. Thus (i) follows from Proposition 16 a). Secondly $\ker(\tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t)) = \ker(\tilde{A}_{tt}^{st}(q, \dot{q}, t))$ [12, Theorem 2.4.3]. As we have seen above, under the conditions for (53) we have $A_{tt}^{st}(q) \succ 0$. Therefore $\ker(\tilde{A}_{tt}^{st}(q, \dot{q}, t)) = \ker(\text{blockdiag} \left[ \begin{pmatrix} \lambda_{n,u,i}^{st} & 0 \\ 0 & \lambda_{n,u,i}^{st} \end{pmatrix} \right])$. Thus (ii) follows from Proposition 16 c), using once again that $A_{tt}^{st}(q)$ is full rank. (iii) is a consequence of Proposition 16 d). ■
Notice that the uniqueness condition in Proposition 16 e) cannot hold. Indeed $\langle \tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{F}(q, \dot{q}, t), e \rangle = \tilde{F}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t)e$ and this is equal to zero for all $e \in \ker(\tilde{A}_{tt}^{st}(q, \dot{q}, t)^{\top} \tilde{A}_{tt}^{st}(q, \dot{q}, t)) = \ker(\tilde{A}_{tt}^{st}(q, \dot{q}, t))$. We could have followed a similar path without facetization and working with $D_{\mu}$. However the facetization paves the way to the use of efficient numerical solvers for MLCPs as (58).

**Recapitulation:** $z^*$ solution of (53) $\Leftrightarrow$ $z^*$ solution of (55) $\Leftrightarrow$ if $\min = 0$ $z^*$ solution of (57) $\Leftrightarrow$ $z^*$

solution of (58) $\Leftrightarrow$ $z^*$ solution of (56) $\Leftrightarrow$ $z^*$ solution of (59). Thus in practice Corollary 4 provides sufficient conditions for the solvability of (59), and the QP (57) can be numerically solved. If the obtained solutions imply that the quadratic objective function vanishes, the all-sticking problem is solved.

What about problem (52) in Proposition 9 (i)? Proposition 6 can be used to check whether the contact LCP (37) (38) is solvable with uniqueness. Getting an analytical expression for one solution $\lambda_{n,u}^{st}$ with $A_c^{st}(q) \ge 0$ only, is not doable in general. One can nevertheless numerically compute one solution and then inject its numerical value in the problem. The only missing step compared with (53), is that the projection is not available.

The most general problem with no rank assumption, consists of the MLCP (43) with the constraint $\lambda_{t,u}^{st} \in D(\mu, \lambda_{n,u}^{st})$. One thus obtains the following problem, where no assumptions are made a priori on the rank of matrices:

$$ \left\{ \begin{array}{l} (a) M(q)\ddot{q} + F(q, \dot{q}, t) - \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} - H_{t,u}^{st}(q)\lambda_{t,u}^{st} = 0 \\ (b) H_{t,u}^{st}(q)^{\top}\ddot{q} + \frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q} = 0 \\ (c) 0 \le \lambda_{n,u}^{st} \perp a_{n,u} = \nabla h_{n,u}^{st}(q)^{\top}\ddot{q} + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \ge 0 \\ (d) \lambda_{t,u}^{st} \in \text{blockdiag} \begin{pmatrix} \lambda_{n,u,i}^{st} & 0 \\ 0 & \lambda_{n,u,i}^{st} \end{pmatrix} D_{\mu}. \end{array} \right. \quad (60) $$

This is a mixed nonlinear CP, with unknowns $\ddot{q}, \lambda_{n,u}^{st}, \lambda_{t,u}^{st}$, and data $\mu, q, \dot{q}$, and external actions. The major discrepancy with respect to (52) and (53), is that we have no means to first calculate $\lambda_{n,u}^{st}$ in (60).

**Proposition 11** Assume that the conditions of solvability in Proposition 8 (i) or (ii) are satisfied. Then there exists $\mu_{\min} < +\infty$ such that problem (60) (a)-(d) has at least a solution for all $\mu \ge \mu_{\min}$.

**Proof 11** The MLCP (43) (which is (60) (a) (b) (c)) being solvable, one can find at least one triplet $(\ddot{q}(t), \lambda_{t,u}^{st}, \lambda_{n,u}^{st})$ which solves (60) (a) (b) (c). The set $D_{\mu}$ can be enlarged increasing $\mu_i$'s so that the inclusion (60) (d) is satisfied. ■

In particular, $\mu = 0$ implies $\lambda_{t,u}^{st} = 0$, so that solving the MLCP (60) (a) (b) (c) is a test for sticking contact existence with no friction (called the frictionless stability in [52]).

## 3.4 Another way to tackle the problem

If we don't want to allow for the possibility of $a_{n,u,i}(t^+) > 0$ for some $i$ (a possible outcome of the foregoing problem, which involves complementarity conditions in (43) and in (60)), then we may
state right from the beginning that $a_{n,u}(t^+) = 0$ and construct a linear algebraic problem similar to (14) (a), adapting the matrices to cope with the unilateral constraint case. One big difference with the bilateral constraints case, is however that we have to check that the solution satisfies $\lambda_{n,u}^{st}(t) \ge 0$ (or $> 0$). The complete problem thus takes the form:

$$ (a) \quad \begin{pmatrix} M(q) & -\nabla h_{n,u}^{st}(q) & -H_{t,u}^{st}(q) \\ \nabla h_{n,u}^{st}(q)^{\top} & 0 & 0 \\ H_{t,u}^{st}(q)^{\top} & 0 & 0 \end{pmatrix} \begin{pmatrix} \ddot{q}(t) \\ \lambda_{n,u}^{st}(t) \\ \lambda_{t,u}^{st}(t) \end{pmatrix} = \begin{pmatrix} -F(q, \dot{q}, t) \\ -\frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \\ -\frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q} \end{pmatrix} \quad (61) $$
$$ (b) \quad \lambda_{n,u}^{st}(t) \ge 0 $$

$$ (c) \quad \lambda_{t,u}^{st} \in \text{blockdiag} \begin{pmatrix} \lambda_{n,u,i}^{st} & 0 \\ 0 & \lambda_{n,u,i}^{st} \end{pmatrix} \quad \mathcal{D}_{\mu}. $$

All the results obtained in section 2 can be repeated for (61) (a) and (c). One can invert (61) (a) as done in (30) (31) and check whether or not (61) (b) and (c) hold. In general this will be possible only numerically. The great advantage of working with complementarity conditions, is that it allows one to formulate conditions for the existence of solutions. This is not the case of problem (61) which is basically a linear algebra problem whose solutions have to satisfy a signum constraint. Notice that it is always possible to solve the problem allowing for $a_{n,u,i}(t^+) > 0$, and then eliminate unwanted situations where some contacts detach, using a suitable numerical test. It is noteworthy that we have no general criterion for the existence of solutions to (61) (a) (b).

**Example 5** Let us consider the manipulation problem depicted in Figure 5, where three “fingers” with prismatic joints are used. All coordinates are in the Galilean frame $\mathcal{L}_0 = (0, i, j)$. The gravity center G coordinates are $(x, y)$, each finger tip $A_i$ has coordinates $(x_i, y_i)$, $1 \le i \le 3$, with $y_1, x_2$ and $y_3$ constant, $B_i$ is the projection of $A_i$ onto the closest the segment of the rectangular object with width L and height l, A is the right bottom corner with coordinates $(x_A, y_A)$, the local frames $\mathcal{L}_i$ at the potential contact points are $(B_i, \mathbf{n}_i, \mathbf{t}_i)$, to simplify we denote $\mathbf{n} = (\cos(\theta), \sin(\theta))^{\top}$ and $\mathbf{t} = (\sin(\theta), -\cos(\theta))^{\top}$ for the normal and tangential unit vectors as indicated in the Figure. The mass matrix is $M = \text{diag}(m, m, I_G, m_1, m_2, m_3)$, and $q = (x, y, \theta, x_1, y_2, x_3)^{\top} \in \mathbb{R}^6$, $n = 6$, $m = m_u = 3$, $d = 2$. We also let $\beta(x, y, \theta) \triangleq -\sin(\theta)y_A - \cos(\theta)x_A$. The gap functions $h_{n,u,i}(q)$ are computed from $h_{n,u,i}(q) = (A_i B_i)^{\top} \mathbf{n}_i$, thus $h_{n,u,3}(q) = (A_3 B_3)^{\top} \mathbf{n}$. The projection $B_3$ may be calculated from the KKT conditions which yield $\min_{a,b} \frac{1}{2}(a - x_3)^2 + \frac{1}{2}(b - y_3)^2$ subject to $\sin(\theta)b + \cos(\theta)a + \beta(x, y, \theta) = 0$, so that $B_3 = (a, b)^{\top}$ with $x_{B_3} = a = x_3 - \cos^2(\theta)x_3 - \beta(x, y, \theta)\cos(\theta) - \cos(\theta)\sin(\theta)y_3$, $y_{B_3} = b = y_3 - \sin^2(\theta)y_3 - \beta(x, y, \theta)\sin(\theta) - \sin(\theta)\cos(\theta)x_3$. Consequently $h_{n,u,3}(q) = (A_3 B_3)^{\top} \mathbf{n} = (x_3 - x_A(x, y, \theta))\cos(\theta) + (y_3 - y_A(x, y, \theta))\sin(\theta)$, $x_A = x + \frac{1}{2}\sqrt{L^2 + l^2}\sin(\alpha + \theta)$, $y_A = y - \frac{1}{2}\sqrt{L^2 + l^2}\cos(\alpha + \theta)$, thus $h_{n,u,3}(q) = (x_3 - x + \frac{1}{2}\sqrt{L^2 + l^2}\sin(\alpha + \theta))\cos(\theta) + (y_3 - y + \frac{1}{2}\sqrt{L^2 + l^2}\cos(\alpha + \theta))\sin(\theta)$. Performing rotations of the Galilean frame it is inferred that $h_{n,u,2}(q) = (y_2 - y + \frac{1}{2}\sqrt{L^2 + l^2}\sin(\alpha + \theta))\cos(\theta) + (-x_2 + x + \frac{1}{2}\sqrt{L^2 + l^2}\cos(\alpha + \theta))\sin(\theta)$, $h_{n,u,1}(q) = (-x_1 + x + \frac{1}{2}\sqrt{L^2 + l^2}\sin(\alpha + \theta))\cos(\theta) + (-y_1 + y + \frac{1}{2}\sqrt{L^2 + l^2}\cos(\alpha + \theta))\sin(\theta)$. The gradient $\nabla h_{n,u}(q)$ in the right-hand side of (1a) is deduced easily. Notice that $\frac{d}{dt}h_{n,u,i}(q) = (\frac{d}{dt}(A_i B_i)^{\top}) \mathbf{n} + (A_i B_i)^{\top} \frac{dn}{dt}$. However $\frac{dn}{dt}|_{\mathcal{L}_0} = \frac{dn}{dt}|_{\mathcal{L}_i} + \Omega_{\mathcal{L}_0/\mathcal{L}_i} \times \mathbf{n}$, hence $(A_i B_i)^{\top} \frac{dn}{dt} = 0$ because $\mathbf{n}$ is fixed in $\mathcal{L}_i$, while $\Omega_{\mathcal{L}_0/\mathcal{L}_i} \times \mathbf{n}$ is orthogonal to $\mathbf{n}$ and $A_i B_i$ is colinear
to $\mathbf{n}^5$. Hence $\frac{d}{dt}h_{n,u,i}(q) = (\frac{d}{dt}(A_i B_i)^\top)\mathbf{n} = v_{n,u,i} = \nabla h_{n,u,i}(q)^\top \dot{q}$ (all vectors being expressed in $\mathcal{L}_0$). In the same way $a_{n,u,i} = (\frac{d^2}{dt^2}(A_i B_i)^\top)\mathbf{n} + (\frac{d}{dt}(A_i B_i)^\top)(\Omega_{\mathcal{L}_0/\mathcal{L}_i} \times \mathbf{n}) = \nabla h_{n,u,i}(q)^\top \ddot{q} + \frac{d}{dt}(\nabla h_{n,u,i}(q)^\top)\dot{q}$. Let us now deal with relative tangential velocities at contact points when the contacts are activated. The relative velocity at contact 3 is given by $V_{r,3} \triangleq V_{A_3} - V_{B_3}$ where both vectors are expressed in $\mathcal{L}_0$. We have $V_{A_3} = (\dot{x}_3, \dot{y}_3, 0)^\top$, $V_{B_3} = (\dot{x}_3 + 2\theta \cos(\theta) \sin(\theta)x_3 - \cos^2(\theta)\dot{x}_3 - \dot{\beta} \cos(\theta) + \dot{\theta}\beta \sin(\theta) + \dot{\theta}\sin^2(\theta)y_3 - \dot{\theta}\cos^2(\theta)y_3 - \cos(\theta)\sin(\theta)\dot{y}_3, \dot{y}_3 - 2\dot{\theta}\cos(\theta)\sin(\theta)y_3 - \sin^2(\theta)\dot{y}_3 - \dot{\beta}\sin(\theta) - \beta\dot{\theta}\cos(\theta) - \dot{\theta}\cos^2(\theta) + \dot{\theta}\sin^2(\theta)x_3 - \sin(\theta)\cos(\theta)\dot{x}_3, 0)^\top$, with $\dot{y}_3 = 0$. One obtains the following: $v_{r,t,3} = V_{r,3}^\top \mathbf{t} = (0, 0, -\cos(\theta)x_3 - \sin(\theta)y_3 - \beta(x, y, \theta), 0, 0, 0) \dot{q}$. The same can be redone at contacts 1 and 2, and $\triangleq H_{t,u,3}^\top$

$H_{t,u}(q) = (H_{t,u,1}(q) \ H_{t,u,2}(q) \ H_{t,u,3}(q)) \in \mathbb{R}^{6 \times 3}.$
The external action at $G$ on the object is given by $F_{\text{ext/ob}}(t) = (F_{\text{ext/ob},x}(t), F_{\text{ext/ob},y}(t), C_{\text{ext/ob},\theta}(t), 0, 0, 0)^\top$, where $C_{\text{ext/ob},\theta}(t)$ is a torque which works on $\theta$. Each finger $i$ is acted upon by a force $F_{\text{ext}/i}(t)$, giving $F_{\text{ext},f}(t) = (0, 0, 0, F_{\text{ext}/1}(t), F_{\text{ext}/2}(t), F_{\text{ext}/3}(t))^\top$. Notice that due to the constant $M$, one has $F(q, \dot{q}, t) = F_{\text{ext/ob}}(t) + F_{\text{ext}/f}(t)$. The MLCP in (43) can be constructed and analysed using the above results. The first step in the analysis of the problem, is to apply Proposition 8. Since rank($M$) = 6, this boils down to checking whether or not the inclusion $y \in (H_{t,u}^{st}(q))^{-1}(\mathcal{N}_{S_{n,u}^{st}}(q))$ has the 6-vector $(0, 0, \dots, 0)^\top$ as its unique solution (item (i)). If not, check item (ii). Using the material in Appendix C, one sees that this inclusion can be transformed in a problem of the form: find $y$ and $\lambda_i \ge 0$ such that $H_{t,u}^{st}(q)y = - \sum_{i=1}^3 \lambda_i \nabla h_{n,u,i}(q)$. One feature of the MLCP in (43) is that $M_u^{st}(q)$ is a sparse matrix, which can easily be rendered symmetric for the purpose of numerical solvers. In case of constraints redundancy that imply loss of rank, efficient interior point algorithms can be used to solve the MLCP.

# 4 Unilaterally/bilaterally constrained systems (all-sticking contacts)

Let us now deal with the more general case with both bilateral and unilateral constraints, which gives rise to the mixed LCP in (10)-(13). As above we first study the case without Coulomb's cone constraints, then these constraints are inserted.

## 4.1 Augmented contact LCP

Consider (b) (c) (d) in (10) as well as (11) (12) and (13). It is possible to rewrite equivalently the three equations (b) (c) (d) in (10) as:

$$ \mathcal{H}_1(q) \begin{pmatrix} \lambda_{t,b}^{st} \\ \lambda_{t,u}^{st} \\ \lambda_{n,b}^{st} \end{pmatrix} + \mathcal{F}(q, \dot{q}, t) + \mathcal{H}_2(q) \lambda_{n,u}^{st} = 0. \quad (62) $$

<footnote n="5">⁵Vector product is calculated adding 0 as the third component of vectors.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 5: A robotic planar manipulation task.</figcaption></figure>

The matrix $\mathcal{H}_1(q)$ in the right-hand side of (62) has the three matrices:

$$H_{t,b}^{st}(q)^{\top} M(q)^{-1} H_{t,b}^{st}(q), H_{t,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q), \nabla h_{n,b}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,b}^{st}(q)$$

on its diagonal. The off-block-diagonal sub-matrices are made of dynamic couplings between the three multipliers $\lambda_{t,b}^{st}$, $\lambda_{t,u}^{st}$ and $\lambda_{n,b}^{st}$, i.e., terms of the form $H_{t,b}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q)$, $H_{t,b}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,b}^{st}(q)$, etc. One has

$$\mathcal{F}(q, \dot{q}, t) \triangleq \begin{pmatrix} -H_{t,b}^{st}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) + \frac{d}{dt} (H_{t,b}^{st}(q)^{\top}) \dot{q} \\ -H_{t,u}^{st}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) + \frac{d}{dt} (H_{t,u}^{st}(q)^{\top}) \dot{q} \\ -\nabla h_{n,b}^{st}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) + \frac{d}{dt} (\nabla h_{n,b}^{st}(q)^{\top}) \dot{q} \end{pmatrix}$$

$$\mathcal{H}_2(q) \triangleq \begin{pmatrix} H_{t,b}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q) \\ H_{t,u}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q) \\ \nabla h_{n,b}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q) \end{pmatrix}.$$

**Assumption 3** One has:

$$H_{t,b}^{st}(q)^{\top} M(q)^{-1} H_{t,b}^{st}(q) \succ 0, H_{t,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q) \succ 0, \nabla h_{n,b}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,b}^{st}(q) \succ 0.$$

**Lemma 2** Let Assumption 3 hold true. Consider the matrix $\mathcal{H}_1(q)$ in (62). If the dynamical couplings between the three multipliers $\lambda_{t,b}^{st}$, $\lambda_{t,u}^{st}$ and $\lambda_{n,b}^{st}$ are small enough, then $\mathcal{H}_1(q)$ is invertible.
The proof follows from diagonal dominance [12, Fact 4.10.15]. If the conditions of Lemma 2 are satisfied, (a) and (e) in (10) can be used to deduce that the contact LCP matrix is given this time by (compare with (35) (36) (37) (38)):

$$A_c^{st}(q) = \nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} [\nabla h_{n,u}^{st}(q) - (H_{t,b}^{st}(q) H_{t,u}^{st}(q) \nabla h_{n,b}^{st}(q)) \mathcal{H}_1(q)^{-1} \mathcal{H}_2(q)] \quad (63)$$

We may now use Corollary 5 to state conditions under which $A_c^{st}(q) > 0$.

**Proposition 12** Suppose that $M(q) > 0$, $A_{n,u}^{st}(q) \triangleq \nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q) > 0$. Assume that:

$$\sigma_{\max}(\nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} (H_{t,b}^{st}(q) H_{t,u}^{st}(q) \nabla h_{n,b}^{st}(q))) \sigma_{\max}(\mathcal{H}_2(q)) < \sigma_{\min}(A_{n,u}^{st}(q)) \sigma_{\min}(\mathcal{H}_1(q)). \quad (64)$$

Then the matrix $A_c^{st}(q)$ in (63) is positive definite, consequently the contact LCP that results from (10)-(13) always has a unique solution.

**Proof 12** The proof is a consequence of Corollary 5. ■

One sees that the inequality (64) boils down to designing systems with small enough crossed normal/tangential inertial couplings.

## 4.2 Analysis of acceleration inclusions

Starting from (9) and including Coulomb's cone constraints one obtains a direct extension of (60), as the following problem with unknowns $\ddot{q}(t)$, $\lambda_{n,b}^{st}(t)$, $\lambda_{t,b}^{st}(t)$, $\lambda_{t,u}^{st}(t)$, $\lambda_{n,u}^{st}(t)$:

$$\left\{
\begin{array}{l}
(a) \quad \underbrace{\begin{pmatrix} M(q) & -\nabla h_{n,b}^{st}(q) & -H_{t,b}^{st}(q) & -H_{t,u}^{st}(q) \\ \nabla h_{n,b}^{st}(q)^{\top} & 0 & 0 & 0 \\ H_{t,b}^{st}(q)^{\top} & 0 & 0 & 0 \\ H_{t,u}^{st}(q)^{\top} & 0 & 0 & 0 \end{pmatrix}}_{\triangleq M_{u,b}^{st}(q) \in \mathbb{R}^{(n+(d-1)m+m_b) \times (n+(d-1)m+m_b)}} \underbrace{\begin{pmatrix} \ddot{q}(t) \\ \lambda_{n,b}^{st}(t) \\ \lambda_{t,b}^{st}(t) \\ \lambda_{t,u}^{st}(t) \end{pmatrix}}_{\triangleq z(t)} + \underbrace{\begin{pmatrix} F(q, \dot{q}, t) \\ \frac{d}{dt}(\nabla h_{n,b}^{st}(q)^{\top})\dot{q} \\ \frac{d}{dt}(H_{t,b}^{st}(q)^{\top})\dot{q} \\ \frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q} \end{pmatrix}}_{\triangleq F_{u,b}^{st}} \\
\\
(b) \quad 0 \le \lambda_{n,u}^{st}(t) \perp a_{n,u}(t) = \nabla h_{n,u}^{st}(q)^{\top} \ddot{q}(t) + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \ge 0 \\
\\
(c) \quad \lambda_{t,u}^{st}(t) \in \text{blockdiag} \begin{pmatrix} \lambda_{n,u,i}^{st}(t) & 0 \\ 0 & \lambda_{n,u,i}^{st}(t) \end{pmatrix} \mathcal{D}_{\mu,u} \\
\\
(d) \quad \lambda_{t,b}^{st}(t) \in \text{blockdiag} \begin{pmatrix} \lambda_{n,b,i}^{st}(t) & 0 \\ 0 & \lambda_{n,b,i}^{st}(t) \end{pmatrix} \mathcal{D}_{\mu,b}.
\end{array}
\right.
\quad (65)$$
Let us recall that it is assumed that the problem is analysed at a time $t$ with active unilateral contacts, i.e., $h_{n,u}(q(t)) = 0$ and $\nabla h_{n,u}(q(t))^{\top} \dot{q}(t) = 0$, so that we cannot exclude $a_{n,u,i}(t^+) > 0$ for some contacts $i$. In the tangential direction, the system (65) (a) is constructed under the assumption that $v_{t,i}(t) = 0$ and $a_{t,i}(t) = 0$, for both bilateral and unilateral contacts: the constraints (65) (c) (d) only check the admissibility of the tangential multipliers which solve (65) (a) (b). The contacts keep tangential sticking provided the acceleration is continuous in a right-neighborhood of $t$, see Example 6 for an illustration. Without the Coulomb's cone constraints in (c) and (d), the problem (65) is an MLCP which can be analysed in a similar way as we did for (43) in Propositions 6, 7 and 8 and Corollaries 2 and 3. It suffices to replace $-H_{t,u}^{st}(q)$ in (43) by $(-\nabla h_{n,b}^{st}(q) - H_{t,b}^{st}(q) - H_{t,u}^{st}(q))$, and $\lambda_{t,u}^{st}$ by $\begin{pmatrix} \lambda_{n,b}^{st}(t) \\ \lambda_{t,b}^{st}(t) \\ \lambda_{t,u}^{st}(t) \end{pmatrix}$, etc. One obtains the variational inequality equivalent to (65) (a) (b): Find $z(t) \in \tilde{K}$ such that

$$ \langle M_{u,b}^{st}(q)z(t) + F_{u,b}^{st}(q, \dot{q}, t), y - z(t) \rangle \ge 0 \text{ for all } y \in \tilde{K}, \quad (66) $$

where $\tilde{K}$ is defined as in (45) with appropriate modifications. Adding the Coulomb's cone constraint (65) (c) (d) we can proceed as in Propositions 9 and 10. The following algorithm may be designed to solve (65):

**All-Sticking Bilateral/Unilateral Contact Algorithm:**

1. **data:** $q(t), \dot{q}(t), M(q), F(q, \dot{q}, t), h_{n,u}^{st}(q(t)), h_{n,b}^{st}(q(t)), H_{t,b}^{st}(q(t)), H_{t,u}^{st}(q(t)), \mu_i, 1 \le i \le m$.

2. **Compute** the set $S_z(t) = \{z(t) \in \mathbb{R}^{n+dm_b+(d-1)m_u} | z(t)$ solves the VI in (66)}$.

3. **Compute** the set $S_{\lambda_{n,u}^{st}}(t) = \{\lambda_{n,u}^{st} \in \mathbb{R}^{m_u} | \nabla h_{n,u}^{st}(q(t))\lambda_{n,u}^{st} = (M(q) - \nabla h_{n,b}^{st}(q) - H_{t,b}^{st}(q) - H_{t,u}^{st}(q)) z(t)$ for all $z(t) \in S_z(t)\}$.

4. **Compute** the sets $S_{\lambda_{t,u}^{st}}(t)$ of multipliers $\lambda_{t,u}^{st}(t)$ from (65) (c), and $S_{\lambda_{t,b}^{st}}(t)$ of multipliers $\lambda_{t,b}^{st}(t)$ from (65) (d), for all $\lambda_{n,u}^{st} \in S_{\lambda_{n,u}}(t)$ and for all $\lambda_{n,b}^{st}(t)$ in $S_z(t)$, respectively.

5. **Find** the sets $\bar{S}_{\lambda_{t,u}^{st}}(t)$ and $\bar{S}_{\lambda_{t,b}^{st}}(t)$ of multipliers $\lambda_{t,u}^{st}(t)$ and $\lambda_{t,b}^{st}(t)$, respectively, which belong to $S_z(t)$ and to $S_{\lambda_{t,u}^{st}}(t)$ and $S_{\lambda_{t,b}^{st}}(t)$, respectively.

6. If $\bar{S}_{\lambda_{t,u}^{st}}(t) = \emptyset$ or $\bar{S}_{\lambda_{t,b}^{st}}(t) = \emptyset$, **then** modify data $q(t)$, $F_{ext}$, $\mu_i$, and **go to** step 1.

7. **Else if** $\bar{S}_{\lambda_{t,u}^{st}}(t) \neq \emptyset$ and $\bar{S}_{\lambda_{t,b}^{st}}(t) \neq \emptyset$, **then**:

8. **Terminate**.

If we assume that conditions for existence of solutions hold, then the sets $S_z(t)$, $S_{\lambda_{n,u}^{st}}(t)$ are guaranteed to be both non empty. Some steps of this algorithm may be non trivially solved. Especially some sets may be unbounded. However one may restrict the search within contact forces which are
technologically feasible and admissible (hence adding non mathematical constraints stemming from the actuators technological features). Step 3 may be already solved in step 2, see for instance Example 6. An additional step that eliminates contacts *i* which detach ($a_{n,u,i}(t^+) > 0$) or such that the reaction force lies on the boundary of the friction cone, can be added for the sake of Lyapunov stability study [72] (different from the stability notions defined in [52]).

**Example 6** Let us illustrate the All-Sticking Bilateral/Unilateral Contact Algorithm on a simple example, as depicted in Figure 6. For this system we have: $q = (x y x_1, x_2)^\top$ ($y_1$ and $y_2$ are supposed to be constants), $n = 4$, $m_b = 1$, $m_u = 2$, $m = 3$, $d = 2$, $h_{n,b}(q) = x - a$, $h_{n,u,1}(q) = x - x_1 - L$, $h_{n,u,2}(q) = x_2 - x - L$, $F(q, \dot{q}, t) = (0 \ \mathbf{m}g \ F_1 \ F_2)^\top$, $M(q) = M = \text{diag}(\mathbf{m}, \mathbf{m}, \mathbf{m}_1, \mathbf{m}_2)$, $H_{t,u}^{st} =$

$$ \begin{pmatrix} 0 & 0 \\ 1 & 1 \\ 0 & 0 \\ 0 & 0 \end{pmatrix}, M_{u,b}^{st} = \begin{pmatrix} \mathbf{m} & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & \mathbf{m} & 0 & 0 & 0 & -1 & -1 \\ 0 & 0 & \mathbf{m}_1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \mathbf{m}_2 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}. \text{ To simplify the notation we will assume that} $$

all masses satisfy $\mathbf{m} = \mathbf{m}_1 = \mathbf{m}_2 = 1$ kg. Since $v_{t,1} = v_{t,2} = \dot{y}$, we assume that $a_{t,1}(t) = a_{t,2}(t) = \ddot{y}(t) = 0$. **Step 1:** $F_1(t) < 0$ and $F_2(t) > 0$, $\mu_1 \ge 0$, $\mu_2 \ge 0$. **Step 2:** using the system's dynamics and the bilateral constraint, it follows that (65) (b) is $0 \le (\lambda_{n,u,1}^{st} \ \lambda_{n,u,2}^{st}) \perp (\lambda_{n,u,1}^{st} \ \lambda_{n,u,2}^{st}) + (\begin{pmatrix} F_1 \\ -F_2 \end{pmatrix}) \ge 0$. Thus $\lambda_{n,u,1}^{st} = -F_1(t)$ and $\lambda_{n,u,2}^{st} = F_2(t)$, $\lambda_{n,b}^{st} = F_1(t) + F_2(t)$. It follows that $S_z(t) = \{z(t) \in \mathbb{R}^7 \mid z(t) = (0 \ 0 \ 0 \ 0 \ F_1 + F_2 \ \lambda_{t,u,1}^{st} \ \lambda_{t,u,2}^{st})^\top\}$. **Step 3:** done at step 2. **Step 4:** $S_{\lambda_{t,u}^{st}}(t) = \{(\lambda_{t,u,1}^{st} \ \lambda_{t,u,2}^{st}) \mid \lambda_{t,u,1}^{st} \in -F_1(t)\mathcal{D}_{\mu_1} = [\mu_1 F_1(t), -\mu_1 F_1(t)], \lambda_{t,u,2}^{st} \in -F_2(t)\mathcal{D}_{\mu_2} = [-\mu_2 F_2(t), \mu_2 F_2(t)]\}$. **Step 5:** one calculates that $S_{\lambda_{t,u}^{st}}(t) = \{(\lambda_{t,u,1}^{st} \ \lambda_{t,u,2}^{st}) \mid \lambda_{t,u,1}^{st} \in [\mu_1 F_1(t), -\mu_1 F_1(t)], \lambda_{t,u,2}^{st} \in [-\mu_2 F_2(t), \mu_2 F_2(t)], \text{ and } \lambda_{t,u,1} + \lambda_{t,u,2} = g\}$. **Step 6:** If $-\mu_1 F_1(t) + \mu_2 F_2(t) \ge g$, then terminate. Else go to step 1 and modify the data.

Notice that if $F_1(t) > 0$ and $F_2(t) < 0$ then the normal multipliers are both zero and conditions for detachment of both fingers hold. It is also noteworthy that the sticking mode persists in a right neighborhood of $t$, provided the data chosen at step 1, $F_1(t)$ and $F_2(t)$, are continuous in a right neighborhood of $t$. But, if these two external forces jump at $t$, the deactivation of the unilateral contacts may occur, or contact may be sustained but sliding occurs because the normal components are no longer large enough.

The system can be modified by allowing for varying $y_1(t)$ and $y_2(t)$, with 2-dimensional external forces on both fingers, so that $\ddot{y}(t) \ne 0$ (the mass has a vertical motion with sticking contacts at $A_1$ and $A_2$).

$$ \text{This time } q = (x \ y \ x_1 \ y_1 \ x_2 \ y_2)^\top \in \mathbb{R}^6, F(q, \dot{q}, t) = (0 \ g \ F_{1,x} \ F_{1,y} \ F_{2,x} \ F_{2,y})^\top, H_{t,u}^{st} = \begin{pmatrix} 0 & 0 \\ 1 & 1 \\ 0 & 0 \\ -1 & 0 \\ 0 & 0 \\ 0 & -1 \end{pmatrix}. $$

The constraints are the same as above, but this time $v_{t,1} = \dot{y} - \dot{y}_1$, $v_{t,2} = \dot{y} - \dot{y}_2$. Again we will assume that all masses satisfy $\mathbf{m} = \mathbf{m}_1 = \mathbf{m}_2 = 1$ kg. **Step 1:** $F_{1x}(t) < 0$, $F_{2x}(t) > 0$, $F_{1y}(t)$, $F_{2y}(t)$,
$\mu_1 \ge 0, \mu_2 \ge 0$. **Step 2:** similarly as above it is found that $\lambda_{n,u,1}^{st}(t) = -F_{1x}(t)$ and $\lambda_{n,u,2}^{st}(t) = F_{2x}(t)$, while $\lambda_{t,u,1}^{st}(t) = \lambda_{t,u,1}^{st}(t) = \frac{1}{3}(g - F_{1y}(t)), \ddot{y}(t) = \ddot{y}_1(t) = \ddot{y}_2(t) = -\frac{1}{3}(g + 2F_{1y}(t)).$ It follows that $S_z(t) = \{z(t) \in \mathbb{R}^8 | z(t) = (0 - \frac{1}{3}(g + 2F_{1y}(t)) \ 0 - \frac{1}{3}(g + 2F_{1y}(t)) \ 0 - \frac{1}{3}(g + 2F_{1y}(t)) \ \frac{1}{3}(g - F_{1y}(t)) \ \frac{1}{3}(g - F_{1y}(t)))^\top\}$. **Step 3:** done at step 2. **Step 4:** $S_{\lambda_{t,u}^{st}}(t) = \{(\lambda_{t,u,1}^{st}, \lambda_{t,u,2}^{st}) | \lambda_{t,u,1}^{st} \in [\mu_1 F_{1x}(t), -\mu_1 F_{1x}(t)], \lambda_{t,u,2}^{st} \in [-\mu_2 F_{2x}(t), \mu_2 F_{2x}(t)]\}$. **Step 5:** $\bar{S}_{\lambda_{t,u}^{st}}(t) = \{(\lambda_{t,u,1}^{st}, \lambda_{t,u,2}^{st}) | \lambda_{t,u,1}^{st} = \frac{1}{3}(g - F_{1y}(t)) \in [\mu_1 F_{1x}(t), -\mu_1 F_{1x}(t)], \lambda_{t,u,2}^{st} = \frac{1}{3}(g - F_{1y}(t)) \in [-\mu_2 F_{2x}(t), \mu_2 F_{2x}(t)]\}$. **Step 6:** If $|\frac{1}{3}(g - F_{1y}(t))| \le \min[-\mu_1 F_{1x}, \mu_2 F_{2x}]$, then terminate. Else go to step 1 and modify the data.

When **Step 6** fails, some intuition from the problem's geometry may be necessary to find suitable values of the data. Notice that Assumption 1 is satisfied in this example, so that Proposition 9 could be applied also.

<figure><img src="image_9.png" /><figcaption>Figure 6: A simple hard-finger manipulation robotic task.</figcaption></figure>

**Remark 4** As alluded to in the introduction of the article, constant singular mass matrices often occur when the natural coordinates are used [37], and are necessarily accompanied by nonlinear bilateral constraints which do not stem from mechanical contacts or joints. It is therefore of interest to consider the case of singular $M(q) \ge 0$, unilateral constraints with friction, and frictionless bilateral constraints, all constraints being potentially redundant. The frictionless case is easily obtained by setting vanishing friction coefficients where needed in (65) (d).

**Example 7** Let us outline a manipulation problem with hard-finger contact points with friction, as depicted in Figure 7, which can be analysed through (65) using an extension of Proposition 8 as indicated above. This example encompasses Example 5. The mass matrix of this system is $M(q) =$
blockdiag($M_{ob}(q_{ob})$, $M_1(q_1)$, ..., $M_m(q_m)$), $q_{ob} \in \mathbb{R}^{n_{ob}}$, $q_i \in \mathbb{R}^{n_i}$, $n = n_{ob} + \sum_{i=1}^{m} n_i$, $m = m_u$, $A_i$ are the contact points between the "fingers" and the object, $1 \le i \le m$. In full 3D generality, one may have $n_{ob} = 6$ but also $n_{ob} > 6$ if natural coordinates are used<fnref n="6" />. Then $M_{ob}(q_{ob}) \ge 0$ and frictionless bilateral constraints $h_{ob}(q_{ob}) = 0$ have to be added to the problem. The dynamics of the object is $M_{ob}(q_{ob})\ddot{q}_{ob} + F_{ob}(q_{ob}, \dot{q}_{ob}, t) = \sum_{i=1}^{m} M_{i/ob}^{n}(q_{ob}, q_i)\lambda_{n,u,i}^{st} + \nabla h_{ob}(q_{ob})\lambda_{ob} + \sum_{i=1}^{m} M_{i/ob}^{t}(q_{ob}, q_i)\lambda_{t,u,i}$. The dynamics of each arm $i$ is $M_i(q_i)\ddot{q}_i + F_i(q_i, \dot{q}_i, t) = M_{ob/i}^{n}(q_{ob}, q_i)\lambda_{n,u,i}^{st} + M_{ob/i}^{t}(q_{ob}, q_i)\lambda_{t,u,i}^{st}$. Recall that $\lambda_{n,u,i}^{st} = F_{n,i} \ge 0$ where $F_{n,i}$ is the force normal component in each local frame at contact $i$, while $\lambda_{t,u,i} \in \mathbb{R}^{d-1}$ is the local tangential force at contact $i$. Recall also that if $v_{r,t,i}$ denotes the relative tangential velocity at contact $i$, then $\sum_{i=1}^{m} v_{r,t,i}\lambda_{t,u,i} = v_{r,t}^{\top}\lambda_{t,u} = \dot{q}^{\top}H_{t,u}(q_{ob}, q_1, \dots, q_m)\lambda_{t,u} = \dot{q}^{\top}\Lambda_{t,u}$ from the virtual work principle. Thus:

$$ \nabla h_{n,u}^{st}(q_{ob}, q_1, \dots, q_m) = \begin{pmatrix} M_{1/ob}^{n} & M_{2/ob}^{n} & \dots & \dots & M_{m/ob}^{n} \\ M_{ob/1}^{n} & 0 & \dots & \dots & 0 \\ 0 & M_{ob/2}^{n} & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & \dots & \dots & 0 & M_{ob/m}^{n} \end{pmatrix} \in \mathbb{R}^{(m+1) \times m}, \quad (67) $$

and

$$ H_{t,u}^{st}(q_{ob}, q_1, \dots, q_m) = \begin{pmatrix} M_{1/ob}^{t} & M_{2/ob}^{t} & \dots & \dots & M_{m/ob}^{t} \\ M_{ob/1}^{t} & 0 & \dots & \dots & 0 \\ 0 & M_{ob/2}^{t} & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & \dots & \dots & 0 & M_{ob/m}^{t} \end{pmatrix} \in \mathbb{R}^{(m+1) \times (d-1)m}, \quad (68) $$

for some matrices $M_{i/ob}^{t}$, $M_{ob/i}^{t}$, $M_{i/ob}^{n}$, $M_{ob/i}^{n}$, which depend on the object's geometry, $1 \le i \le m$. The matrix $M_{u,b}^{st}(q)$ in (65) (a) therefore possesses a specific, sparse structure, which should be taken into account to design efficient solvers for (66), as alluded to in Example 5. It is clear however, that excepted in very simple cases (like manipulation of spherical objects with prismatic hard fingers), the gap functions and their gradients are nonlinear functions of $q$ and so is $M_{u,b}^{st}(q)$. This may in general complicate step 6 of the above Contact Algorithm.

**Remark 5 (Stick → slip transitions)** In section 5 the case of both sticking and sliding contacts will be analyzed. Obviously the intermediate case represented by stick → slip transitions is of interest. Active contacts with $v_{t,i}(t) = 0$ which can potentially undergo such transition, must have $\lambda_{t,u,i}^{st}(t) \in \lambda_{n,u,i}^{st}(t)\text{bd}(D_{\mu_i})$ and $\lambda_{t,b,i}^{st}(t) \in \lambda_{n,b,i}^{st}(t)\text{bd}(D_{\mu_i})$, where bd denotes the boundary. Coulomb's law at the acceleration in (5) and its sticking mode formulation in (6) or (7), is useful, because it allows one to calculate $a_{t,i}(t)$ without assuming a priori that $a_{t,i}(t) = 0$ as in (65). Consequently the obtained

<footnote n="6">⁶Which is the case in many multibody commercial software packages.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 7: A hard-finger manipulation robotic task.</figcaption></figure>

problem is different from (65). It gives rise to the following problem:

$$
\left\{
\begin{array}{l}
(a) \quad \underbrace{\begin{pmatrix} M(q) & -\nabla h_{n,b}^{st}(q) \\ \nabla h_{n,b}^{st}(q) & 0 \end{pmatrix}}_{\triangleq M_b^{st/sl}(q)} \begin{pmatrix} \ddot{q}(t) \\ \lambda_{n,b}^{st}(t) \end{pmatrix} + \begin{pmatrix} F(q, \dot{q}, t) \\ \frac{d}{dt}((\nabla h_{n,b}^{st})^{\top})\dot{q} \end{pmatrix} \\
= \begin{pmatrix} \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} + H_{t,u}^{st}(q)\lambda_{t,u}^{st} + H_{t,b}^{st}(q)\lambda_{t,b}^{st} \\ 0 \end{pmatrix} \\
(b) \quad \lambda_{n,u}^{st} \in -\mathcal{N}_{\mathbb{R}_+^{m_u}}(\nabla h_{n,u}^{st}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}((\nabla h_{n,u}^{st})^{\top})\dot{q}) \\
(c) \quad \Lambda_{t,b,i}^{st} \in \partial\psi_{\mathcal{D}(\lambda_{n,b,i}^{st}, \mu_i)}^{\star}(H_{t,b,i}^{st}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}(H_{t,b,i}^{st}(q)^{\top})\dot{q}), \quad 1 \le i \le m_b, \\
(d) \quad \Lambda_{t,u,i}^{st} \in \partial\psi_{\mathcal{D}(\lambda_{n,u,i}^{st}, \mu_i)}^{\star}(H_{t,u,i}^{st}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}(H_{t,u,i}^{st}(q)^{\top})\dot{q}), \quad m_b+1 \le i \le m.
\end{array}
\right.
\tag{69}
$$

The variables $\ddot{q}(t)$ and $\lambda_{n,b}^{st}(t)$ are the primary unknowns of the generalized equation in (69) (the other multipliers are “eliminated” inserting the inclusions (b) (c) (d) into the right-hand side of (a)). The major difficulty in analyzing and solving (69) lies in the fact that the convex sets $\mathcal{D}(\lambda_{n,u,i}^{st}, \mu_i)$ and $\mathcal{D}(\lambda_{n,b,i}^{st}, \mu_i)$ depend on the unknowns $\lambda_{n,u,i}^{st}$ and $\lambda_{n,b,i}^{st}$. It is noteworthy that due to the virtual work principle, all three terms in the right-hand side of (69) possess a structure which allows one to apply the chain rule of Convex Analysis, provided some constraints qualifications hold. If the normal components of the reaction force are known, this allows one to rewrite (69) in the form of a generalized equation with unknowns $\ddot{q}(t)$ and $\lambda_{n,b}^{st}(t)$, involving maximal monotone set-valued terms as follows:

$$ \begin{pmatrix} M(q) & -\nabla h_{n,b}^{st}(q) \\ \nabla h_{n,b}^{st}(q) & 0 \end{pmatrix} \begin{pmatrix} \ddot{q}(t) \\ \lambda_{n,b}^{st}(t) \end{pmatrix} + \begin{pmatrix} F(q, \dot{q}, t) \\ \frac{d}{dt}((\nabla h_{n,b}^{st})^{\top})\dot{q} \end{pmatrix} \in \left( -\mathcal{N}_{\Phi_{n,u}^{st}(q,\dot{q})} (\ddot{q}(t)) + \partial f_{q,\dot{q}}^{t,u}(\ddot{q}(t)) + \partial f_{q,\dot{q}}^{t,b}(\ddot{q}(t)) \right) \quad (70) $$

where $\Phi_{n,u}^{st}(q, \dot{q}) = \{w \in \mathbb{R}^n \mid \nabla h_{n,u}^{st}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}((\nabla h_{n,u}^{st})^{\top})\dot{q} \in \mathbb{R}_+^{m_u}\}$, $f_{q,\dot{q},i}^{t,u}(\ddot{q}(t)) = \psi_{\mathcal{D}(\lambda_{n,u,i}^{st}, \mu_i)}^{\star} \circ (H_{t,u,i}^{st}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}(H_{t,u,i}^{st}(q)^{\top})\dot{q})$ and $f_{q,\dot{q},i}^{t,b}(\ddot{q}(t)) = \psi_{\mathcal{D}(\lambda_{n,b,i}^{st}, \mu_i)}^{\star} \circ (H_{t,b,i}^{st}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}(H_{t,b,i}^{st}(q)^{\top})\dot{q})$ (using
[35, Theorem 4.2.1]). When the normal components are known, this generalized equation could be analyzed using the material in section B. The kernel of $M_b^{st/sl}(q)$ is once again central in the analysis, similarly as in the previous results (see Propositions 1 and 8). We do not analyze further the stick/slip problem, see [53] for a complete study (with unilateral constraints only).

# 5 Systems with 3D sticking and sliding unilateral contacts

Let us outline the problem which consists of 3-dimensional contacts, some of which are sticking, the other ones are sliding. To simplify the presentation we restrict ourselves to unilateral contacts, i.e., $\mathcal{I}_b = \emptyset$, $\mathcal{I}_u = \mathcal{I}_u^{\mu,st} \cup \mathcal{I}_u^{\mu,sl}$, and $\mathcal{I}_u^0 = \emptyset$. The contact problem in its index-1 formulation is:

$$ egin{cases} 
(a) M(q)\ddot{q}(t) + F(q, \dot{q}, t) = \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} + H_{t,u}^{st}(q)\lambda_{t,u}^{st} + \nabla h_{n,u}^{sl}(q)\lambda_{n,u}^{sl} + H_{t,u}^{sl}(q)\lambda_{t,u}^{sl} \\ 
(b) H_{t,u}^{st}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q} = 0 \\ 
(c) 0 \le \begin{pmatrix} \lambda_{n,u}^{st} \\ \lambda_{n,u}^{sl} \end{pmatrix} \perp \begin{pmatrix} \nabla h_{n,u}^{st}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \\ \nabla h_{n,u}^{sl}(q)^{\top}\ddot{q}(t) + \frac{d}{dt}(\nabla h_{n,u}^{sl}(q)^{\top})\dot{q} \end{pmatrix} \ge 0 \\ 
(d) \Lambda_{t,u,i}^{st} \in \mathcal{D}(\mu_i, \lambda_{n,u,i}^{st}), i \in \mathcal{I}_u^{\mu,st}. 
\end{cases} 
\quad (71) $$

Let us write $\lambda_{t,u}^{sl} = -[\mu_u^{sl}][\operatorname{sgn}(v_{t,u}^{sl})]\lambda_{n,u}^{sl}$, where the equality holds because sliding is assumed for these contact points, and for $x \in \mathbb{R}^n$, $[x] \triangleq \operatorname{diag}(x_i)$. The problem in (71) (a) (b) (c) is equivalent to the generalized equation (time argument is dropped):

$$ egin{gathered} 
\begin{pmatrix} M(q) & -H_{t,u}^{st}(q) \\ H_{t,u}^{st}(q)^{\top} & 0 \end{pmatrix} \begin{pmatrix} \ddot{q} \\ \lambda_{t,u}^{st} \end{pmatrix} + \begin{pmatrix} F(q, \dot{q}, t) \\ 0 \end{pmatrix} = \begin{pmatrix} (\nabla h_{n,u}^{st}(q) \quad \nabla h_{n,u}^{sl}(q) - H_{t,u}^{sl}(q)[\mu_u^{sl}][\operatorname{sgn}(v_{t,u}^{sl})]) & \begin{pmatrix} \lambda_{n,u}^{st} \\ \lambda_{n,u}^{sl} \end{pmatrix} \\ 0 & 0 \end{pmatrix} \\ 
\begin{pmatrix} \lambda_{n,u}^{st} \\ \lambda_{n,u}^{sl} \end{pmatrix} \in -\mathcal{N}_{\mathbb{R}_+^{mu}} \left( \begin{pmatrix} (\nabla h_{n,u}^{st}(q)^{\top}) \ddot{q} \\ (\nabla h_{n,u}^{sl}(q)^{\top}) \ddot{q} \end{pmatrix} + \begin{pmatrix} \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \\ \frac{d}{dt}(\nabla h_{n,u}^{sl}(q)^{\top})\dot{q} \end{pmatrix} \right), 
\end{gathered} 
\quad (72) $$

which is obviously more complex than previous problems, due to sliding friction which introduces a distortion in the gradient in the right-hand side of the first equality in (72). This problem is close to the frictional bilaterally/unilaterally constrained contact problem studied in [15, section 3.4, Proposition 13], and could be analysed in a similar way. Due to space limitations this is not tackled here. It is interesting to observe how the basic KKT problem gets more complex when starting from (26), then (43), (61), (65), (69) and finally (72).

The problem may be simplified as follows. From (71), assuming $M(q) > 0$ and disregarding for the moment the friction disk inclusion in (71) (d), one eliminates $\ddot{q}(t)$ using (71) (a) and constructs the
following MLCP with unknowns $\lambda_{n,u}^{st}$, $\lambda_{t,u}^{st}$, $\lambda_{n,u}^{sl}$:

$$ (a) \quad W_{t,u}^{st}(q)\lambda_{t,u}^{st} + H_{t,u}^{st}(q)M(q)^{-1}\nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} + A_{t,u}^{st,sl}(q)\lambda_{n,u}^{sl} - F_{t,u}^{st}(q, \dot{q}, t) = 0 $$

$$ (b) \quad 0 \le \begin{pmatrix} \lambda_{n,u}^{st} \\ \lambda_{n,u}^{sl} \\ \lambda_{t,u}^{st} \end{pmatrix} \perp \begin{pmatrix} W_{n,u}^{st}(q) & A_{n,u}^{st,sl}(q) & \nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q) \\ \nabla h_{n,u}^{sl}(q)^{\top} M(q)^{-1} \nabla h_{n,u}^{st}(q) & A_{n,u}^{sl,sl}(q) & \nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q) \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \lambda_{n,u}^{st} \\ \lambda_{n,u}^{sl} \\ \lambda_{t,u}^{st} \end{pmatrix} \\ - \begin{pmatrix} F_{n,u}^{st}(q, \dot{q}, t) \\ F_{n,u}^{sl}(q, \dot{q}, t) \\ 0 \end{pmatrix} \ge 0, \tag{73} $$

with:

$$ \begin{align} F_{t,u}^{st}(q, \dot{q}, t) &= H_{t,u}^{st}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) - \frac{d}{dt} (H_{t,u}^{st}(q)^{\top}) \dot{q} \nonumber \\ F_{n,u}^{sl}(q, \dot{q}, t) &= \nabla h_{n,u}^{sl}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) - \frac{d}{dt} (\nabla h_{n,u}^{sl}(q)^{\top}) \dot{q} \nonumber \\ F_{n,u}^{st}(q, \dot{q}, t) &= \nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} F(q, \dot{q}, t) - \frac{d}{dt} (\nabla h_{n,u}^{sl}(q)^{\top}) \dot{q} \nonumber \\ A_{t,u}^{st,sl}(q) &= H_{t,u}^{st}(q)^{\top} M(q)^{-1} (\nabla h_{n,u}^{sl}(q) - H_{t,u}^{sl}(q)[\mu_u^{sl}]\text{sgn}(v_{t,u}^{sl})) \tag{74} \\ W_{t,u}^{st}(q) &= H_{t,u}^{st}(q)^{\top} M(q)^{-1} H_{t,u}^{st}(q) \ge 0 \nonumber \\ W_{n,u}^{st}(q) &= \nabla h_{n,u}^{st}(q) M(q)^{-1} \nabla h_{n,u}^{st}(q) \ge 0 \nonumber \\ A_{n,u}^{st,sl}(q) &= \nabla h_{n,u}^{st}(q)^{\top} M(q)^{-1} (\nabla h_{n,u}^{sl}(q) - H_{t,u}^{sl}(q)[\mu_u^{sl}]\text{sgn}(v_{t,u}^{sl})) \nonumber \\ A_{n,u}^{sl,sl}(q) &= \nabla h_{n,u}^{sl}(q)^{\top} M(q)^{-1} (\nabla h_{n,u}^{sl}(q) - H_{t,u}^{sl}(q)[\mu_u^{sl}]\text{sgn}(v_{t,u}^{sl})). \nonumber \end{align} $$

Let us state the following propositions.

**Proposition 13** Let $M(q) \succ 0$ and $\frac{d}{dt}(H_{t,u}^{st}(q)^{\top})\dot{q} \in \ker^{\perp}(H_{t,u}^{st}(q))$. Then the MLCP in (71) (a) (b) (c) with the Coulomb's disk constraint reduces to the problem of existence of $\lambda_{n,u}^{st}$, $\lambda_{n,u}^{sl}$ and $\lambda_{t,u}^{st}$ such that:

$$ (a) \quad 0 \le \begin{pmatrix} \lambda_{n,u}^{st} \\ \lambda_{n,u}^{sl} \end{pmatrix} \perp \begin{pmatrix} \nabla h_{n,u}^{st}(q)^{\top} (z+y) + \frac{d}{dt} (\nabla h_{n,u}^{st}(q)^{\top}) \dot{q} \\ \nabla h_{n,u}^{sl}(q)^{\top} (z+y) + \frac{d}{dt} (\nabla h_{n,u}^{sl}(q)^{\top}) \dot{q} \end{pmatrix} \ge 0 $$

$$ (b) \quad \lambda_{t,u}^{st} \in \mathcal{D}(\mu_u^{st}, \lambda_{n,u}^{st}) $$

$$ (c) \quad H_{t,u}^{st}(q)\lambda_{t,u}^{st} + \nabla h_{n,u}^{st}(q)\lambda_{n,u}^{st} + (\nabla h_{n,u}^{sl}(q) - H_{t,u}^{sl}(q)[\mu_u^{sl}]\text{sgn}(v_{t,u}^{sl}))\lambda_{n,u}^{sl} - F(q, \dot{q}, t) + M(q)(z + y) = 0 \tag{75} $$
for some $z \in \mathbb{R}^n$ and $y \in \text{Im}^{\perp}(H_{t,u}^{st}(q))$.

**Proof 13** The proof follows from [15, Proposition 8], noting the similarity between (71) (a) (b) (c) and [15, equation (16)]. ■
The degenerated complementarity problem in (75) (a) has a solution if and only if

$$ \left( \begin{array}{c} \nabla h_{n,u}^{st}(q)^{\top}[z+y] + \frac{d}{dt}(\nabla h_{n,u}^{st}(q)^{\top})\dot{q} \\ \nabla h_{n,u}^{sl}(q)^{\top}[z+y] + \frac{d}{dt}(\nabla h_{n,u}^{sl}(q)^{\top})\dot{q} \end{array} \right) \geq 0 \qquad (76) $$

for some $z \in \mathbb{R}^n$ and $y \in \text{Im}^{\perp}(H_{t,u}^{st}(q))$. This can be used as a test to eliminate systems which do not verify such necessary condition. Using the proposition's assumption and (75) (b) (c), one finds that $\lambda_{n,u}^{st}$ and $\lambda_{n,u}^{sl}$ have to satisfy:

$$ \begin{gathered} (\nabla h_{n,u}^{st}(q) \ \nabla h_{n,u}^{sl}(q) - H_{t,u}^{sl}(q)[\mu_u^{sl}][\text{sgn}(v_{t,u}^{sl})]) \begin{pmatrix} \lambda_{n,u}^{st} \\ \lambda_{n,u}^{sl} \end{pmatrix} - F(q, \dot{q}, t) - M(q)[z+y] \qquad (77) \\ \in -H_{t,u}^{st}(q)\mathcal{D}(\mu_u^{st}, \lambda_{n,u}^{st}) \end{gathered} $$

for the same $z$ and $y$. Intuitively, small $\mu_i^{sl}$ and large $\mu_i^{st}$ could be needed for (77) to be solvable, which complies with the fact that sliding friction usually has to be small enough to guarantee existence of solutions [15] [19, Section 5.6]. Let us now focus on the MLCP in (73). The first thing to look at is the well-posedness of the LCP in (73) (b). Various sufficient conditions may be stated to guarantee that this LCP possesses solutions. Conditions may be imposed that guarantee that the LCP matrix is positive semi-definite (that is, it is also copositive), using Corollary 5. Then Theorem 1 in Appendix C can be applied. Roughly speaking, this will imply that sliding contacts have small enough friction coefficients so that $A_{n,u}^{sl,sl}(q) \succcurlyeq 0$, while friction coefficients at sticking contact be large enough so that (71) (d) is verified.

**Remark 6** The above problems in (65) and (71) (or (73)), are related to the so-called force and form closure in dexterous manipulation and grasp analysis [14, 62, 10, 33, 68]. As Example 6 shows, the All-Sticking Bilateral/Unilateral Contact Algorithm could be used to determine if force or form closure hold. The criteria that are proposed all along this article, supply general ways to numerically solve the all-sticking or mixed sticking/sliding (Proposition 13) problems, see the algorithm in section 4.2. It seems difficult to get more precise results in terms of simple-to-use analytical criteria using the above generic formulations, without specifying further the data.

**Remark 7** Several results and criteria presented in this article, use full-rank assumptions and matrix inversion. It is noteworthy that the same problem (consider for instance (14) or (43) with $M(q) \succcurlyeq 0$ or $M(q) \succcurlyeq 0$) can be solved numerically either after a transformation using matrix inversion (see the LCP in (37)), or solving directly the MLCP in (43) (or its bilateral counterpart in (14)) with suitable methods like interior point algorithms. Nevertheless matrix inversion can be used if dimensions are low, or just for the sake of analysis (as in Propositions 6 and 7, Corollary 2, Corollary 4, Lemma 2). Whether or not the proposed criteria could be used in on-line implementations for checking stability or form closure, is yet unanswered to.
# 6 Conclusion

This article analyses multibody systems undergoing possibly redundant bilateral and/or unilateral contacts with Coulomb's set-valued friction, with a possibly singular mass matrix and redundant constraints, when the contacts are sticking in both the normal and tangential directions. Various results are given, which state under which conditions the problem has at least one solution. The cases with pure bilateral, pure unilateral, and mixed bilateral/unilateral constraints, are treated, as well as mixed sticking/sliding contacts. Complementarity theory (in particular mixed linear complementarity problems) as well as variational inequalities, are used for the analysis. The results provide criteria which should, in principle, be verifiable numerically. An algorithm is proposed which paves the way to force and form closure verification, for object manipulation. Several examples illustrate the theoretical findings. One major future task is to develop complete, non-academic examples, and efficient numerical solvers.

# A Useful results from Linear Algebra

The first part of this lemma is taken from [59], the second part is [12, Fact 6.4.29], see also [67, equation (1.5)].

**Lemma 3** Let $M = \begin{pmatrix} A & C \\ C^\top & B \end{pmatrix}$ be symmetric and positive semidefinite. Assume that $Q = B - C^\top A^\dagger C$ is nonsingular. Then the Moore-Penrose generalized inverse of $M$ is given by:

$$ M^\dagger = \begin{pmatrix} A^\dagger + A^\dagger C Q^\dagger C^\top A^\dagger & -A^\dagger C Q^\dagger \\ -Q^\dagger C^\top A^\dagger & Q^\dagger \end{pmatrix}. \qquad (78) $$

When $B = 0$ one obtains with $E = A + CC^\top$ and $D = C^\top E^\dagger C$:

$$ M^\dagger = \begin{pmatrix} E^\dagger - E^\dagger C D^\dagger C^\top E^\dagger & E^\dagger C D^\dagger \\ (E^\dagger C D^\dagger)^\top & D D^\dagger - D^\dagger \end{pmatrix}. \qquad (79) $$

**Proposition 14** [12, Proposition 6.1.7] Let $A \in \mathbb{R}^{n \times m}$ and $b \in \mathbb{R}^n$. Then the two statements are equivalent:

* (i) There exists a vector $x \in \mathbb{R}^m$ satisfying $Ax = b$.

* (ii) $AA^\dagger b = b$.

If (i) or (ii) is satisfied, then for all $y \in \mathbb{R}^m$, $x = A^\dagger b + (I - A^\dagger A)y$ satisfies $Ax = b$, and $y = 0$ minimizes $x^\top x$.
## B Well-posedness of Variational Inequalities

The next results use the notions of recession functions and cones, which we briefly introduce now (see [20, 32] for illustrating examples). Let $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be a proper convex and lower semi-continuous function, we denote by $\text{dom}(f) \stackrel{\Delta}{=} \{x \in \mathbb{R}^n | f(x) < +\infty\}$ the domain of the function $f(\cdot)$. The epigraph of $f(\cdot)$ is the set $\text{epi}(f) \stackrel{\Delta}{=} \{(x, \alpha) \in \mathbb{R}^n \times \mathbb{R} | \alpha \ge f(x)\}$. The Fenchel transform $f^*(\cdot)$ of $f(\cdot)$ is the proper, convex and lower semi-continuous function defined by

$$ (\forall z \in \mathbb{R}^n) |f^*(z)| = \sup_{x \in \text{dom}(f)} \{\langle x, z \rangle - f(x)\}. \qquad (80) $$

The subdifferential $\partial f(x)$ of $f(\cdot)$ at $x \in \mathbb{R}^n$ is defined by

$$ \partial f(x) = \{\omega \in \mathbb{R}^n | f(v) - f(x) \ge \langle \omega, v - x \rangle, \forall v \in \mathbb{R}^n\}. $$

We denote by $\text{Dom}(\partial f) \stackrel{\Delta}{=} \{x \in \mathbb{R}^n | \partial f(x) \ne \emptyset\}$ the domain of the subdifferential operator $\partial f: \mathbb{R}^n \to \mathbb{R}^n$. Recall that (see e.g. Theorem 2, Chapter 10, Section 3 in [8]): $\text{Dom}(\partial f) \subset \text{dom}(f)$.

Let $x_0$ be any element in the domain $\text{dom}(f)$ of $f(\cdot)$, the recession function $f_\infty(\cdot)$ of $f(\cdot)$ is defined by

$$ (\forall x \in \mathbb{R}^n) : f_\infty(x) = \lim_{\lambda \to +\infty} \frac{1}{\lambda} f(x_0 + \lambda x). $$

The function $f_\infty: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ is a proper convex and lower semi-continuous function which describes the asymptotic behavior of $f(\cdot)$.

Let $K \subset \mathbb{R}^n$ be a nonempty closed convex set. Let $x_0$ be any element in $K$. The recession cone of $K$ is defined by

$$ K_\infty = \bigcap_{\lambda>0} \frac{1}{\lambda} (K - x_0) = \{u \in \mathbb{R}^n | x + \lambda u \in K, \forall \lambda \ge 0, \forall x \in K\}. $$

The set $K_\infty$ is a nonempty closed convex cone that is described in terms of the directions which recede from $K$. The indicator function of a set $K \subset \mathbb{R}^n$ is $\psi_K(x) = 0$ if $x \in K$, $\psi_K(x) = +\infty$ if $x \notin K$. If $K$ is closed non empty convex, we have $\partial \psi_K(x) = N_K(x)$, the so-called normal cone to $K$ at $x$, defined as $N_K(x) = \{v \in \mathbb{R}^n | v^\top(s-x) \le 0 \text{ for all } s \in K\}$. When $K$ is finitely represented, i.e., $K = \{x \in \mathbb{R}^n | k_i(x) \ge 0, 1 \le i \le m\}$, and if the functions $k_i(\cdot)$ satisfy some constraint qualification (like, independency, or extensions like the MFCQ, see Appendix C), then $N_K(x)$ is generated by the outwards normals at the active constraints $k_i(x) = 0$, i.e., $N_K(x) = \{v \in \mathbb{R}^n | v = -\lambda_i \nabla k_i(x), k_i(x) = 0, \lambda_i \ge 0\}$.

Let us here recall some important properties of the recession function and recession cone (see e.g. [13, Proposition 1.4.8]):

**Proposition 15** The following statements hold:

a) Let $f_1: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ and $f_2: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be two proper, convex and lower semi-continuous functions. Suppose that $f_1 + f_2$ is proper. Then for all $x \in \mathbb{R}^n$: $(f_1 + f_2)_\infty(x) = (f_1)_\infty(x) + (f_2)_\infty(x)$.
b) Let $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be a proper, convex and lower semi-continuous function and let $K$ be a nonempty closed convex set, such that $f + \Psi_K$ is proper (equivalently $\text{dom}(f) \cap K$ is non empty). Then for all $x \in \mathbb{R}^n$: $(f + \Psi_K)_\infty(x) = f_\infty(x) + (\Psi_K)_\infty(x)$.

c) Let $K \subset \mathbb{R}^n$ be a nonempty, closed and convex set. Then for all $x \in \mathbb{R}^n$: $(\Psi_K)_\infty(x) = \Psi_{K_\infty}(x)$. Moreover for all $x \in K$ and $e \in K_\infty$: $x + e \in K$.

d) If $K \subset \mathbb{R}^n$ is a nonempty closed and convex cone, then $K_\infty = K$.

e) Let $K = P(a, b) \triangleq \{x \in \mathbb{R}^n | Ax \ge b\}$ for $A \in \mathbb{R}^{m \times n}$ and $b \in \mathbb{R}^m$. If $K \ne \emptyset$ then $K_\infty = P(A, 0) = \{x \in \mathbb{R}^n | Ax \ge 0\}$.

f) $K \subset \mathbb{R}^n$ is a non-empty closed convex bounded set if and only if $K_\infty = \{0_n\}$.

Let us now concatenate [3, Theorem 3, Corollaries 3 and 4]. They concern variational inequalities of the form: Find $u \in \mathbb{R}^n$ such that

$$ \langle \mathbf{M}u + \mathbf{q}, v - u \rangle + \varphi(v) - \varphi(u) \ge 0, \forall v \in \mathbb{R}^n \quad (81) $$

where $\mathbf{M} \in \mathbb{R}^{n \times n}$ is a real matrix, $\mathbf{q} \in \mathbb{R}^n$ a vector and $\varphi: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ a proper convex and lower semicontinuous function. The analogy with the generalized equation in (46) is clear taking $\varphi(\cdot) = \Psi_{\tilde{K}}(\cdot)$, thus restricting the variation of $v$ to $\tilde{K}$, $\mathbf{q} = F(q, \dot{q}, t)$ and $\mathbf{M} = M(q)$.

The problem in (81) is denoted as $VI(\mathbf{M}, \mathbf{q}, \varphi)$ in the next proposition. We also set:

$$ \mathcal{K}(\mathbf{M}, \varphi) = \{x \in \mathbb{R}^n | \mathbf{M}x \in (\text{dom}(\varphi_\infty))^*\}. \quad (82) $$

Note that $(\text{dom}(\varphi_\infty))^*$ is the dual cone of the domain of the recession function $\varphi_\infty$ while $(\text{dom}(\varphi))^*$ is the recession cone of $\text{dom}(\varphi)$.

**Proposition 16** [3] Let $\varphi: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be a proper, convex and lower semicontinuous function with closed domain, and suppose that $\mathbf{M} \in \mathbb{R}^{n \times n}$ is positive semi-definite (not necessarily symmetric).

a) If $(\text{dom}(\varphi))^* \cap \ker\{\mathbf{M} + \mathbf{M}^\top\} \cap \mathcal{K}(\mathbf{M}, \varphi) = \{0\}$ then for each $\mathbf{q} \in \mathbb{R}^n$, problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$ has at least one solution.

b) Suppose that $(\text{dom}(\varphi))^* \cap \ker\{\mathbf{M} + \mathbf{M}^\top\} \cap \mathcal{K}(\mathbf{M}, \varphi) \ne \{0\}$. If there exists $x_0 \in \text{dom}(\varphi)$ such that

$$ \langle \mathbf{q} - \mathbf{M}^\top x_0, v \rangle + \varphi_\infty(v) > 0, \forall v \in \text{dom}(\varphi)_\infty \cap \ker\{\mathbf{M} + \mathbf{M}^\top\} \cap \mathcal{K}(\mathbf{M}, \varphi), v \ne 0, \quad (83) $$

then problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$ has at least one solution.

b') If $\mathbf{M} = \mathbf{M}^\top$ then one can take $x_0 = 0$ in b).

c) If $u_1$ and $u_2$ denote two solutions of problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$ then $u_1 - u_2 \in \ker\{\mathbf{M} + \mathbf{M}^\top\}$.

d) If $\mathbf{M} = \mathbf{M}^\top$ and $u_1$ and $u_2$ denote two solutions of problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$, then $\langle \mathbf{q}, u_1 - u_2 \rangle = \varphi(u_2) - \varphi(u_1)$.

e) If $\mathbf{M} = \mathbf{M}^\top$ and $\varphi(x + z) = \varphi(x)$ for all $x \in \text{dom}(\varphi)$ and $z \in \ker\{\mathbf{M}\}$ and $\langle \mathbf{q}, e \rangle \ne 0$ for all $e \in \ker\{\mathbf{M}\}$, $e \ne 0$, then problem $VI(\mathbf{M}, \mathbf{q}, \varphi)$ has at most one solution.
Notice that the function $\varphi(\cdot)$ will never be strictly convex in our case (it is an indicator function) so that the strict convexity argument of [3, Theorem 5] which applies when $\mathbf{M}$ is a $P_0$-matrix never holds.

## C Some Convex Analysis and Complementarity Theory tools

**Theorem 1** [25, Theorem 3.8.6] Let $M \in \mathbb{R}^{n \times n}$ be copositive and let $q \in \mathbb{R}^n$ be given. If the implication: $0 \le v \perp Mv \ge 0 \Rightarrow v^\top q \ge 0$ is valid, then the LCP($M, q$) is solvable.

Let $Q_M$ denote the solution set of the homogeneous LCP. This theorem can be restated equivalently as: $v \in Q_M \Rightarrow q \in Q_M^*$.

The next corollary is proved in [18], and is a consequence of results in [23].

**Corollary 5** Let $D = P+N$, where $D, P$ and $N$ are $n \times n$ real matrices, and $P > 0$, not necessarily symmetric. If

$$ ||N||_2 < \frac{1}{||\left(\frac{P+P^\top}{2}\right)^{-1}||_2} \quad (84) $$

then $D > 0$.

If $K \subset \mathbb{R}^n$ is a set, then $K^* = \{z \in \mathbb{R}^n | \langle z, x \rangle \ge 0 \text{ for all } x \in K\}$ is its dual cone. Let $K$ be a closed convex cone, then:

$$ K^* \ni x \perp y \in K \Leftrightarrow x \in -\mathcal{N}_K(y) \Leftrightarrow y \in -\mathcal{N}_{K^*}(x). \quad (85) $$

Let $M = M^\top > 0$, $x$ and $y$ two vectors, then

$$ M(x - y) \in -\mathcal{N}_K(x) \Leftrightarrow x = \operatorname{proj}_M[K; y] \Leftrightarrow x = \min_{z \in K} \frac{1}{2}(z - y)^\top M(z - y). \quad (86) $$

We note that this is a particular case of (81), so that Proposition 16 can be considered as the characterization of a generalized projection operator $VI(\mathbf{M}, \mathbf{q}, \varphi)$.

In this paper we deal with non convex sets, for which it is needed to define suitable notions of normal and tangent cones. The so-called Mangasarian-Fromovitz constraint qualification (MFCQ) [26, p.17, 252] is also used.

**Definition 1 (The MFCQ)** Let $K$ be a finitely represented set, i.e. $K = \{x \in \mathbb{R}^n | h_i(x) \ge 0, 1 \le i \le m\}$, and let $\mathcal{J}(x) = \{i \in \{1, m\} | h_i(x) = 0\}$ be the set of active constraints indices. The continuously differentiable functions $h_i: \mathbb{R}^n \to \mathbb{R}$, satisfy the MFCQ at $x$ if there exists a vector $v \in \mathbb{R}^n$ such that $\nabla h_i(x)^\top v > 0$ for all $i \in \mathcal{J}(x)$.

Under the MFCQ, Clarke's normal cone of prox-regular sets which are finitely represented by inequalities, can be expressed in the so-called linearized form, using the normals to the constraint at the active points, as follows. Let $K = \{x \in \mathbb{R}^n | h(x) \ge 0\}$. Suppose that the functions $h_i: \mathbb{R}^n \to \mathbb{R}$, are continuously differentiable and satisfy the MFCQ. Then Clarke's normal cone to $K$ at $x$ is equal to $\mathcal{N}_K = \{w \in \mathbb{R}^n | w = -\sum_{i \in \mathcal{J}(x)} \lambda_i \nabla h_i(x), \lambda_i \ge 0\} = -(\mathcal{T}_K(x))^*$, with Clarke's tangent cone equal to $\mathcal{T}_K = \{z \in \mathbb{R}^n | z^\top \nabla h_i(x) \ge 0, \text{ for all } i \in \mathcal{J}(x)\}$. In the case of the set $S_{n,u}^{st}$ in section 3.2, $\mathcal{N}_{S_{n,u}^{st}}(q)$ is the cone generated by $-\nabla h_{n,u,i}^{st}(q)$ for active constraints $h_{n,u,i}^{st}(q) = 0$. If $K$ is closed convex this coincides with the definitions of Convex Analysis.
# D Dynamics of the RB+pendulum system

The kinetic energy is given by:

$$T(q, \dot{q}) = \frac{1}{2}m\dot{x}^2 + \frac{1}{2}m\dot{y}^2 + \frac{1}{2}I_G\dot{\theta}^2 + \frac{1}{2}I_p\dot{\alpha}^2 + \frac{1}{2}\bar{m}\dot{x}_p^2 + \frac{1}{2}\bar{m}\dot{y}_p^2, \quad (87)$$

where $I_p$ is the pendulum moment of inertia at its gravity center (supposed to be at the middle of the bar), $\bar{m}$ is its mass. From the Lagrange equations (the potential energy is supposed to be null, hence the Lagrangian function is the kinetic energy), the term $\frac{d}{dt}(\frac{\partial T}{\partial \dot{q}})$ yields the inertia matrix:

$$M(q) = \begin{pmatrix} m + \bar{m} & 0 & -\frac{\bar{m}}{2}(a \sin(\alpha + \theta) + l \cos(\theta)) & -\frac{\bar{m}a}{2} \sin(\alpha + \theta) \\ 0 & m + \bar{m} & \frac{\bar{m}}{2}(a \cos(\theta + \alpha) - l \sin(\theta)) & \frac{\bar{m}a}{2} \cos(\alpha + \theta) \\ -\frac{\bar{m}}{2}(a \sin(\alpha + \theta) + l \cos(\theta)) & \frac{\bar{m}}{2}(a \cos(\theta + \alpha) - l \sin(\theta)) & I_G + \frac{a^2+l^2}{2} + \frac{al}{2} \cos(\alpha) & \frac{\bar{m}a}{4}(a + l \cos(\alpha)) \\ -\frac{\bar{m}a}{2} \sin(\theta + \alpha) & \frac{\bar{m}a}{2} \cos(\theta + \alpha) & \frac{\bar{m}a}{4}(a + l \cos(\alpha)) & \frac{\bar{m}a^2}{4} \end{pmatrix} \quad (88)$$

The Coriolis and centripetal forces $C(q, \dot{q})\dot{q}$ can be computed from the remaining terms in $\frac{d}{dt}(\frac{\partial T}{\partial \dot{q}})$ and in $-\frac{\partial T}{\partial q}$, such that $C(q, \dot{q}) + C^T(q, \dot{q}) = \frac{d}{dt}M(q) = (\frac{\partial m_{ij}(q)}{\partial q})_{ij}$ if the Christoffel's symbols are used. One has:

$$\frac{\partial T}{\partial x} = 0, \quad \frac{\partial T}{\partial y} = 0$$

$$\frac{\partial T}{\partial \theta} = \frac{\bar{m}}{2} \left\{ -a(\dot{\alpha} + \dot{\theta}) (\dot{x} \cos(\alpha + \theta) + \dot{y} \sin(\theta + \alpha)) + l\dot{\theta} (\dot{x} \sin(\theta) - \dot{y} \cos(\theta)) \right\} \quad (89)$$

$$\frac{\partial T}{\partial \alpha} = \frac{\bar{m}}{2} \left\{ -a(\dot{\alpha} + \dot{\theta}) (\dot{x} \cos(\alpha + \theta) + \dot{y} \sin(\theta + \alpha)) + \frac{al}{2} \dot{\theta} (\dot{\alpha} + \dot{\theta}) \cos(\alpha) \right\}.$$

We therefore deduce the vector of inertial plus external generalized forces:

$$F^x(q, \dot{q}, t) = -\frac{\bar{m}}{2} \left\{ a(\dot{\alpha} + \dot{\theta}) \cos(\alpha + \theta) - l\dot{\theta} \sin(\theta) \right\} \dot{\theta} - \frac{\bar{m}a}{2}(\dot{\alpha} + \dot{\theta}) \dot{\alpha} \cos(\alpha + \theta) + F_{ext}^x$$

$$F^y(q, \dot{q}, t) = \frac{\bar{m}}{2} \left\{ -a(\dot{\alpha} + \dot{\theta}) \sin(\alpha + \theta) - l\dot{\theta} \cos(\theta) \right\} \dot{\theta} - \frac{\bar{m}a}{2}(\dot{\alpha} + \dot{\theta}) \dot{\alpha} \sin(\alpha + \theta) + F_{ext}^y$$

$$F^\theta(q, \dot{q}, t) = -\frac{\bar{m}}{2} \left\{ a(\dot{\alpha} + \dot{\theta}) \cos(\alpha + \theta) - l\dot{\theta} \sin(\theta) \right\} \dot{x} - \frac{\bar{m}}{2} \left\{ -a(\dot{\alpha} + \dot{\theta}) \sin(\alpha + \theta) - l\dot{\theta} \cos(\theta) \right\} \dot{y} \\ -\frac{\bar{m}al}{4} \dot{\alpha} \dot{\theta} \sin(\alpha) - \frac{\bar{m}al}{4} \dot{\alpha}^2 \sin(\alpha) \\ -\frac{\bar{m}}{2} \left\{ -a(\dot{\alpha} + \dot{\theta}) (\dot{x} \cos(\alpha + \theta) + \dot{y} \sin(\theta + \alpha)) + l\dot{\theta} (\dot{x} \sin(\theta) - \dot{y} \cos(\theta)) \right\}$$

$$F^\alpha(q, \dot{q}, t) = -\frac{\bar{m}a}{2}(\dot{\alpha} + \dot{\theta})\dot{x} \cos(\alpha + \theta) - \frac{\bar{m}a}{2}(\dot{\alpha} + \dot{\theta})\dot{y} \sin(\alpha + \theta) - \frac{\bar{m}al}{4}\dot{\alpha}\dot{\theta}\sin(\alpha) \quad (90)$$
References

[1] Abadie, M.: Dynamic simulation of rigid bodies: Modelling of frictional contacts. In: B. Brogliato (ed.) Impacts in Mechanical Systems. Analysis and Modelling, no. 551 in Lecture Notes in Physics, pp. 61–144. Springer Verlag (2000)

[2] Acary, V., Brogliato, B.: Numerical Methods for Nonsmooth Dynamical Systems, Lecture Notes in Applied and Computational Mechanics, vol. 35. Springer Verlag, Berlin Heidelberg (2008)

[3] Addi, K., Brogliato, B., Goeleven, D.: A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems. Applications in electronics. Mathematical Programming A **126**(1), 31–67 (2011)

[4] Adly, S., Nacry, F., Thibault, L.: Preservation of prox-regularity of sets and application to constrained optimization. SIAM Journal on Optimization **26**(1), 448–473 (2016)

[5] Akhadkar, N., Acary, V., Brogliato, B.: Multibody systems with 3D revolute joints with clearances: an industrial case study with an experimental validation. Multibody System Dynamics **42**(3), 249–282 (2018)

[6] Anitescu, M., Cremer, J.F., Potra, F.A.: Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems. Nonlinear Dynamics **24**, 405–437 (1997)

[7] Anitescu, M., Cremer, J.F., Potra, F.A.: On the existence of solutions to complementarity formulations of contact problems with friction. In: M.C. Ferris, J.S. Pang (eds.) Complementarity and Variational Problems. State of the Art, pp. 12–21. SIAM (1997)

[8] Aubin, J.P.: Applied Functional Analysis. Wiley, New York (1979)

[9] Audren, H., Kheddar, A.: 3-D robust stability polyhedron in multicontact. IEEE Transactions on Robotics **34**(2), 388–403 (2018)

[10] Balkcom, D., Trinkle, J.: Computing wrench cones for planar rigid body contact tasks. The International Journal of Robotics Research **21**(2), 1053–1066 (2002)

[11] Baraff, D.: Issues in computing contact forces for non-penetrating rigid bodies. Algorithmica **10**(2-4), 292–352 (1993)

[12] Bernstein, D.S.: Matrix Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory. Princeton University Press (2005)

[13] Bertsekas, D.P.: Convex Optimization Theory. Athena Scientific, Belmont (2009)

[14] Bicchi, A.: Hands for dexterous manipulation and robust grasping: a difficulty road toward simplicity. IEEE Transactions on Robotics and Automation **16**(6), 652–662 (2000)
[15] Blumentals, A., Brogliato, B., Bertails-Descoubes, F.: The contact problem in Lagrangian systems subject to bilateral and unilateral constraints, with or without sliding Coulomb's friction: a tutorial. Multibody System Dynamics **38**, 43–76 (2016)

[16] Bretl, T., Lall, S.: Testing static equilibirum for legged robots. IEEE Transactions on Robotics **24**(4), 794–807 (2008)

[17] Brogliato, B.: Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems. Multibody System Dynamics **29**, 289–325 (2013)

[18] Brogliato, B.: Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction. Multibody System Dynamics **32**, 175–216 (2014)

[19] Brogliato, B.: Nonsmooth Mechanics. Models, Dynamics and Control, 3rd edn. Communications and Control Engineering. Springer International Publishing, Switzerland (2016). Erratum/addendum at https://hal.inria.fr/hal-01331565v2

[20] Brogliato, B., Goeleven, D.: Singular mass matrix and redundant constraints in unilaterally constrained Lagrangian and Hamiltonian systems. Multibody System Dynamics **35**, 39–61 (2015)

[21] Brogliato, B., Thibault, L.: Existence and uniqueness of solutions for non-autonomous complementarity dynamical systems. Journal of Convex Analysis **17**(3-4), 961–990 (2010)

[22] Caron, S., Pham, Q.C., Nakamura, Y.: ZMP support areas for multicontact mobility under frictional constraints. IEEE Transactions on Robotics **33**(1), 67–80 (2017)

[23] Chen, X., Xiang, S.: Perturbation bounds of P-matrix linear complementarity problems. SIAM Journal on Optimization **18**(4), 1250–1265 (2007)

[24] Choudhury, D., Horn, R., Pierce, S.: Quasi-positive definite operators and matrices. Linear Algebra and its Applications **99**, 161–176 (1988)

[25] Cottle, R., Pang, J., Stone, R.: The Linear Complementarity Problem. Computer Science and Scientific Computing. Academic Press (1992)

[26] Facchinei, F., Pang, J.S.: Finite-Dimensional Variational Inequalities and Complementarity Problems, vol. I. Operations Research. Springer Verlag, New York (2003)

[27] Fraczek, J., Wojtyra, M.: On the unique solvability of a direct dynamics problem for mechanisms with redundant constraints and Coulomb friction. Mechanism and Machine Theory **46**, 312–334 (2011)

[28] Gholami, F., Nasri, M., Kovacs, J., Teichmann, M.: A linear complementarity formulation for contact problems with regularized friction. Mechanism and Machine Theory **105**, 568–582 (2016)

[29] Glocker, C.: Set-Valued Force Laws: Dynamics of Non-Smooth Systems. Springer Verlag (2001)
[30] Glocker, C., Pfeiffer, F.: Dynamical systems with unilateral contacts. Nonlinear Dynamics (1992)

[31] Glocker, C., Pfeiffer, F.: Complementarity problems in multibody systems with planar friction. Archive of Applied Mechanics (1993)

[32] Goeleven, D.: Complementarity and Variational Inequalities in Electronics. Mathematical Analysis and its Applications. Academic Press (2017)

[33] Han, L., Trinkle, J., Li, Z.X.: Grasp analysis as linear matrix inequality problems. IEEE Transactions on Robotics and Automation **16**(6), 663–674 (2000)

[34] Higashimori, M., Kimura, M., Ishii, I., Kaneko, M.: Dynamic capturing strategy for a 2-D stick-shaped object based on friction independent collision. IEEE Transactions on Robotics (2007)

[35] Hiriart-Urruty, J.B., Lemaréchal, C.: Fundamentals of Convex Analysis. Grundlehren Text Editions. Springer-Verlag, Berlin Heidelberg (2001)

[36] de Jalón, J.G., Gutteriez-Lopez, M.D.: Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces. Multibody System Dynamics **30**(3), 311–341 (2013)

[37] de Jalón, J.G., Unda, J., Avello, A.: Natural coordinates for the computer analysis of multibody systems. Computer Methods in Applied Mechanics and Engineering **56**(3), 309–327 (1986)

[38] Klepp, H.J.: Existence and Uniqueness of Solutions for Accelerations for Multibody Systems with Friction. ZAMM, Z. angew. Math. Mech. **75**(1), 679–689 (1995)

[39] Laulusa, A., Bauchau, O.A.: Review of classical approaches for constraint enforcement in multi-body systems. ASME J. of Comput. and Nonlinear Dynamics **3**(1), 011,004 (2008)

[40] Leine, R.I., Glocker, C.: A set-valued force law for spatial Coulomb-Contensou friction. European Journal of Mechanics A/Solids **22**(2), 193–216 (2003)

[41] Leine, R.I., van de Wouw, N.: Stability and Convergence of Mechanical Systems with Unilateral Constraints, *Lecture Notes in Applied and Computational Mechanics*, vol. 36. Springer-Verlag, Berlin Heidelberg (2008)

[42] Lötstedt, P.: Coulomb friction in two-dimensional rigid body systems. Zeitschrift für Angewandte Mathematik and Mechanik **61**, 605–615 (1981)

[43] Lötstedt, P.: Mechanical systems of rigid bodies subject to unilateral constraints. SIAM Journal on Applied Mathematics **42**(2), 281–296 (1982)

[44] Ma, S., Wang, T.: Planar multiple-contact problems subject to unilateral and bilateral kinetic constraints with static Coulomb friction. Nonlinear Dynamics **94**, 99–121 (2018)
[45] Moreau, J.J.: Application of convex analysis to some problems of dry friction. In: H. Zorski (ed.) Trends in Applications of Pure Mathematics to Mechanics, vol. 2, pp. 263–280. Pitman Publishing Ltd. (1979)

[46] Negrut, D., Serban, R., Tasora, A.: Posing multibody dynamics with friction and contact as a differential complementarity problem. ASME Journal of Computational and Nonlinear Dynamics **13**, 014,503 (2018)

[47] Nikolic, M., Borovac, B., Rakovic, M.: Dynamic balance preservation and prevention of sliding for humanoid robots in the presence of multiple spatial contacts. Multibody System Dynamics **42**, 197–218 (2018)

[48] Or, Y., Rimon, E.: Analytic characterization of a class of three-contact frictional equilibrium postures in three-dimensional gravitational environments. The International Journal of Robotics Research **29**(1), 3–22 (2010)

[49] Or, Y., Rimon, E.: Characterization of frictional multi-legged equilibrium postures on uneven terrains. The International Journal of Robotics Research **36**(1), 105–128 (2017)

[50] Panagiotopoulos, P.D.: Inequality Problems in Mechanics and Applications. Convex and Non-convex Energy Functions. Birkhäuser Verlag, Basel, CH (1985)

[51] Panagiotopoulos, P.D., Al-Fahed, A.M.: Robot hand grasping and related problems: optimal control and identification. The International Journal of Robotics Research **13**(2), 127–136 (1994)

[52] Pang, J., Trinkle, J.: Stability characterization of rigid body contact problems with Coulomb friction. ZAMM Z. Angew. Math. Mech. **80**(10), 643–663 (2000)

[53] Pang, J.S., Trinkle, J.C.: Complementarity formulation and existence of solutions of dynamic rigid-body contact problems with Coulomb friction. Mathematical Programming **73**(2), 199–226 (1996)

[54] Pang, J.S., Trinkle, J.C., Lo, G.: A complementarity approach to a quasistatic multi-rigid-body contact problem. Journal of Computational Optimization and Applications **5**(2), 139–154 (1996)

[55] Pekal, M., Fraczek, J.: Comparison of selected formulations for multibody systems dynamics with redundant constraints. Archive of Mechanical Engineering **LXII**(1), 93–112 (2016)

[56] Pfeiffer, F.: Non-smooth engineering dynamics. Meccanica **51**(12), 3167–3184 (2016)

[57] Pfeiffer, F., Glocker, C.: Multibody Dynamics with Unilateral Contacts. Nonlinear Science. Wiley (1996)

[58] Rockafellar, R.: Convex Analysis. Princeton University Press (1970)

[59] Rohde, C.A.: Generalized inverses of partitioned matrices. J. Soc. Indust. Appl. Math. **13**(4), 1033–1035 (1965)
[60] Seifried, R.: Dynamics of Underactuated Multibody Systems. Solid Mechancs and Its Applica-
tions. Springer International Publishing Switzerland (2014)

[61] Seon, J.A., Dahmouche, R., Gauthier, M.: Enhance in-hand dexterous micromanipulation by exploiting adhesion forces. IEEE Transactions on Robotics 34(1), 113–125 (2018)

[62] Shi, J., Woodruff, J., Umbanhowar, P., Lynch, K.: Dynamic in-hand sliding manipulation. IEEE Transactions on Robotics 33(4), 778–795 (2017)

[63] Simeon, B.: Computational Flexible Multibody Dynamics. A Differential-Algebraic Approach.
Springer Verlag (2013). Differential-Algebraic Equations Forum

[64] Studer, C., Glocker, C.: Representation of normal cone inclusion problems in dynamics via non-linear equations. Arch. Appl. Mech. 76(5), 327–348 (2006)

[65] Tasora, A., Anitescu, M.: A matrix-free cone complementarity approach for solving large-scale, nonsmooth, rigid body dynamics. Computer Methods in Applied Mechanics and Engineering 200(5-8), 439–453 (2011)

[66] Tasora, A., Anitescu, M.: A complementarity-based rolling friction model for rigid contacts. Meccanica 48(7), 1643–1659 (2013)

[67] Tian, Y., Takane, Y.: Schur complements and Banachiewicz-Schur forms. Electronic Journal of Linear Algebra 13, 405–418 (2005)

[68] Trinkle, J.C.: On the stability and instantaneous velocity of grasped frictionless objects. IEEE Transactions on Robotics and Automation 8(5), 560–572 (1992)

[69] Trinkle, J.C., Pang, J.S., Sudarsky, S., Lo, G.: On dynamic multi-rigid-body contact problems with Coulomb friction. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 77(4), 267–279 (1997)

[70] Trinkle, J.C., Tzitzouris, J.A., Pang, J.S.: Dynamic multi-rigid-body systems with concurrent distributed contacts: Theory and examples. Philosophical Transactions: Mathematical, Physical, and Engineering Sciences 359(1789), 2575–2593 (2001)

[71] Vanderbei, R., Carpenter, T.: Symmetric indefinite systems for interior point methods. Mathematical Programming 58(1-3), 1–32 (1993)

[72] Varkonyi, P.L., Or, Y.: Lyapunov stability of a rigid body with two frictional contacts. Nonlinear Dynamics 88, 363–393 (2017)

[73] Wojtyra, M.: Joint reactions in rigid body mechanisms with dependent constraints. Mechancsm and Machine Theory 44, 2265–2278 (2009)

[74] Wojtyra, M.: Modeling of static friction in closed-loop kineatic chains–uniqueness and parametric sensitivity. Multibody System Dynamics 39, 337–361 (2017)
[75] Wojtyra, M.: The Moore-Penrose inverse approach to modeling of multibody systems with redundant constraints. In: T. Uhl (ed.) Advances in Mechanisms and Machine Science, *Mechanisms and Machine Science*, vol. 73, pp. 3087–3096. Springer Nature Switzerland AG (2019)

[76] Wojtyra, M., Fraczek, J.: Comparison of selected methods of handling redundant constraints in multibody systems simulations. ASME Journal of Computational and Nonlinear Dynamics **8**, 021,007 (2013)