Bernard Brogliato
email: bernard.brogliato@inria.fr

Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction

Keywords: Bilateral holonomic constraints, Unilateral constraints, Complementarity conditions, Coulomb friction, Tangential restitution, Painlevé paradox, Quasi-velocities, Quasi-Lagrange dynamics, Kinematic impact law, Moreau's impact law, Kinetic angles, Kinetic metric

Kinetic quasi-velocities in unilaterally constrained Lagrangian mechanics with impacts and friction

Introduction

Finite-dimensional Lagrangian mechanical systems subject to bilateral, or unilateral constraints with or without Coulomb's friction, have received a considerable attention. It is known that the choice of the generalized coordinates may be a crucial step for either feedback control design, or numerical simulation. Many different ways of transforming the Lagrange dynamics into more suitable "canonical" forms have been proposed. Among these some are based on quantities known as quasi-velocities (or non-holonomic velocities, or generalized velocities, or generalized speeds, or pseudo-velocities, or kinematic characteristics) [START_REF] Gantmacher | Lectures in Analytical Mechanics[END_REF][START_REF] Herman | A survey of equations of motion in terms of inertial quasi-velocities for serial manipulators[END_REF][START_REF] Kane | On the derivation of equations of motion[END_REF][START_REF] Kane | The use of Kane's dynamical equations in robotics[END_REF][START_REF] Kovecses | Finite and impulsive motion of constrained mechanical systems via Jourdain's principle: discrete and hybrid parameter models[END_REF][START_REF] Loduha | On first-order decoupling of equations of motion for constrained dynamical systems[END_REF][START_REF] Lurie | Analytical Mechanics[END_REF][START_REF] Neimark | Dynamics of Non-holonomic Systems[END_REF]. Others use generalized coordinate transformations [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]. The coordinate partitioning methods have almost always been introduced for systems with a set of bilateral (holonomic or non-holonomic) constraints. The case of unilateral constraints and impacts deserves, as is shown in this paper, a specific analysis because of the varying set of active constraints.

The issues related to impacts in multibody systems constitute an important field of research because of the many applications and challenges; see [START_REF] Khulief | Modeling of impact in multibody systems: an overview[END_REF][START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF] and [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF]Chaps. 4,[START_REF] Appel | Traité de Mécanique Rationnelle, Tome 2: Dynamique des Systèmes-Mécanique Analytique[END_REF][START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] for reviews of single and multiple collisions modeling approaches. Among these, so-called multiple impacts are particularly tough to model. Frémond [START_REF] Frémond | Rigid bodies collisions[END_REF] briefly analyzed the rocking block and pointed out the necessity of "distance effects" in multiple impacts, dealing later on with chains of balls [START_REF] Caselli | Collision of three balls on a plane[END_REF][START_REF] Cholet | Chocs de solides rigides[END_REF]. The kinetic quasi-velocity change of variables described in this paper was introduced in [START_REF] Brogliato | Nonsmooth Impact Mechanics: Models, Dynamics and Control[END_REF], inspired by ideas in [START_REF] Ivanov | Impacts in a system with certain unilateral couplings[END_REF][START_REF] Ivanov | On multiple impact[END_REF][START_REF] Kozlov | Billiards: A Genetic Introduction to the Dynamics of Systems with Impacts[END_REF]. A generalized impact kinematic impact law based on this particular state space transformation was also pointed out in [START_REF] Brogliato | Nonsmooth Impact Mechanics: Models, Dynamics and Control[END_REF]Eq. (6.86)]; see also [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF]Eq. (6.96)]. Glocker [START_REF] Glocker | An introduction to impacts[END_REF] later interpreted Frémond's idea as a restitution matrix, recovering (as far as chains of balls are concerned) the normal part of the generalized law presented in the sequel of this paper: this restitution matrix should not in general be diagonal, but should contain off-diagonal entries representing the "distance effects". More recently, it has been shown that even more may sometimes be needed [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF], because the normal part of the restitution law may not be sufficient to describe some motions like perfect rocking of a planar block. Recent results on multiple impacts with or without friction have also been obtained in [START_REF] Bowling | Energetically consistent simulation of simultaneous impacts and contacts in multibody systems with friction[END_REF][START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF][START_REF] Modarres Najafabadi | Generalization of the energetic coefficient of restitution for contacts in multibody systems[END_REF][START_REF] Modarres Najafabadi | Impacts in multibody systems: modeling and experiments[END_REF][START_REF] Payr | Oblique frictional impact of a bar: analysis and comparison of different impact laws[END_REF][START_REF] Rodriguez | Solution to indeterminate multipoint impact with frictional contact using constraints[END_REF][START_REF] Seghete | Variational solutions to simultaneous collisions between multiple rigid bodies[END_REF][START_REF] Seghete | Conditions for uniqueness in simultaneous impact with application to mechanical design[END_REF][START_REF] Zhang | Study of the planar rocking-block dynamics with Coulomb friction: critical kinetic angles[END_REF].

In this paper, it is shown that the kinetic quasi-velocity change of variables is very well suited for unilaterally constrained Lagrangian systems, and allows one to analyze these systems as a generalized multi-constrained particle with which two Riemannian metrics are associated: one in a normal space and one in a tangent space. In the frictionless case, the analogy with a particle is clear. In the case with friction, the analogy is less obvious because of normal/tangential coupling terms. It is nevertheless still possible to derive a general criterion for energy consistency at impacts, involving "small couplings" and "small friction". This criterion (see Theorem 2 below) is not just an existence result but is a constructive result since all the bounds can be computed for a given system.

The paper is organized as follows: Sect. 2 presents the quasi-Lagrangian dynamics in the frictionless (Sect. 2.1) and frictional (Sect. 2.2) cases, and ends with several examples in Sect. 2.3. Section 3 is dedicated to the analysis of kinematic impact laws, in the frictionless (Sect. 3.1) and frictional (Sect. 3.2) cases. Issues of energetical consistency, uniqueness of coefficients of restitution, link with Moreau's impact law, tangential restitution and Coulomb's friction at impacts are analyzed. In passing, we provide a correct version of a result stated in [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF], where some crucial assumptions are missing. Also a general and constructive criterion for the energetical consistency when normal restitution and Coulomb's friction are used, is proposed. In Sect. 4, the contact LCP well-posedness is studied. Conclusions end the paper in Sect. 5, and some useful mathematical results are recalled in Appendix A.

Notation and definitions [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF][START_REF] Cottle | In: The Linear Complementarity Problem[END_REF][START_REF] Lancaster | The Theory of Matrices[END_REF][START_REF] Rockafellar | Convex Analysis[END_REF] w ∈ R n + if all components of w are non-negative. I generically denotes the identity matrix of appropriate dimension. A real matrix A ∈ R n×n is positive definite if x T Ax > 0 for all x = 0, it is non-negative if all its entries are nonnegative, its Moore-Penrose pseudo-inverse is denoted as A † , it is a P-matrix if all its principal minors are positive, it is a Z-matrix if its off-diagonal entries are all non-positive, it is a K-matrix if it is a Z-matrix and also a P-matrix. It is a co-positive matrix if x T Ax 0 for all x ∈ R n + , and a strictly co-positive matrix if x T Ax > 0 for all x ∈ R n + , x = 0. Let A ∈ R n×m be a real matrix, then its Euclidean norm is

A 2 = ( n i=1 m j =1 A 2 ij ) 1 2 .

Let

A ∈ R n×n be a real matrix, then

Q A = {z ∈ R n |0 z ⊥ Az 0}. Its dual is Q * A = {z ∈ R n |z T x
0 for all x ∈ Q A }. Eigenvalues are denoted as λ i (A). The smallest eigenvalue of A is denoted λ min (A), the largest one as λ max (A), whereas σ max (A) is its largest singular value (σ 2 i (A) = λ i (AA T ) = λ i (A T A)). The multi-valued sign function R → [-1, 1] is defined as sgn(x) = 1 if x > 0, sgn(x) = -1 if x < 0, sgn(0) = [-1, 1]. Let R n×n M = M T > 0, then x M = √ x T Mx, and x, y M = x T My, for all x, y ∈ R n . [a ij ] denotes the matrix with entries a ij . R(A) = {z|∃y such that z = Ay} is the range of the matrix A. For a given convex, non-empty set K ⊆ R n , and a matrix M = M T > 0, its normal cone at x in the metric defined by M is denoted N M K (x) = {w ∈ R n | w, zx M 0 for all z ∈ K}. If the metric is the Euclidean one (i.e. M = I ), then we denote the normal cone as N K (•), and

N M K (•) = M -1 N K (•).
The indicator function of a set K is defined as ψ K (x) = 0 if x ∈ K, ψ K (x) = +∞ if x / ∈ K. Let f (•) be a convex, proper, lower semi-continuous function. Then ∂f (•) is its sub-differential of convex analysis. In particular, for a convex non-empty K, one has ∂ψ K (•) = N K (•).

The kinetic quasi-velocities

Let us consider a Lagrangian system with n degrees of freedom, configuration space C, n-dimensional generalized coordinate vector q with independent coordinates q i , i ∈ {1, n}, inertia matrix M(q) = M T (q) > 0, and potential energy U(q). The tangent space to C at q is denoted as T q C q, its dual space is denoted as T * q C F , where F are the forces that work on q dt. In this paper, two kinds of time-independent constraints are considered: holonomic bilateral constraints, and unilateral constraints. The first take the form f (q) = 0, f : C → R m , while the second take the form h(q) 0, h : C → R p . It is understood that the coordinates q i are independent when the system is unconstrained, i.e., the system may be treated as a generalized particle, evolving freely on the configuration manifold C. Obviously, the addition of bilateral holonomic constraints may add some dependency between some coordinates. Generally speaking, quasi-velocities v are defined as

v = A(q) q (1)
where v has dimension n, and A(q) is invertible, but not necessarily integrable. In other words, there does not necessarily exist any quasi-position q = g(q) such that d q dt = ∂g ∂q (q) q, so that A(q) is not the Jacobian of any mapping g(q). It is clear that v may correspond to some non-holonomic constraints, hence the name non-holonomic velocities that is sometimes given to quasi-velocities. A very generic form of the unconstrained dynamics follows [START_REF] Blajer | A geometric unification of constrained system dynamics[END_REF][START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF]:

q = A -1 (q)v M(q) v = H (q, v, t) (2)
where H (q, v, t) gathers inertial forces (centrifugal, Coriolis), forces that derive from the potential energy (gravity, elasticity), and external forces (control inputs, disturbances). Notice that the mass matrix M(q) in ( 2) is not necessarily equal to M(q). This is the case for instance if (2) represents the Newton-Euler dynamics of an unconstrained 3-dimensional rigid body with rotation parameterization using Euler angles. In other cases, M(q) may be the identity, showing in passing that in general the dynamics in (2) is not in the Lagrange formalism. Quasi-velocities may also be defined as the projection of the instantaneous angular and gravity center velocities on the axes of some suitable frame [START_REF] Kane | The use of Kane's dynamical equations in robotics[END_REF]. Usually quasivelocities are applied to decouple and simplify, in a way, the dynamics, with diagonalized Lagrangian dynamics and simplified kinetic energy form [START_REF] Buttazzo | On the approximation of the elastic bounce problem on Riemannian manifolds[END_REF][START_REF] Herman | A survey of equations of motion in terms of inertial quasi-velocities for serial manipulators[END_REF]. This is also the case, as shown below, of the quasi-velocity that is analyzed in this paper. Clearly, one does not have in general v ∈ T q C for some generalized coordinate q, except if A(q) is a Jacobian in which case v = q.

The basic idea for defining the kinetic quasi-velocities as introduced in [START_REF] Brogliato | Nonsmooth Impact Mechanics: Models, Dynamics and Control[END_REF][START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF] for the case m = 0 and p = 1, is first to embed the configuration manifold C with the kinetic metric that is the metric defined with the inertia matrix M(q): given two vectors x and y in R n , their inner product is x, y q = x T M(q)y. Thus equipped, the configuration space is a Riemannian manifold. This idea is quite natural since the kinetic energy is, together with C and U(q), a fundamental ingredient of the Lagrange dynamics under study [START_REF] Lewis | Simple mechanical control systems with constraints[END_REF], so that M(q) appears as the metric of the system. It has been used in several analysis [START_REF] Blajer | A geometric unification of constrained system dynamics[END_REF][START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF][START_REF] Brauchli | Mass-orthogonal formulation of equations of motion for multibody systems[END_REF][START_REF] Buttazzo | On the approximation of the elastic bounce problem on Riemannian manifolds[END_REF][START_REF] Dzonou | A convergence result for a vibro-impact problem with a general inertia operator[END_REF][START_REF] Lee | Passive decomposition of mechanical systems with coordinate requirements[END_REF][START_REF] Liu | A unified geometric approach to modeling and control of constrained mechanical systems[END_REF][START_REF] Modarres Najafabadi | Generalization of the energetic coefficient of restitution for contacts in multibody systems[END_REF]. This is not to be confused with so-called inertial quasi-velocities defined from some factorization of M(q) = m(q)m T (q) and A(q) = m T (q) in (1), or variants of this [START_REF] Herman | A survey of equations of motion in terms of inertial quasi-velocities for serial manipulators[END_REF][START_REF] Jain | Diagonalized Lagrangian robot dynamics[END_REF].

The kinetic quasi-velocities are designed as follows. To simplify the notation and without loss of generality, we assume that the m bilateral constraints are denoted as f i (q) = h i (q) = 0 for i ∈ {1, . . . , m}, and the p unilateral constraints are h i (q) 0 for i ∈ {m + 1, . . . , m + p}, with m + p n. We also assume that all the constraints h i (q) are functionally independent at any q ∈ C, that is the (m + p) × n gradient matrix ∇h(q) has full column rank m + p. This in particular precludes that the gradients vanish in the domain of interest on C. The m + p normal unitary vectors to the co-dimension 1 constraints manifolds Σ i = {q ∈ C|h i (q) = 0} equipped with the kinetic metric are defined as

n q,i = M -1 (q)∇h i (q) ∇h T i (q)M -1 (q)∇h i (q) (3)
Clearly, the normal vectors n q,i are independent. If m + p < n we have to complete the set (n q,m+1 , . . . , n q,m+p ) by nmp mutually independent vectors t q,i in order to make a basis. The t q,i vectors are chosen such that t q,i , n q,j M = t T q,i M(q)n q,j = 0 for all i ∈ {1, . . . , nm -p}, j ∈ {1, . . . , m + p}. We notice that t T q,i M(q)n q,j = t T q,i ∇h i (q) √ ∇h T i (q)M -1 (q)∇h i (q) so that the vectors t T q,i are orthogonal to the kinetic gradients n q,j in the kinetic metric, and orthogonal to the Euclidean gradients ∇h i (q) in the Euclidean metric. One may choose unitary vectors t q,i , i.e. t T q,i M(q)t q,i = 1. Therefore, the vectors t q,i span T q C whereas the vectors n q,i span the normal cone N Φ (q) to the admissible domain Φ of C. This admissible domain for q is defined as follows:

Φ = Φ b × Φ u with Φ b = {q ∈ C|h i (q) = 0, i ∈ {1, . . . , m}} and Φ u = {q ∈ C|h i (q) 0, i ∈ {m + 1, . . . , m + p}}. Thus, Φ b is the bilateral holonomic constraints manifold with co-dimension m, Φ b = m i=1 Σ i , whereas Φ u is the admissible domain defined by the unilateral constraints, Φ u = m+p i=m+1 Φ u,i , with Φ u,i = {q ∈ C|h i (q) 0, i ∈ {m + 1, . . . , m + p}}.
For obvious reasons, we assume that Φ u contains a ball of radius > 0. One has N Φ (q) = N Φ b (q) × N Φu (q), where

N Φ b (q) = {w ∈ R n |w = m i=1 α i n q,i , α i ∈ R}.
It is in fact possible to use more geometrical arguments to characterize the bilateral constraints following, e.g. the developments for passive decomposition in [53, Eq. (10) [START_REF] Bowling | Energetically consistent simulation of simultaneous impacts and contacts in multibody systems with friction[END_REF]], see also [START_REF] Liu | A unified geometric approach to modeling and control of constrained mechanical systems[END_REF] and [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF]. Concerning the geometrical interpretation in case of unilateral constraints, see [START_REF] Glocker | An introduction to impacts[END_REF].

Let us now define the set of active unilateral constraints as I a (q) = {i ∈ {m + 1, . . . , m + p}|h i (q) = 0}. Then the normal cone to Φ u at q in the kinetic metric is given by: 1 1 This is the set denoted T ⊥ C in [START_REF] Glocker | An introduction to impacts[END_REF]Eq. (5.11)].

N M Φu (q) = w ∈ R n | w, q -q M 0 for all q ∈ Φ u = w ∈ R n |w = i∈Ia (q)
-α i n q,i , α i 0 (4)

If I a (q) = ∅ (the system evolves strictly inside Φ u ), then following [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] we adopt the convention that N M Φu (q) = {0}. The minus sign pre-multiplying α i comes from the fact that one considers the outwards normal cone, which is generated by the vectors -n q,i . In (4), the dual space T * q C is identified with R n . The expression in (4) holds when Φ u is for instance a prox-regular set (including convex sets), and constraints are independent at all q ∈ Φ; see Theorem 6.14 in [START_REF] Rockafellar | Variational Analysis[END_REF]. This in turn is closely related to the so-called Clarke or tangential regularity of Φ u , which implies the absence of re-intrant corners [79, §6.B]. One classically defines the polar cone to N M Φu (q), which is the tangent cone defined as [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]:

T Φu (q) = w ∈ R n | w, n q,i M = 1 h i (q) M -1
w T ∇h i (q) 0, for all i ∈ I a (q) [START_REF] Appel | Traité de Mécanique Rationnelle, Tome 2: Dynamique des Systèmes-Mécanique Analytique[END_REF] where this time T q C is identified with R n , and ∇h i (q) M -1 = ∇h T i (q)M -1 (q)∇h i (q). Thus, N M Φu (q) ⊂ T * q C while T Φu (q) ⊂ T q C. Actually T Φu (q) is generated by p vectors t q,i such that t q,i , n q,i q = 0. The normal and tangent cones replace, for unilaterally constrained systems, the classical constraint subspace and its orthogonal subspace. They are both convex polyhedral cones [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF].

Let us define the

n × n matrix Ξ(q) = n T q t T q
, where n q = (n q,m+1 , . . . , n q,m+p ) and t q = (t q,1 , . . . , t q,m ). The kinetic quasi-velocities are defined as

v = qnorm qtan = Ξ (q)M(q) q (6)
where the notation norm and tan come from the fact that v in ( 6) is the Euclidean projection of the generalized momentum p = M(q) q on the basis n q and t q (equivalently the projection of q on n q and t q in the kinetic metric). One could therefore call the kinetic quasi-velocities, the mass-projected momentum. From ( 6), qnorm = n T q M(q) q has dimension m + p and qtan = t T q M(q) q has dimension nmp. Notice that the (m + p) × n matrix n T q M(q) has rows

∇h T i (q)
∇h i (q) M -1 . Thus, qnorm,i = ∇h T i (q) q ∇h i (q) M -1 , and n q = M -1 (q)∇h(q) diag(

1 ∇h i (q) M -1 ).
Remark 1 The use of the kinetic metric for the study of multiple impacts was perhaps first advocated in [START_REF] Ivanov | Impacts in a system with certain unilateral couplings[END_REF]. It is also implicitly present in Moreau's works [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques. Frictionless unilateral constraints and inelastic shocks[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] where the tangent and normal cones are defined in a generic way, independently of the metric; see also [START_REF] Glocker | An introduction to impacts[END_REF]Sect. 4] for a detailed analysis. It is also used in mathematical proofs for convergence of numerical schemes [START_REF] Dzonou | A convergence result for a vibro-impact problem with a general inertia operator[END_REF]. It is clear that as far as one analyses the system at a fixed q (like for impacts), then M(q) is constant and the metric is Euclidean.

Frictionless Lagrangian systems

The frictionless Lagrange dynamics is given by

M(q) q + F (q, q, t) = ∇h(q)λ n = m i=1 ∇h i (q)λ n,i + m+p i=m+1 ∇h i (q)λ n,i (7)
where other ingredients like complementarity conditions and impact law are disregarded for the moment, and we recall that h i (q) = 0 for all i ∈ {1, . . . , m}. The m + p vector λ n groups Lagrange multipliers associated with the constraints. Let us denote F u (q) = m+p i=m+1 ∇h i (q)λ n,i the generalized contact force associated with the unilateral constraints. Notice that for a given q one has F u (q) ∈ -M(q)N M Φu (q) = N Φu (q) due to the nonnegativity of the multipliers, in particular F u (q) = 0 if q is in the interior of Φ u . Let us now perform the kinetic quasi-velocity transformation of the Lagrange dynamics [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]. First notice that qnorm qtan = Ξ (q)M(q) q + d dt Ξ (q)M(q) q (8)

Pre-multiplying both sides of (7) by Ξ(q) one obtains qnorm qtan + Ξ (q)F (q, q, t) -d dt Ξ (q)M(q) q = n T q ∇h(q)λ n t T q ∇h(q)λ n [START_REF] Blajer | A geometric unification of constrained system dynamics[END_REF] Let us define λn such that λn,i = ∇h i (q) M -1 λ n,i , i.e. λn = diag( ∇h i (q) M -1 )λ n . 2 From the definition of t q,i , it follows that t T q ∇h(q)λ n = 0, therefore, (9) becomes:

qnorm + F norm (q, qnorm , qtan , t) = n T q M(q)n q λn qtan + F tan (q, qnorm , qtan , t) = 0 [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF] with obvious definitions for F norm (q, qnorm , qtan , t) and F tan (q, qnorm , qtan , t). The terms indexed by tan are not affected by the contact force and may be thought of as some kind of tangential dynamics. We may choose to call the first line of (10) the quasi-normal dynamics and the second line the quasi-tangential dynamics. The dynamics in (6) [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF] is consequently a particular case of (2). It is clear that [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF] usually is not a Lagrange dynamics since M is constant (the identity) whereas non-linear inertial forces do not vanish (such dynamics are sometimes called Lagrange's equations in quasi-velocities, or Boltzmann-Hamel equations [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF], and they may be written in a Lagrangian-like form [START_REF] Duindam | Singularity-free dynamic equations of open-chain mechanisms with general holonomic and nonholonomic joints[END_REF][START_REF] From | An explicit formulation of singularity-free dynamic equations of mechanical systems in Lagrangian form. Part one. Single rigid bodies. Part two[END_REF]). The so-called Delassus' matrix defined when m = 0 (only unilateral constraints) is equal to ∇h T (q)M -1 (q)∇h(q) [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]. The matrix n T q M(q)n q may be seen as a normalized Delassus' matrix, 3

n T q M(q)n q = n b,T q M(q)n b q n b,T q M(q)n u q n u,T q M(q)n b q n u,T q M(q)n u q (11)
Thus, the first line in [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF] can be rewritten as

qb norm + F b norm (q, qnorm , qtan , t) = n b,T q M(q)n b q λb n + n b,T q M(q)n u q λu n qu norm + F u norm (q, qnorm , qtan , t) = n u,T q M(q)n b q λb n + n u,T q M(q)n u q λu n ( 12
)
2 Notice that the assumption that the constraints are functionally independent, guarantees that the norms ∇h i (q) M -1 never vanish, so diag( ∇h i (q) M -1 ) is positive definite. 3 The Delassus' operator is sometimes called the fundamental matrix [START_REF] Brauchli | Mass-orthogonal formulation of equations of motion for multibody systems[END_REF].

where λb n collects the first m multipliers λn,i and λu n collects the last p multipliers λn,i . Since qb norm = 0 at all times because the system evolves on the co-dimension m manifold {q ∈ C|h i (q) = 0, ∇h T i (q) q = 0, i ∈ {1, . . . , m}}, the first equation in ( 12) is equal to F b norm (q, qu norm , qtan , t) = n b,T q M(q)n b q λb n + n b,T q M(q)n u q λu n . If the m × m matrix n b,T q M(q)n b q is invertible one may obtain λ b n from this equation and insert it into the second equation in [START_REF] Brauchli | Mass-orthogonal formulation of equations of motion for multibody systems[END_REF] to obtain a dynamics that no longer depends on λb n . This modifies the unilateral part of the dynamics (and in particular one obtains a new Delassus' matrix in [START_REF] Brogliato | Absolute stability and the Lagrange-Dirichlet theorem with monotone multivalued mappings[END_REF]); see [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF] for detailed analysis of the couplings between unilateral and bilateral constraints.

The kinetic energy

The kinetic energy of our system is T (q, q) = 1 2 qT M(q) q. Clearly qb norm does not play any role in the kinetic energy, being zero. We will see later that the same applies to qtan when one considers the kinetic energy variation at an impact. Let us assume that Ξ(q) has full rank n. One has

Ξ (q)M(q)Ξ T (q) = n T q t T q M(q)(n q t q )
= n T q M(q)n q 0 0 t T q M(q)t q [START_REF] Braun | Simulation of constrained mechanical systems. Part I. An equation of motion[END_REF] from which one deduces the inverse matrix:

Ξ -T (q)M -1 (q)Ξ -1 (q) = (n T q M(q)n q ) -1 0 0 (t T q M(q)t q ) -1 (14)
which holds provided the normalized Delassus' matrix has full rank n. Now we have:

T (q, q) = 1 2 qT M(q) q = 1 2 qT M(q)Ξ T (q)Ξ -T (q)M -1 (q)Ξ -1 (q)Ξ (q)M(q) q = 1 2 v T (n T q M(q)n q ) -1 0 0 (t T q M(q)t q ) -1 v = 1 2 qT norm n T q M(q)n q -1 qnorm + 1 2 qT tan t T q M(q)t q -1 qtan = T (q, v) (15)
Now one may use [START_REF] Bowling | Energetically consistent simulation of simultaneous impacts and contacts in multibody systems with friction[END_REF] and the Schur complement [START_REF] Brogliato | Dissipative Systems Analysis and Control. Communications and Control Engineering[END_REF]Sect. A.5] to deduce:

T (q, q) = 1 2 qu,T norm G -1 (q) qu norm + 1 2 qT tan t T q M(q)t q -1 qtan (16)
with G(q) = n u,T q M(q)n u qn u,T q M(q)n b q n b,T q M(q)n b q -1 n b,T q M(q)n u q [START_REF] Brogliato | Absolute stability and the Lagrange-Dirichlet theorem with monotone multivalued mappings[END_REF] If m = 0 (no bilateral constraints) and p = 1, then qu norm = qnorm and one recovers the result in [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF]Eq. (6.11)] that T (q, q) = 1 2

q2 norm + 1 2 qT
tan (t T q M(q)t q ) -1 qtan .

It is noteworthy that the basis (n q , t q ) is not ortho-normal, because the vectors n q,i , i ∈ {1, . . . , m + p}, and t q,i , i ∈ {1, . . . , nm -p} are not necessarily orthogonal to one another (except if the constraints are orthogonal). Thus, despite the quasi-mass matrix M(q) in ( 10) is the identity, the kinetic energy in [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF] does not have the simple form 2T (q, v) = v T v as is for instance the case in [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF]Eq. (20)].

Remark 2 Let h u (q) = (h m+1 (q), . . . , h m+p (q)) T . The multiplier vector λ u n satisfies the complementarity conditions λ u n 0, h u (q) 0, (λ u n ) T h u (q) = 0, written compactly as 0 λ u n ⊥ h u (q) 0. If qnorm is integrable, i.e. there exists q norm such that d dt q norm = qnorm , then the complementarity conditions can be written as 0 q norm ⊥ λ u n 0. It is, however, always true that on a contact time interval where h i (q) = 0 persistently, then qnorm,i (t + ) 0 and qnorm,i (t + ) 0, where we implicitly suppose that velocities and accelerations are rightcontinuous functions of time. This will allow us to construct the so-called contact linear complementarity problem (contact LCP) to calculate λ u n .

The normal contact forces power

Let us denote F c (q) = ∇h(q)λ n . Let us investigate now the power performed by the generalized contact force: P n = F T c (q) q where q is assumed to be compatible with the bilateral and the unilateral constraints (i.e. we consider virtual velocities q such that the virtual displacement δq = qdt is compatible with the constraints and such that the virtual work is W n = P n dt). Then from the above developments we obtain:

P n = F T c (q) q = λ T n ∇h T (q) q = λ T n diag ∇h i (q) M -1 n T q M(q) q = λT n n T q M(q) q = λT n qnorm = λu,T n qu norm ( 18
)
where we used that qb norm = 0 always. Now one has 0 λ u n ⊥ h u (q) 0, therefore, if the system lies in the interior of the admissible domain Φ u one has λ u n = 0 and P n = 0. If the system evolves smoothly on a part of the boundary bd(Φ u ) that is finitely represented by the active constraints indexed in I a , one has 0 qnorm,i ⊥ λ u n,i 0 for all i ∈ I a . Consequently, in this case also P n = 0. Since the constraints are all perfect, the power developed by the contact forces outside possible impacts is always zero, as expected. The interest of ( 18) is to highlight the fact that the "forces" that perform work on the quasi-velocities qu norm are the multipliers λu n . This may be useful when dealing with impacts where the work done by the normal force during the collision has to be calculated (compare [START_REF] Brogliato | Dissipative Systems Analysis and Control. Communications and Control Engineering[END_REF] and ( 47)).

Let us denote F c norm (q) = n T q M(q)n q λn and D n (q) = (n T q M(q)n q ) -1 . Then from (18), one gets

P n = λT n qnorm = F c norm (q), qnorm Dn ( 19
)
Let us also denote D t (q) = (t T q M(q)t q ) -1 . 4 As a result, one finds that the frictionless Lagrangian system [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF] with a set of holonomic bilateral and unilateral constraints is equivalently represented as a generalized particle with dynamics qnorm + F norm (q, qnorm , qtan , t) = F c norm (q) qtan + F tan (q, qnorm , qtan , t) = 0 [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF] and the kinetic metric D(q) = diag(D n (q), D t (q)) (see [START_REF] Brogliato | On the control of finite-dimensional mechanical systems with unilateral constraints[END_REF]), while qnorm dt performs work on F c norm (q) in the metric of D n (q).

The kinetic angles

The kinetic angles play a crucial role in systems with multiple unilateral constraints [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF][START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF][START_REF] Gale | An indeterminate problem in classical mechanics[END_REF][START_REF] Paoli | Continuous dependence on data for vibro-impact problems[END_REF]. The kinetic angle between two active constraints i and j at q is given by

θ ij (q) = π -arccos ∇f i (q) T M -1 (q)∇f j (q) ∇f i (q) T M -1 (q)∇f i (q) ∇f j (q) T M -1 (q)∇f j (q) ( 21
)
Kinetic angles are quantities that reflect the couplings between the inertial properties and the geometrical properties of the system with unilateral constraints. It readily follows that

n T q M(q)n q = cos(π -θ ij ) = -cos(θ ij ) (22)
In particular, the diagonal entries arecos(θ ii ) = 1. One sees from (10) (11) that the normal part of the dynamics is strongly influenced by the kinetic angles matrix. This is known to play a crucial role in existence and uniqueness of solutions [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF], continuous dependence on initial data [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF][START_REF] Gale | An indeterminate problem in classical mechanics[END_REF][START_REF] Paoli | Continuous dependence on data for vibro-impact problems[END_REF], couplings between unilateral and bilateral constraints [START_REF] Brogliato | Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems[END_REF], and dynamics at impacts [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF].

Systems with Coulomb's friction

To simplify the analysis, from now on we assume that m = 0, that is there are no bilateral constraints, and all the unilateral constraints are rough with Coulomb's frictional effects. The Lagrange dynamics with tangential effects at the contacts have the generic form:

M(q) q + F (q, q, t) = ∇h u (q)λ u n + H T (q)λ t ( 23
)
where ∇h u (q)λ u n = p i=1 ∇h i (q)λ n,i , λ t collects the p multipliers that represent tangential contact forces at the contact points, and H T (q) is an n × p transformation matrix from the local frames at the contact points to the configuration space [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]Chap. 3] [34,[START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF][START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF]. Premultiplying both sides of (23) by Ξ(q), one obtains:

qnorm - d dt n T q M(q) q + n T q F (q, q, t) = n T q M(q)n q λu n + n T q H T (q)λ t qtan - d dt t T q M(q) q + t T q F (q, q, t) = t T q H T (q)λ t ( 24
)
It is remarkable in [START_REF] Cholet | Chocs de solides rigides[END_REF] that there is no reason in general that n T q H T (q) = 0, i.e. n q is not in general an annihilator of H T (q). This means that the quasi-tangential dynamics may influence the quasi-normal dynamics, but the reverse never holds since by construction of the basis (n q , t q ) one has t T q ∇h(q) = 0. We will see later that this explains some of the difficulties encountered when facing Painlevé paradoxes.

The contact LCP

Let us construct the contact complementarity problem for [START_REF] Cholet | Chocs de solides rigides[END_REF]. This is done using the complementarity conditions: 0 λ u n ⊥ qnorm 0 on phases of persistent contact with h u (q) = 0 and qnorm = 0:

0 λu n ⊥ n T q M(q)n q λu n + n T q H T (q)λ t + d dt n T q M(q) q -n T q F (q, q, t) 0 (25) Inserting Coulomb's model λ t,i ∈ -μ i λ n,i sgn(v t,i
), where v t,i = H T T,i (q) q is the relative tangential velocity at the contact point i, one obtains a complementarity problem with unknown λn . Here, H T,i is the ith column of H T and sgn is the multi-valued sign function.

To obtain [START_REF] Cottle | In: The Linear Complementarity Problem[END_REF], we also used the fact that since λ u n and λu n are related with a diagonal positive definite matrix, stating the orthogonality and the non-negativity conditions with one or the other vector is equivalent. We may rewrite the contact LCP with Coulomb's friction as

0 λu n ⊥ n T q M(q)n q -n T q H T (q) diag μ i ξ i ∇h i (q) M -1 λu n + d dt n T q M(q) q -n T q F (q, q, t) 0 ( 26
)
with ξ i ∈ sgn(v t,i ). The existence and uniqueness of solutions to this problem depends on the matrix D(q, μ, v t ) = (n T q M(q)n qn T q H T (q) diag(

μ i ξ i ∇h i (q) M -1
)) being or not a P-matrix. In the following, we shall denote b(q, q, t) = d dt (n T q M(q)) qn T q F (q, q, t).

Remark 3 One has qtan = t T q M(q) q and v t = H T T (q) q. As alluded to above, t q spans T q C so that qtan is a projection of the generalized momentum on T q C. The vector v t collects the tangential velocities in local Euclidean frames at the contact points. There is in general no reason that qtan and v t have the same physical units. There are even examples where only normal local velocities are defined, while the quasi-tangent dynamics exists (see the aligned chain of balls example below).

The tangential contact forces power

Let us compute the virtual power developed by the tangential forces. Let us denote

F t (q) = Ξ(q)H T (q)λ t = F t norm (q) F t
tan (q) . Then:

P t = qT H T (q)λ t = v T Ξ -T (q)M -1 (q)Ξ -1 (q)Ξ (q)H T (q)λ t = v, F t (q) D = qnorm , F t norm (q) Dn + qtan , F t tan (q) Dt (27)
Thus, the total virtual power of the contact forces of the dynamics in ( 24) is equal to

P = F c norm (q), qnorm Dn + qnorm , F t norm (q) Dn + qtan , F t tan (q) Dt ( 28
)
The positive definite matrices D n (q) and D t (q) define natural metrics for the system analyzed in kinetic quasi-velocities. The coupling between normal and tangential directions appears in the second term in [START_REF] Frémond | Rigid bodies collisions[END_REF]. There is no orthogonality of the quasi-generalized contact forces

F t norm (q) F t tan (q)
and F c norm (q) 0 in the inner product defined by the metric D(q). This is in contrast with what happens at the local kinematics level [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF] at the p contact points. Let us denote the ortho-normal local frame at contact point i as (n i , t i,1 , t i,2 ), with n i ∈ R 3 , t i,j ∈ R 3 . One has n i , t i,j = 0 in the Euclidean metric. Each contact force can be denoted as

F c i = F c i,n + F c i,t with F c i,n = f i,n n i and F c i,t = f i,t,1 t i,1 + f i,t,2 t i,2 . The p Coulomb's cones are denoted as C i , with F c i ∈ C i . Let U i ∈ R 3 be the local velocity, decomposed naturally as U i = U i,n + U i,t = u i,n n i + u i,t,1 t i,1 + u i,t,2 t i,2
. We may thus define virtual local velocities that are compatible with the constraints, and the virtual power at contact i is given by

P i = U i , F c i = U i,n , F c i,n + U i,t , F c i,t , while P = p i=1 P i,n + P i,t = P n + P t (29)
Thus, in the local kinematics there is a decoupling between tangential and normal virtual powers, which does not transport very well into generalized frameworks, because of the term n T q H T (q) in ( 24). Notice that if u i,n = ∇h T i (q) q, then the multipliers λ n in (7) satisfy λ n,i = f i,n , and thus P n in [START_REF] Brogliato | Dissipative Systems Analysis and Control. Communications and Control Engineering[END_REF] and P n in ( 29) are the same.

Examples

In this section, several simple systems are presented and recast into the kinetic quasi-velocity framework.

Chains of balls

Let us consider a chain of N aligned balls with masses m i , radii R and positions q i . There are N -1 unilateral constraints h i (q) = q i+1q i -2R 0. Thus ∇h i (q) = (0 . . . 0 -1 1 0 . . . 0) T , where the -1 is at the ith place and the 1 is at the (i + 1) place,

∇h T i M -1 ∇h i = m i +m i+1 m i m i+1 , n q,i = m i m i+1 m i +m i+1 (0 . . . 0 -1 m i 1 m i+1 0 . . . 0) T , t q = 1 N i=1 m 2 i × (m 1 m 2 . . . m N ) T . One finds qnorm,i = m i m i+1
m i +m i+1 ( qi+1 -qi ), and qtan =

1 N i=1 m 2 i N
i=1 m i qi . The normalized Delassus matrix has entries (i, i + 1) and

(i + 1, i) given by n T q,i Mn q,i+1 = -1 m i m i m 2 i+1 m i+2 (m i +m i+1 )(m i+1 +m i+2 )
, with diagonal entries equal to 1, and other entries zero.

Rocking block

Let us consider the planar system in Fig. 1(a). The block is a three degrees of freedom planar homogeneous solid, with generalized coordinates q = (x y θ) T , where x and y and the horizontal and vertical positions of the center of gravity G, θ is the angular position. When tangential effects act at the contact points, the dynamics is [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF]:

m ẍ(t) = λ t,1 (t) + λ t,2 (t) m ÿ(t) = λ n,1 (t) + λ n,2 (t) -mg I G θ(t) = λ n,1 (t) l 2 sin θ(t) + L 2 cos θ(t) + λ n,2 (t) l 2 sin θ(t) - L 2 cos θ(t) + l 2 cos(θ ) - L 2 sin(θ ) λ t,1 + l 2 cos(θ ) + L 2 sin(θ ) λ t,2 (30)

Fig. 1 Various systems

There are two unilateral constraints:

h 1 (q) = y - l 2 cos(θ ) + L 2 sin(θ ) 0 h 2 (q) = y - l 2 cos(θ ) - L 2 sin(θ ) 0 (31)
One computes

qnorm,1 = ẏ + ( l 2 sin(θ ) + L 2 cos(θ )) θ 1 m + 1 4I G (l sin(θ ) + L cos(θ )) 2 qnorm,2 = ẏ + ( l 2 sin(θ ) -L 2 cos(θ )) θ 1 m + 1 4I G (l sin(θ ) + L cos(θ )) 2 , qtan = √ m ẋ (32)
The calculation of the dynamics in ( 24) is lengthy, however, we are not interested in the explicit calculation of the non-linear terms but rather of the right-hand side. To simplify the presentation, we adopt the notation cθ for cos(θ ) and sθ for sin(θ ), and 2 . We obtain:

f (θ) = 1 m + 1 4I G (l sin(θ ) + L cos(θ ))
qnorm,1 + F norm,1 = 1 √ f (θ) 1 + l 2 sθ + L 2 cθ 2 1 + l 2 4 s 2 θ - L 2 4 c 2 θ λ n + 1 √ f (θ) 1 4 l 2 -L 2 cθsθ + lL 4 c 2 θ -s 2 θ 1 4 (l 2 -L 2 )cθ sθ + lL 4 λ t qnorm,2 + F norm,2 = 1 √ f (θ) 1 + l 2 4 s 2 θ - L 2 4 c 2 θ 1 + l 2 sθ - L 2 cθ 2 λ n + 1 √ f (θ) 1 4 l 2 + L 2 sθcθ - lL 4 1 4 l 2 -L 2 cθsθ + lL 4 s 2 θ -c 2 θ λ t qtan = √ m(λ t,1 + λ t,2 ) (33)
The various terms appearing in the right-hand side of ( 24) can be easily identified from [START_REF] Génot | New results on Painlevé paradoxes[END_REF].

Rod on a plane

We now deal with the system depicted in Fig. 1(b). It has three degrees of freedom and q = (x y θ) T , with x andy the coordinate of the gravity center. Its dynamics is given by:

m ẍ = λ t m ÿ = -mg + λ n I θ = -l cos(θ )λ n + l sin(θ )λ t (34)
The unilateral constraint is the signed distance between the tip A and the x-axis: h(q) = yl sin(θ ) 0. One obtains:

qnorm = m 1 + 3 cos 2 (θ ) ẏ -l cos(θ ) θ , qtan,1 = √ m ẋ, qtan,2 = 1 I + ml 2 cos 2 (θ ) l cos(θ ) ẏ + θ (35)
Also, H T = (1 0 l sin(θ )) T .

Planar box over an inclined plane

Let us consider the system depicted in Fig. 1(c). One has q = (x y ϕ) T , h(q) = yl cos(ϕ)r sin(ϕ) 0, so that ∇h(q) = (0 1l sin(ϕ)r cos(ϕ)) T . The mass matrix is M = diag(m, m, J ) with J = k 2 m, where k is the radius of gyration of the block. One has

n q = 1 1 m + (l sin(ϕ)+r cos[ϕ)) 2 J 0 1 m -l sin(ϕ) -r cos(ϕ) J T t q,1 = 1 √ m (1 0 0) T , t q,2 = 1 √ J +m(l sin(ϕ)+r cos(ϕ)) 2 (0 l sin(ϕ) + r cos(ϕ) 1) T , qnorm = ẏ-l sin(ϕ)-r cos(ϕ) 1 m + (l sin(ϕ)+r cos[ϕ)) 2 J , qtan,1 = √ m ẋ, qtan,2 = 1 √ J +m(l sin(ϕ)+r cos(ϕ)) 2 (m(l sin(ϕ) + r cos(ϕ)) ẏ + φ). Also, H T (q) = (x 0 -l) T .

Disk over a plane

Let us consider a planar system made of a disk constrained as in Fig. 1(d). The disk has mass m, moment of inertia I , radius R > 0, and q = (x y θ ) T . The unilateral constraint is

h(q) = y -R 0. One has n q = 1 √ m (0 1 0) T , t q,1 = 1 √ m (1 0 0) T , t q,2 = 1 √ I (0 0 1) T ,

and

H T = 0 1 R . Thus, qnorm = √ m ẏ, qtan,1 = √ m ẋ, qtan,2 = √ I θ . The dynamics is given by qnorm = 1 √ m λ n qtan,1 = 1 √ m λ t qtan,2 = R √ I λ t ( 36
)
The tangential relative velocity at the contact point is

v t = ẋ + R θ = 1 √ m qtan,1 + R √ I qtan,2 .

Disk in an angle

Let us consider the system depicted in Fig. 1(e). The disk has mass m, moment of inertia I , radius R > 0, and q = (x y θ ) T . The two unilateral constraints are

h 1 (q) = y -R 0 and h 2 (q) = x -R 0. One has n T q,1 = (0 1 √ m 0), n T q,2 = ( 1 √ m 0 0), and t T q = (0 0 1 √ I )
, and

H T = 1 0 0 1 R R . Thus, qnorm,1 = √ m ẏ, qnorm,2 = √ m ẋ, qtan = √ I θ . The dynamics (24) is qnorm,1 = 1 √ m (λ n,1 + λ t,2 ) qnorm,2 = 1 √ m (λ n,2 + λ t,1 ) qtan = 1 √ I R(λ t,1 + λ t,2 ) (37)

Constrained pendulum

Consider now the system in Fig. 1

(g).

There is one bilateral constraint represented by a prismatic joint (the horizontal position x A of point A) and one unilateral constraint (the vertical position of the tip B), thus n = 3, m = p = 1. Let x and y be the coordinates of the gravity center, supposed to be at the middle of the bar, with mass m. One has q = (x y θ) T ,

x A = x -l 2 sin(θ ), y B = y -l 2 cos(θ ). The mass matrix is M = diag(m, m, ml 2 12 ), h 1 (q) = x -l 2 sin(θ ), h 2 (q) = y -l 2 cos(θ ).
After few calculations, one finds:

n b q = 1 1 + 3 cos 2 (θ ) ⎛ ⎜ ⎝ 1 √ m 0 -6 cos(θ) √ ml ⎞ ⎟ ⎠ , n u q = 1 1 + 3 sin 2 (θ ) ⎛ ⎜ ⎝ 0 1 √ m 6 sin(θ) √ ml ⎞ ⎟ ⎠ t q = 1 √ 2m ⎛ ⎝ 1 -1 0 ⎞ ⎠ , n T q Mn q = ⎛ ⎝ 1 -3 sin(θ) cos(θ) √ 1+3 cos 2 (θ) √ 1+3 sin 2 (θ) -3 sin(θ) cos(θ) √ 1+3 cos 2 (θ) √ 1+3 sin 2 (θ) 1 ⎞ ⎠ qb norm = √ m 1 + 3 cos 2 (θ ) ẋ - l cos(θ ) 2 θ , qu norm = √ m 1 + 3 sin 2 (θ ) ẏ + l sin(θ ) 2 θ qtan = m 2 ( ẋ -ẏ) (38)
The Delassus' matrix shows couplings between the unilateral and the bilateral constraints, which vanish when the bar is either horizontal (θ = π 2 ), or vertical (θ = 0). Friction may act at both the bilateral and the unilateral contacts, in the tangential directions v t,2 = ẋB = ẋ -l 2 sin(θ ) θ and v t,1 = ẏA = ẏ + l 2 cos(θ ) θ . Thus,

H T (q) = ⎛ ⎝ 0 1 1 0 l 2 cos(θ ) -l 2 sin(θ ) ⎞ ⎠
It is easy to see that the term n T q H T (q) is non-zero, showing that quasi-tangential forces act on quasi-normal dynamics.

Discussion and comparison with other transformations

Most of the works reviewed below deal with the bilateral (holonomic or non-holonomic) constraints only. This makes a big difference since the normal cone in (4), which accounts for unilateral constraints, and hence topology changes, does not appear. The constraint manifold is therefore always constant of co-dimension m.

Comparing [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF] with the kinetic energy expressions when inertial quasi-velocities analyzed in [START_REF] Aghili | Dynamics and control of constrained mechanical systems in terms of reduced quasi-velocities[END_REF][START_REF] Herman | A survey of equations of motion in terms of inertial quasi-velocities for serial manipulators[END_REF][START_REF] Jain | Diagonalized Lagrangian robot dynamics[END_REF][START_REF] Junkins | An instantaneous eigenstructure quasivelocity formulation for nonlinear multibody dynamics[END_REF] are used, shows that they are not equivalent to the kinetic quasivelocities because they all yield T (q, q) = 1 2 v T v. It is known also that Kane's generalized speeds are only functions of q, q and time (like projections of instantaneous angular velocities on some frame axes), whereas the inertial and kinetic quasi-velocities also depend on the dynamical parameters of the system. It is also noteworthy that the kinetic quasi-velocities are fundamentally generalized quantities since they are defined from any set of generalized coordinates and velocities. Due to this it is not possible in general to assign physical units to the kinetic quasi-velocities, which have to be understood from a geometrical point of view.

Various methods to achieve coordinate partitioning, which roughly speaking allows one to split the dynamics into a part that depends on the contact force (a "normal" part) and one part that does not (a "tangential" part), have been proposed. The transformation proposed in [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF] decouples the system with holonomic bilateral constraints. This is a quite convenient framework for control purposes [START_REF] Morarescu | Trajectory tracking control of multiconstraint complementarity Lagrangian systems[END_REF]. However, it is a generalized coordinate transformation (hence it preserves the Lagrange structure of the dynamics), usually defined only locally on C, and whose velocity splitting is different from the one in (6): indeed the "tangential" velocity may jump at impacts because of couplings with the "normal" dynamics [15, Eq. ( 63)], which is not the case for the quantity qtan .

There also exists many techniques using direct projections of the Lagrange dynamics (7) on T q C based on annihilation matrices, or null space matrices [START_REF] Aghili | A unified approach for inverse and direct dynamics of constrained multibody systems based on linear projection operator: applications to control and simulation[END_REF][START_REF] Aghili | Control of redundant mechanical systems under equality and inequality constraints on both input and constraint forces[END_REF][START_REF] Betsch | The discrete null space method for the energy consistent integration of constrained mechanical systems. Part I. Holonomic constraints[END_REF][START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF][START_REF] Loduha | On first-order decoupling of equations of motion for constrained dynamical systems[END_REF][START_REF] Yen | An efficient Newton-type iteration for the numerical solution of highly oscillatory constrained multbody dynamic systems[END_REF]. It consists of finding a matrix P (q) (the annihilator) such that P (q)∇h(q) = 0 for all q ∈ C in the bilateral manifold {q|h i (q) = 0, i ∈ {1, . . . , m}}. P (q) is constructed from a basis of T q C [8], just as t q is. Other methods [START_REF] Aghili | A unified approach for inverse and direct dynamics of constrained multibody systems based on linear projection operator: applications to control and simulation[END_REF][START_REF] Aghili | Control of redundant mechanical systems under equality and inequality constraints on both input and constraint forces[END_REF] use the projector I -( ∂h(q) ∂q ) † ∂h(q) ∂q to obtain dynamics independent of the contact force. This may be seen as a particular case of more general projectors as defined in [12, Eq. ( 31) [START_REF] Gantmacher | Lectures in Analytical Mechanics[END_REF]], who proposes to use projections in the kinetic metric for the analysis of systems with bilateral constraints. It can be verified that the constrained velocity u β in [12, Eq. ( 27)] (called H c q in [START_REF] Modarres Najafabadi | Generalization of the energetic coefficient of restitution for contacts in multibody systems[END_REF][START_REF] Modarres Najafabadi | Impacts in multibody systems: modeling and experiments[END_REF]) is equal to n q diag( ∇h i (q) M -1 )(∇h T (q)M -1 (q)∇h(q)) -1 diag( ∇h i (q) M -1 ) qnorm . Hence, it coincides with our kinetic quasi-velocity qnorm if there is a single constraint (i.e. u β = qnorm n q ) or if the constraints are mutually orthogonal in the kinetic metric so that the matrix (∇h T (q)M -1 (q)∇h(q)) -1 is diagonal. Similar arguments and projectors are used in [START_REF] Braun | Simulation of constrained mechanical systems. Part I. An equation of motion[END_REF][START_REF] Modarres Najafabadi | Generalization of the energetic coefficient of restitution for contacts in multibody systems[END_REF][START_REF] Modarres Najafabadi | Impacts in multibody systems: modeling and experiments[END_REF]; see also [START_REF] Liu | Another form of equations of motion for constrained multibody systems[END_REF]. Actually these papers consider orthogonal (in the kinetic metric) projections onto the co-dimension m constrained manifold, whereas we rather consider projections on each vector of the basis (n q , t q ). In [65, Sect. 5], the same projection of the generalized momentum as in ( 6) is used (without normalization of n q,i ). It is pointed out the advantage of such a projection over the orthogonal projection used in [65, Sects. 2, 3] or [START_REF] Brauchli | Mass-orthogonal formulation of equations of motion for multibody systems[END_REF] for impact modeling. As clarified here, the key is the existence of the normal cone in (4) for unilaterally constrained Lagrangian systems.

Application to kinematic impact law design and analysis

Let us assume for simplicity that there are no bilateral constraints (i.e. m = 0). Thus, G(q) = D -1 n (q) and in the sequel we shall use both notations equally. It is nevertheless possible to consider m > 0, using the matrix G(q) in ( 17) and assuming that it is positive definite. We also assume that qnorm is constructed from the active constraints at the impact time t , i.e. with the constraints whose index belongs to I a (q(t)). We denote p = card(I a (q(t))). 5 From classical arguments in frictionless non-smooth mechanics, the impact dynamics at an instant t such that there is at least one i ∈ {1, . . . , p } such that qnorm,i (t -) < 0 and h i (q(t)) = 06 is given by (using (10)): qnorm t + -qnorm t -= n T q M(q)n q pn (t) qtan t + -qtan t -= 0, (39) 5 For simplicity of notation the dependence of p on q is not recalled.

where pn,i (t) = ∇h i (q) M -1 p n,i (t), i.e. pn (t) = diag( ∇h i (q) M -1 )p n (t), and p n,i (t) is the impulse of the contact force multiplier λ n,i at the impact instant t . More rigorously λ n,i is a measure at t and p n,i (t) is its density with respect to the Dirac measure at the atom t . The role played by the projection of the generalized momentum on the basis t q clearly appears in [START_REF] Hiriart Urruty | A variational approach to copositive matrices[END_REF]: the quasi-velocities qtan are conserved at the impacts when friction is absent (the constraints are said perfect). It is important to notice that despite there may be qnorm,i (t -) = 0 for some i ∈ I a (q(t)), all the terms qnorm,i , i ∈ {1, . . . , m} may undergo a jump because of the inertial couplings between the constraints, as reflected by the normalized Delassus' matrix n T q M(q)n q , which is not diagonal in general. This is for instance the case for chains of balls, or the planar rocking block (see [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF] for detailed analysis of these two cases).

It readily follows from the impact dynamics in ( 39) and ( 16) that the kinetic energy loss T L (t) = T (q, q(t + )) -T (q, q(t -)) at a time t of impact is given by:

T L (t) = 1 2 qu,T norm t + G -1 (q) qu norm t + - 1 2 qu,T norm t -G -1 (q) qu norm t - (40)
where q denotes q(t), and due to the absence of bilateral constraints G(q) = n u,T q M(q)n u q . From now on, we will drop the upper-script u since there are no bilateral constraints. The framework in [START_REF] Hiriart Urruty | A variational approach to copositive matrices[END_REF] is suitable to formulate a kinematic impact law as [START_REF] Brogliato | Nonsmooth Impact Mechanics: Models, Dynamics and Control[END_REF][START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF][START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF]:

v t + = qnorm (t + ) qtan (t + ) = -E qnorm (t -) qtan (t -) ( 41
)
where E is a generalized n × n restitution matrix. Its entries will be named the coefficients of restitution. Let us decompose it as

E = E nn E nt E tn E tt ( 42
)
with obvious dimensions of the four submatrices: -p ) . In the frictionless case, one has qtan (t + ) = qtan (t -) for any pre-impact velocity qnorm (t -), so necessarily E tn = 0 and E tt = -I . The restitution law in [START_REF] Ivanov | Impacts in a system with certain unilateral couplings[END_REF] is very general as the next result shows.

E nn ∈ R p ×p , E tt ∈ R (n-p )×(n
Proposition 1 Suppose that at least one component of qnorm (t -) or of qtan (t -) is non-zero. Then given any post-impact kinetic quasi-velocity, there exists a value of E such that (41) is satisfied. If at least one component of qnorm (t -) is negative, then there exists a value of E nn such that qnorm

(t + ) = -E nn qnorm (t -).
Proof Without loss of generality suppose that v 1 (t -) = 0 while v i (t -) = 0 for all i 2. Then it suffices to choose i1 = -v i (t + ) v 1 (t -) .

Remark 4

The fact that the quantity qtan is continuous at impacts (implying a constraint on the restitution matrix in (41) [START_REF] Ivanov | On multiple impact[END_REF], like E tn = 0 and E tt = -I ) was first noticed, in the case when a single constraint is given by q 1 0, in [START_REF] Appel | Traité de Mécanique Rationnelle, Tome 2: Dynamique des Systèmes-Mécanique Analytique[END_REF], and in a more general setting in [START_REF] Sun | The equations of motion of a system under the action of the impulsive constraints[END_REF].

Remark 5 Usually the unilateral constraints are constructed from the gap functions at the contact points, yielding the so-called local kinematics [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]Sect. 3.3]. In such a case, the quantities ∇h T i (q) q are the gap functions time-derivatives. Hence, setting E nn = diag(e n,i ) boils down to applying a Newton's restitution law at each impact point.

Frictionless systems

In the following, we examine the relation with Moreau's impact law, the dissipativity conditions, and the issue related to the uniqueness of the restitution coefficients.

Moreau's impact law

This impact law has been proposed in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], and is used in granular matter simulation [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF]. It has been deeply analyzed in [START_REF] Glocker | An introduction to impacts[END_REF]. This law can be written as [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF]:

q t + = q t --(1 + e) proj N M Φu (q); q t - (43)
where the projection is understood here as the Euclidean orthogonal projection in view of the definition of the normal cone in (4). The coefficient e is a global restitution coefficient. Let us write [START_REF] Ivanov | Singularities in the dynamics of systems with non-ideal constraints[END_REF] in terms of the vectors n q,i . For this, let us use Moreau's lemma of the two cones in the kinetic metric. Any vector z ∈ R n can be decomposed uniquely as z = x + y with x = proj[N M Φu (q); q(t -)] ∈ N M Φu (q) and y = proj[T Φu (q); q(t -)] ∈ T Φu (q), with x T M(q)y = 0 [40, Theorem 3.2.5]. The velocity q(t -) may be decomposed this way with y = q(t + )+e q(t -) 1+e and x = q(t -)-q(t + ) 1+e , so that q(t + ) satisfies ( 43) [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF]. Let us define the index set I a (q) ⊆ I a (q) such that proj[N M Φu (q); q(t -)] = i∈I a (q) α i n q,i , α i 0. By the two cones lemma I a (q) is uniquely defined. Using the second line of ( 5) and Moreau's viability lemma [16, Proposition 5.1], which states that q(t -) ∈ -T Φu (q), one deduces that the index set I a (q) can also be computed from the signs of the projections of q(t -) on the vectors n q,i . For instance, let e = 0. If q(t -) ∈ N M Φu (q) then I a (q) = I a (q). If q(t -) / ∈ N M Φu (q), then I a (q) ⊂ I a (q). Next, we define q norm = [ qnorm,i ] T for i ∈ I a (q). The new vector q tan is defined accordingly. Moreau's law in [START_REF] Ivanov | Singularities in the dynamics of systems with non-ideal constraints[END_REF] is equivalently written as

q norm t + = -E nn q norm t -, E nn = diag(e)
q tan t + = q tan t - (44)
and we recover that Moreau's law, under the constraint q(t -) ∈ -T Φu (q) and equal restitution coefficients, is equivalent to Newton's law applied at each contact [START_REF] Glocker | An introduction to impacts[END_REF]. Moreau's law is always kinetically and kinematically consistent, and energetically consistent whenever e ∈ [0, 1] [35, §5]. It is clear that the impact law in [START_REF] Ivanov | Impacts in a system with certain unilateral couplings[END_REF] significantly enlarges the scope of Moreau's law which spans only a half-line of the admissible post-impact velocities subspace. Chains of aligned balls are analyzed in [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF] (the 4-ball chain is treated but the analysis easily generalized to n-ball chains). In such systems, one always has qtan that is continuous at the impact because of linear momentum conservation. Depending on the contact properties (stiffness ratio, elasticity coefficient-linear or non-linear elasticity), the outcome may be quite different from what is predicted by Moreau's law. Extensions of Moreau's law are analyzed in [START_REF] Glocker | An introduction to impacts[END_REF][START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF]. In the case of no tangential effects, the law in [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF] allows for different coefficients at each contact, i.e. E nn = diag(e i ) in [START_REF] Jain | Diagonalized Lagrangian robot dynamics[END_REF] (see Eq. (5.103) in [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF]); see also [41, Eq. (1.3)] for the same idea. The above proof that (44) holds, and is equivalent to Proposition 5.6 in [START_REF] Glocker | An introduction to impacts[END_REF]. Filling-in the restitution matrix in (41) also generalizes Frémond's idea of introducing "distance effects" in multibody systems [START_REF] Cholet | Chocs de solides rigides[END_REF][START_REF] Frémond | Rigid bodies collisions[END_REF], which led Glocker to introduce what he named Frémond matrices. It is, however, noteworthy that the restitution matrix defined in [START_REF] Glocker | An introduction to impacts[END_REF]Eq. (5.36)] corresponds to the normal part E nn only (except for the normalization terms, see ( 46)), which is not sufficient to model the rocking block for which a quasi-tangential restitution E tt is necessary, as shown in [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF] (the rocking block is in fact the example treated in [START_REF] Frémond | Rigid bodies collisions[END_REF], and the necessity to introduce a tangential coefficient was also pointed out in [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF]). In case of a single contact, ( 43) is equivalently rewritten as q(t + ) = q(t -) -(1 + e) qT (t -)M(q)n q n q , that is

q t + = I -(1 + e)
M -1 (q)∇h(q)∇h T (q) ∇h T (q)M -1 (q)∇h(q) q t - (

In a more general setting, since the normal cone is closed convex, the projection in (43) belongs to the face of N M Φu (q) exposed by q(t -)proj[N M Φu (q); q(t -)] [40, p. 50]. The expression in [START_REF] Junkins | An instantaneous eigenstructure quasivelocity formulation for nonlinear multibody dynamics[END_REF] is used in [START_REF] Seghete | Variational solutions to simultaneous collisions between multiple rigid bodies[END_REF][START_REF] Seghete | Conditions for uniqueness in simultaneous impact with application to mechanical design[END_REF] with e = 1, who study a kind of binary collisions operator associating each constraint with such a linear mapping between post and pre-impact velocities.

Example 1 Let us consider a chain of three aligned balls with masses m 1 , m 2 , and m 3 , radius R > 0, coordinates q 1 , q 2 , and q 3 . The initial velocities are q1 (t -) = 1 m/s, q2 (t -) = q3 (t -) = 0 m/s. The two unilateral constraints are h 1 (q) = q 2q 1 + 2R 0, h 2 (q) = q 3q 2 + 2R 0. Initially, the three balls are in contact so that I a (q(t)) = {1, 2}. One calculates that qnorm,1 =

1 1 m 1 + 1 m 2 ( q2 -q1 ), qnorm,2 = 1 1 m 2 + 1 m 3 ( q3 -q2
). The quasi-tangent vector can be chosen as

t T q = 1 m 1 + 1 m 2 + 1 m 3 (1 1 1), so that qtan = 1 m 1 + 1 m 2 + 1 m 3 (m 1 q1 + m 2 q2 + m 3 q3 ). One has also -n T q,1 q(t -) = 1 m 1 + 1 m 2 ( 1 m
1 ) > 0 and -n T q,2 q(t -) = 0, consequently I a = I a = {1, 2}. Thus, applying [START_REF] Jain | Diagonalized Lagrangian robot dynamics[END_REF]

yields q2 (t + ) -q1 (t + ) = -e( q2 (t -) - q1 (t -)), q3 (t + ) -q2 (t + ) = -e( q3 (t -) -q2 (t -)), and m 1 q1 (t + ) + m 2 q2 (t + ) + m 3 q3 (t + ) = m 1 q1 (t -) + m 2 q2 (t -) + m 3 q3 (t -
), the linear momentum conservation. Let e = 1 and m 1 = m 2 = m 3 , then q1 (t + ) = -1 3 m/s, q2 (t + ) = q3 (t + ) = 2 3 m/s. One checks that energy is conserved. This solution is the one given by Moreau's law; see [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF]Example 5.7]. Similarly, e = 0 yields q1 (t + ) = q2 (t + ) = q3 (t + ) = 1 3 m/s, and one checks again that this is Moreau's law solution [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF]Example 5.7]. Moreau's law is unable to make balls 2 and 3 separate, as it spans only a small portion of the post-impact velocities subspace [START_REF] Glocker | An introduction to impacts[END_REF][START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]. From Proposition 1, we infer that given the above pre-impact velocity one can find E nn such that any post-impact qnorm (t + ) can be reached, by adding if necessary off-diagonal coefficients of restitution. This is closely linked with "distance effects", which are due to waves traveling inside the chains during the impact, and responsible for the dispersion of energy.

Example 2 Let us now consider a particle in the plane with m = 1, with two unilateral constraints as depicted in Fig. 1(f). At the impact time which occurs at the corner one has I a (q) = {1, 2}. Let us consider the pre-impact velocity q(t -) = (-1 1 2 ) T . One has -n T q,2 q(t -) > 0 and -n T q,1 q(t -) < 0. Therefore I a (q) = {2}. We can simply define t q as a vector orthogonal to n q,2 as shown on the figure. Then Moreau's law states qnorm,2 (t + ) = -e qnorm,2 (t -) = e while qtan (t + ) = qtan (t -). If e = 0 the particle slides on the boundary defined by h 2 (q) = 0.

Finally, let us notice that the restitution matrix introduced in [35, Eq. (5.36)] is designed from the un-normalized normal vectors ∇h i (q). Let us denote its entries ij , while the entries of E nn are e n,ij . Then ii = e n,ii (diagonal elements are the same), but

ij = ∇h i (q) M -1 ∇h j (q) M -1 e n,ij (46)
This can be shown from the equalities qnorm,i (t + ) = -p j =1 e n,ij qnorm,i (t -), and ∇h T i (q) q(t + ) = -p j =1 ij ∇h T i (q) q(t -).

Energy loss at impacts

There are various equivalent ways to express the quantity T L (t), using ( 39) and ( 40):

T L (t) = 1 2 qnorm t + + qnorm t -T pn (47)
which we may define as being the Thomson and Tait formula [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF], or

T L (t) = 1 2 qT norm t -(E nn -I ) T G -1 (q)(E nn + I ) qnorm t - (48)
or, using the symmetry of G(q) [19]:7

T L (t) = 1 2 qT norm t -E T nn G -1 (q)E nn -G -1 (q) qnorm t - (49)
or, following [START_REF] Glocker | Energetic consistency conditions for standard impacts[END_REF] and with ξ = qnorm (t

+ ) + E nn qnorm (t -): T L (t) = 1 2 pT n 2ξ -(I -E nn )G(q) pn (50)
A first general result may be stated.

Proposition 2 Suppose that -(E nn -I ) T G -1 (q)(E nn + I ) or -(E T nn G -1 (q)E nn -G -1 (q)
) are co-positive matrices. Then T L (t) 0. If they are strictly co-positive, then T L (t) < 0 for any non-zero pre-impact velocity.

Proof Due to the impact conditions, one has qnorm (t -) 0, in other words the p dimensional vector -qnorm (t -) belongs to R p + . From the definition of co-positivity, the results follow.

Notice that in view of the definition of co-positivity, the conditions of the proposition are necessary and sufficient. The major issue with Proposition 2 is that it is rather difficult in general to characterize the co-positivity of a matrix [START_REF] Hiriart Urruty | A variational approach to copositive matrices[END_REF]. Moreover, the co-positivity conditions imply some dependence of E nn on G(q), i.e. the restitution coefficients are not a priori independent on the system's configuration. Independence holds only if the constraints are orthogonal one to each other, in which case G(q) = I . Explicit criteria exist for symmetric matrices and p 4, or when the matrices have special structures.

Proposition 3 Suppose that all the unilateral constraints are orthogonal to one another, in the kinetic metric, and that E nn = diag(e). Then e 2 = T (t + ) T (t -) .

Proof Under the stated conditions, one has G(q) = I , and using (49) the result follows.

This result is obvious in the case p = 1 (a single constraint) using Eq. (6.11) in [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF]; see [START_REF] Modarres Najafabadi | Impacts in multibody systems: modeling and experiments[END_REF] for some experimental validation. Let us now state various characterizations of the kinetic energy loss.

Proposition 4 Suppose that E T

nn G -1 (q)E nn = A(q) for some matrix A(q) = A T (q) and real . Then there exists * such that E T nn G -1 (q)E nn -G -1 (q) < 0 for all ∈ [0, * ).

Proof The identity matrix perturbed by a small enough symmetric matrix is still positive definite [51, Exercise 8, p. 218]. One has for any x ∈ R p :

x T (-E T nn G -1 (q)E nn + G -1 (q))x = z T (I + G -1 2 (q)A(q)G -1 2 (q))z with z = G 1 2 (q)
x, where G 1 2 (q) is the symmetric square root of G(q). Thus, the result follows.

Proposition 5 Suppose that E nn = G(q)E T nn G -1 (q) and G(q) is positive definite. Then a necessary and sufficient condition for T L (t) 0 for any vector qnorm (t -) is that

|λ max (E n )| 1.
Proof From (49), we have T L (t) 0 for any qnorm (t -) if and only if

E T nn G -1 (q)E nn G -1 (q). Let G 1 2 (
q) be the symmetric positive definite square root of G(q). This inequality is equivalent to

G 1 2 (q)E T nn G -1 (q)E nn G 1 2 (
q) I , using Proposition 8.1.2 (xi) and (xiii) in [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]. Let us denote

B(q) = G 1 2 (q)E T nn G -1 2 (q)
. By the assumption of the proposition, we have

G -1 2 (q)E nn G 1 2 (q) = G 1 2 (q)E T nn G -1 2 (
q) so B(q) = B T (q), and since

B T (q) = G -1 2 (q)E nn G 1 2 (
q) we obtain B 2 (q) I . Using [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Lemma 8.4.1], it follows that equivalently λ max (B 2 (q)) 1, because B 2 (q) = B(q)B T (q) is positive semi definite and symmetric. Now we have that

B 2 (q) = G 1 2 (q)(E T nn ) 2 G -1 2 (q)
, and since it is a symmetric matrix one obtains

B 2 (q) = G -1 2 (q)E 2 nn G 1 2 (
q). Therefore, B 2 (q) and E 2 nn are similar matrices so they have the same eigenvalues [51, Proposition 1, p. 152]. Therefore, λ max (E 2 nn (q)) 1. Since the eigenvalues of E 2 nn are the squares of those of E nn , the result follows.

Proposition 6 Let G(q) > 0. Then T L (t) 0 if σ max (E nn ) 1 √ λmax(G(q))λmax(G -1 (q))
, which implies that σ max (E nn ) 1.

Proof The proof begins similarly to the proof of Proposition 5, and we obtain that

T L (t) 0 ⇔ B(q)B T (q) I with B(q) = G 1 2 (q)E T nn G -1 2 (q)
. By [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Lemma 8.4.1], one has equivalently λ max (B(q)B T (q)) = σ 2 max (B(q)) 1. From [7, Corollary 9.6.5], one has σ max (B(q)) σ max (G

1 2 (q))σ max (G -1 2 (q))σ max (E nn ).Therefore, σ max (G 1 2 (q))σ max (G -1 2 (q))× σ max (E nn ) 1 implies that σ 2 max (B(q)) 1.
From the symmetry and positive definiteness of G(q) and of its square root, one has σ max (G

1 2 (q)) = √ λ max (G(q)
), so the result follows. For the last statement notice that I = G -1 2 (q)G 1 2 (q), so again from [7, Corollary 9.6.5]

1 σ max (G 1 2 (q))σ max (G -1 2 (q)) = λ max (G(q))λ max (G -1 (q))
.

If E nn is symmetric, one has σ i (E nn ) = |λ i (E nn )|, 1 i p . Thus, σ max (E nn ) = sprad(E nn )
, where sprad is the spectral radius, i.e. the largest |λ i (E nn )|, 1 i p . One can therefore replace σ max (E nn ) by sprad(E nn ) in Proposition 6. If E nn has only nonnegative eigenvalues, then sprad(E nn ) = λ max (E nn ). The condition imposed in Proposition 5 holds if for instance E nn = diag(e). In fact,

E nn = G(q)E T nn G -1 (q) is equivalent to G -1 (q)E nn = E T nn G -1 (q)
, which allows us to rewrite [START_REF] Kovecses | Finite and impulsive motion of constrained mechanical systems via Jourdain's principle: discrete and hybrid parameter models[END_REF] as

T L (t) = 1 2 qT norm t -E T nn E nn -I G -1 (q) qnorm t - (51)
This shows clearly that the kinetic angles play a crucial role in the multiple impact energetical behavior. The expression in ( 51) is very similar to the classical energy loss expression for two frictionless particles that undergo a single impact (see, for instance, [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF]Eq. (4.44)]. We may therefore see [START_REF] Hiriart Urruty | A variational approach to copositive matrices[END_REF] and ( 51) as the true generalization of the simple frictionless shock between two particles, to multi-constraint frictionless Lagrange systems considered as a point moving on C. The kinetic angles matrix G(q) is a measure of the distortion due not only to the inertia matrix M(q) but also to the unilateral constraints. Another characterization of energy loss is as follows, which starts from [START_REF] Kovecses | Finite and impulsive motion of constrained mechanical systems via Jourdain's principle: discrete and hybrid parameter models[END_REF].

Proposition 7 Assume that G(q) is positive definite. If E nn 2 < 1 √ G -1 (q) 2 G(q) 2
, then the matrix G -1 (q) -E T nn G -1 (q)E nn is positive definite, and T L (t) 0 in [START_REF] Kovecses | Finite and impulsive motion of constrained mechanical systems via Jourdain's principle: discrete and hybrid parameter models[END_REF].

Proof Using [START_REF] Chen | Perturbation bounds of P-matrix linear complementarity problems[END_REF]Theorem 2.11]

(see Theorem 3 in Appendix A) with A = G -1 (q) - E T nn G -1 (q)E nn and M = G -1 (q), one obtains that M -A = E T nn G -1 (q)E nn , from which it follows that provided G(q) 2 E T nn G -1 (q)E nn 2 < 1 then A is positive definite. From [7, Proposition 9.3.5] one has E T nn G -1 (q)E nn 2 E nn 2 2 G -1 (q) 2 . So the condition G(q) 2 E nn 2 2 G -1 (q) 2 < 1 guarantees that A is positive definite.
Obviously, in Proposition 4, the value * may be linked with the entries of E nn . Roughly, it says that if the restitution coefficients are small enough the energy loss is guaranteed. When all constraints are orthogonal to one another (G(q) = I ) and choosing E = diag(e n,i ), then one finds from Proposition 7 that p i=1 e 2 n,i p . If all coefficients are equal to one another, then it follows that |e n | 1. The interest of Propositions 6 and 7 is that they clearly link the restitution matrix to the kinetic angles, i.e. to the inertial/geometrical couplings that exist in the system, with computable upper-bounds. Propositions 5 and 6 state in a correct way Proposition 1 in [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF], which wrongly asserts that |λ max (E nn )| 1 is a sufficient condition for T L (t) 0 under symmetry of E nn . The energy consistency of an extended frictionless Moreau's law with E nn = diag(e i ) is analyzed in [55, §7.1] [START_REF] Glocker | Energetic consistency conditions for standard impacts[END_REF], starting from the Thomson and Tait formula [START_REF] Kane | The use of Kane's dynamical equations in robotics[END_REF], or from [START_REF] Khulief | Modeling of impact in multibody systems: an overview[END_REF], or from [START_REF] Kozlov | Billiards: A Genetic Introduction to the Dynamics of Systems with Impacts[END_REF]. Actually, one may use Propositions 7.1 and 7.2 in [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF] to analyze [START_REF] Lancaster | The Theory of Matrices[END_REF]. The condition of Proposition 5 is quite close to the commuting conditions of [55, p. 159]. The conditions of Proposition 7 may be conservative as shown in Remark 6. Finally, let us remind that in the case Poisson coefficients are used (kinetic impact law), one obtains similar expressions for the loss of kinetic energy (see Eq. ( 43) in [START_REF] Glocker | Multiple impacts with friction rigid multi-body systems[END_REF]). The quadratic forms in ( 48)-( 50) therefore possess a general interest for both kinematic and kinetic impact laws.

Remark 6 Let us compare the criteria obtained in Propositions 2 and 7, when p = 2. We denote -ae n,1 e n,2 + a + (e 2 . Proposition 7, however, proposes a criterion that is easy to check, whereas as alluded to above, co-positivity is in general hard to check.

G(q) = 1 a a 1 hence G -1 (q) = 1 1 -a 2 1 -a -a 1 and |a| < 1. Proposition 7 yields E nn 2 < 1-a 2 2+a 2

Uniqueness of the coefficients or restitution

As shown in [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF]Sect. 7] on a 4-ball chain, a major discrepancy between single and multiple impacts is that for the latter the kinematic, kinetic, and energetic consistency constraints do not yield a unique set of coefficients of restitution for given pre-impact velocity. Proposition 8 Let a frictionless impact occur at t . It is necessary and sufficient that:

• (i) E nn is non-negative (kinematic consistency), • (ii) G -1 (q)(I + E nn ) is non-negative (kinetic consistency), • (iii) G -1 (q) -E T nn G -1 (q)E nn is co-positive (energetic consistency)
, for E nn to be an admissible restitution matrix for any pre-impact velocity qnorm (t -).

Proof (i) assures that qnorm (t + ) 0 for any qnorm (t -) 0 (kinematic consistency), (ii) guarantees that pn 0 (kinetic consistency), and (iii) is the energetical constraint.

Let us note that similar conditions have been stated in [START_REF] Payr | An Experimental and Theoretical Study of Perfect Multiple Contact Collisions in Linear Chains of Bodies[END_REF]Sect. 3.3] and [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF]. Clearly, if G -1 (q) is itself non-negative, then (i) satisfied implies that (ii) is satisfied (kinematic consistency implies kinetic consistency). This is the case for a chain of three aligned balls [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF].

Corollary 1 Suppose that the system is conservative. Then it is necessary and sufficient that:

• (i) E nn is non-negative (kinematic consistency), • (ii) E T nn G -1 (q)E nn (I + E nn ) is non-negative (kinetic consistency), • (iii) G -1 (q) = E T nn G -1 (
q)E nn (energetic consistency), for E nn to be an admissible restitution matrix for any pre-impact velocity qnorm (t -). Then the eigenvalues of E nn have modulus equal to one. Proof The proof of (i) and (ii) is a direct consequence of Proposition 8. (iii) Let B(q) = G 12 (q)E nn G -1 2 (q) . One has B T (q)B(q) = I . Being a unitary (orthogonal) matrix B(q) has its eigenvalues on the unit circle [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Proposition 5.5.25]. Since B(q) and E nn are similar, they have the same eigenvalues [51, Proposition 1].

If E nn has real eigenvalues, then they are equal to one under Corollary 1 conditions. The conditions (i) (ii) (iii) of Corollary 1 make a non-linear equation for E nn . There are p 2 unknowns (the restitution coefficients), with p (p +1) may not hold. Suppose that (i) is satisfied. From [START_REF] Caselli | Collision of three balls on a plane[END_REF], if all the kinetic angles satisfy θ ij π 2 , then (ii) is satisfied. When p = 1 (single impact) then (i) e n 0, (ii) e 2 n (1 + e n ) 0, (iii) e 2 n = 1, so that e n = 1 is the only solution (the same holds if the constraints are orthogonal so that G(q) = I , which implies that e n,i = 1 for all ∈ {1, . . . , p }). Proposition 8 and Corollary 1 are stated for all admissible pre-impact velocities. However, there are typical problems (chains of balls, rocking block) for which some components qnorm,i (t -) = 0. As alluded to in [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF], this may explicitly yield non-uniqueness of the coefficients of restitution. To illustrate this point, let us consider chains of balls (Sect. 2.3.1) with N = 3, q1 (t -) = 1 m/s, q2 (t -) = q3 (t -) = 0 m/s, and equal masses m. One has:

qnorm,1 = m 2 q2 -q1 , qnorm,2 = m 2 q3 -q2 , qtan = 1 √ 3 q1 + q2 + q3 G(q) = 1 -1 2 -1 2 1 , λn = 4 3 1 1 2 1 2 1 e n,1 + 1 e n,21 (52)
It is noteworthy that qnorm

(t + ) = e n,1
e n,21 , i.e. e n,2 and e n,12 do not play any role in this collision. The kinematic and kinetic constraints are satisfied for any non-negative coefficients. The conservation of kinetic energy yields that the two coefficients lie on the ellipse:

e 2 n,1 + e 2 n,21 -e n,1 e n,21 = 1 (53)
The non-negativity implies that the portion of this ellipse in R 2 + contains all (an infinity) the admissible restitution law parameters. It is therefore a peculiarity of multiple impacts that the kinematic, kinetic, and energetic consistencies do not allow one to always uniquely determine the restitution coefficients. This, in passing, raises the following question: What is the mechanical meaning of such coefficients? Such uniqueness issues stem from the fact that multiple impacts intrinsically are a multiscale phenomenon, where the elasticity properties at each contact strongly influence the impact outcome [START_REF] Brogliato | Nonsmooth Impact Mechanics: Models, Dynamics and Control[END_REF][START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF][START_REF] Ivanov | On multiple impact[END_REF][START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]. For a fixed energetical behavior (e.g. conservative systems) varying the stiffnesses ratios and the elasticity coefficient (linear, Hertz elasticity, or else) may induce drastic changes of the post-impact velocity [START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF]Chaps. 5,[START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF]. Kinematic impact laws do not convey enough information on the contact process to a priori detect, which one of the outcomes is the right one. Some parameter-fitting process seems unavoidable.

This analysis motivates the next definition.

Definition 1 One says that the system has the U-property if the set of conditions (i) (ii) (iii) in Corollary 1 yields a unique restitution matrix E nn (equivalently a unique set of restitution coefficients).

The orthogonality conditions n T q,i M(q)n q,j = 0 for all i = j (equivalently G(q) = I ) guarantee the U-property. It would be quite interesting to link the U-property with discontinuity of the trajectories with respect to initial conditions and the study of limits of compliant models as the stiffnesses diverge to infinity. In other words, are the following equivalences true? U-property ⇔ continuity w.r.t. initial data ⇔ uniqueness of limit for infinite stiffnesses Elements of answers may exist in the literature, using the results about a particle hitting a corner [START_REF] Ivanov | On multiple impact[END_REF][START_REF] Kozlov | Billiards: A Genetic Introduction to the Dynamics of Systems with Impacts[END_REF][START_REF] Nguyen | Multiple Impacts in Dissipative Granular Chains[END_REF][START_REF] Towne | One-dimensional collisions and Chebyshev polynomials[END_REF], where explicit formulas giving the particle evolution when it collides the angle are provided. This analysis, however, quickly become quite cumbersome and seem to be limited to simple dynamics with p = 2. A general characterization of continuity w.r.t. initial data is made in [START_REF] Paoli | Continuous dependence on data for vibro-impact problems[END_REF]. Continuity holds when the kinetic angles between the p constraints are all π 2 and the coefficients of restitution with each constraint are zero, or if the constraints are all pairwise orthogonal (in the kinetic metric). In such cases, the first equivalence holds. Indeed by continuation from the neighborhood of the singularities of the boundary of the admissible domain Φ = {q ∈ C|h(q) 0}, one may define a unique post-impact velocity for any pre-impact data. Therefore, E nn is unique.

One may think that [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF] straightforwardly shows that a frictionless Lagrangian system is equivalent, via the kinetic quasi-velocities transformation, to a particle. Even if we disregard the non-linear terms in the left-hand side in [START_REF] Blajer | A geometrical interpretation and uniform matrix formulation of multibody system dynamics[END_REF], the above shows that inertial couplings between the constraints render the impact problem quite different from the single collision case.

Systems with tangential effects

Some systems like chains of balls possess an intrinsic "invariance at impacts" of the quasitangential velocity qtan (in this case, the total linear momentum of the chain). However, in general some tangential effects exist at the contact points, in the locally tangent planes spanned by (t i,1 , t i,2 ) (see Sect. 2.2.2 for the notations).

Tangential restitution

The relation between local tangential coefficients of restitution and Coulomb's friction at the impulse level is discussed in [START_REF] Glocker | Energetic consistency conditions for standard impacts[END_REF]Sect. 4.5]. It is pointed out that despite some few experimental cases show the necessity of introducing tangential restitution coefficients, in general tangential effects are primarily due to Coulomb's friction. The relation between local tangential and quasi-tangential coefficients of restitution is established for the planar rocking block when there is perfect sticking at the two impact points in [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF]Sect. 6]. In fact, sticking imposes some constraints on the tangential velocities, and the tangential coefficients (local and quasi) impose also some constraints on tangential velocities. The mixture of all those constraints create the relations between the various coefficients.

The starting point is to consider the local kinematics at the contact points [1, Sect. 3.3], which dictates the form of H T (q) in ( 23) and [START_REF] Cholet | Chocs de solides rigides[END_REF]. When the local frames consist of a single normal vector, like in aligned chains of balls, then H T (q) = 0 and there are no tangential effects. Otherwise, H T (q) = 0. It is noteworthy that since qtan is an np vector, and since p t is a p vector (if all the contact points are rough, otherwise it has dimension less than p ), in general the tangent quasi-velocity and the local-kinematics tangent velocity do not have a clear, one-to-one, correspondence.

A first condition is the consistency of the impact law (41) (42) and of the shock dynamics qnorm t + -qnorm t -= n T q M(q)n q pn (t) + n T q H T (q)p t qtan t + -qtan t -= t T q H T (q)p t [START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF] Inserting ( 41) into ( 54), one obtains

- I + E nn E nt E tn I + E tt = Ē qnorm (t -) qtan (t -) = G(q) n T q H T (q) 0 t T q H T (q) = Ḡ(q) pn p t (55)
In the frictionless case, pn can always be calculated (independently of its sign) as pn = -G -1 (q)(I + E nn ) qnorm (t -), for any value of qnorm (t -), because G(q) has full rank. We recall that Ē is n × n, while Ḡ(q) is n × 2p since we suppose that friction acts at all the contacts. The equation in [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF] has a solution for the impulse if and only if:

Ē qnorm t -∈ R Ḡ(q) (56)
A sufficient condition is that R( Ē) ⊆ R( Ḡ(q)). Taking into account that the dimensions of p t and of qtan usually are different, there is no systematic way to translate this range condition into conditions on E and the impact geometry in Ḡ(q). Case by case analysis seem unavoidable.

Implicitly defined restitution operator Inspired by Moreau [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] and others later [START_REF] Frémond | Rigid bodies collisions[END_REF][START_REF] Frémond | Non-smooth Thermomechanics[END_REF][START_REF] Glocker | Energetic consistency conditions for standard impacts[END_REF][START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF], we may impose pn ∈ -N Vn(q) (Γ n ( qnorm

(t + ) + Λ n qnorm (t -))) and p t ∈ -N V t(q) (Γ t ( qtan (t + ) + Λ t qtan (t -))
), for some matrices Λ n , Λ t , Γ n , Γ t , and convex sets V n (q) and V t (q). Inserting this into (54), one obtains the generalized equation:

v t + -v t -∈ -Ḡ(q) N Vn(q) (Γ n ( qnorm (t + ) + Λ n qnorm (t -))) N V t(q) (Γ t ( qtan (t + ) + Λ t qtan (t -))) (57)
Defining the convex set W (q) = V n (q) × V t (q) and Λ = diag(Λ n , Λ t ), Γ = diag(Γ n , Γ t ), we get

v t + -v t -∈ -Ḡ(q) N W (q) Γ v t + + Λv t - (58)
The existence and uniqueness of a solution v(t + ) to the generalized equation in [START_REF] Liu | Another form of equations of motion for constrained multibody systems[END_REF] depend on the matrices Ḡ(q), Γ , Λ, and on the convex sets V n (q) and V t (q). Suppose for instance that there exists a symmetric positive definite matrix P such that P Ḡ(q) = Γ T , and let us denote R its symmetric square root: R 2 = P . 9 Then using the chain rule from convex analysis [START_REF] Rockafellar | Convex Analysis[END_REF], one obtains that ( 58) is equivalent to

z -Rv t -∈ -N W (q) z + RΛv t - (59)
where z = Rv(t + ) is the unknown, 10 and W (q) = {x| ḠT Rx ∈ W (q)} is a convex set (which we assume to be non-empty, which is equivalent to R( Ḡ(q)) ∩ W (q) = ∅). Then by classical arguments from convex analysis, one infers that the unique solution of the generalized equation ( 58) is given by:

v t + = -Λv t -+ R -1 proj W (q); R(Λ + I )v t - = -Λv(t -) + R -1 z * with z * = argmin z∈ W (q) 1 2 z -R(Λ + I )v t -T z -R(Λ + I )v t - (60)
where proj denotes the orthogonal projection in the Euclidean metric. The proof for ( 59) and ( 60) is given in Appendix D. Since Γ is free to the designer, fulfilling the above conditions is always possible. However it is another, much more difficult matter to design an impact law which possesses a mechanical meaning. Depending on W (q), Γ , P and v(t -), [START_REF] Lubarda | The bounds on the coefficients of restitution for the frictional impact of rigid pendulum against a fixed surface[END_REF] implicitly defines a restitution operator of the same type as [START_REF] Ivanov | Singularities in the dynamics of systems with non-ideal constraints[END_REF].

Energetic consistency (restitution matrix)

The minimum requirement one may impose on an impact law like the one in (41) [START_REF] Ivanov | On multiple impact[END_REF] is that it satisfies the three fundamental consistencies: kinematic, kinetic, and energetic.

The energetic consistency may be checked using ( 15) and [START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF]. Indeed one has three basic ingredients:

T (q, v) = 1 2 v T D(q)v, v t + = -Ev t -, v t + -v t -= Ḡ(q) p (61)
with p = pn pt and D(q) = diag(D n (q), D t (q)). The main discrepancy between the frictionless case and the case with tangential effects, is that Ḡ(q) may not be positive definite, even if the constraints are independent. Actually, G(q) is always at least positive semi definite, but Ḡ(q) may be indefinite, or negative definite. Similarly to the frictionless case in (47)- [START_REF] Kozlov | Billiards: A Genetic Introduction to the Dynamics of Systems with Impacts[END_REF], several expressions of the kinetic energy loss may be written:

T L (t) = 1 2 pT ḠT (q)(I -E)v t - (62)
or

T L (t) = 1 2 v T t -(I + E) T D(q)(E -I )v t - (63)
or

T L (t) = 1 2 v T t -E T D(q)E -D(q) v t - (64)
Deriving a similar proposition as Proposition 4, it follows that small enough restitution coefficients guarantee energetic consistency. The extension of Proposition 7 is as follows.

Proposition 9 Assume that D(q) is positive definite.

If E 2 < 1 √ D -1 (q) 2 D(q) 2
, then the matrix D(q) -E T D(q)E is positive definite, and T L (t) 0 in [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF].

One can use the block-diagonal form of D(q) to refine the result since one has

D(q) 2 = D n (q) 2 2 + D t (q) 2 2 and D -1 (q) 2 = D -1 n (q) 2 2 + D -1 t (q) 2 2 .
Modeling tangential effects with constant tangential restitution coefficients is in most of the cases a very crude approximation of the real mechanical phenomenon.

Energetic consistency (implicit restitution operator)

Starting from T L (t) = 1 2 (v(t + ) + v(t -)) T D(q)(v(t + ) -v(t -)
), a natural way to deal with its dissipativity is to force the operator defined in [START_REF] Liu | A unified geometric approach to modeling and control of constrained mechanical systems[END_REF] to be maximal monotone. Using [START_REF] Loduha | On first-order decoupling of equations of motion for constrained dynamical systems[END_REF], one obtains

T L (t) ∈ - 1 2 z t + + z t -T R -1 D(q)R N W (q) z t + + RΛR -1 z t - (65)
At this stage , we are left with the analysis of the monotonicity of y → R -1 D(q)R N W (q) (y + (RΛR -1 -I )z(t -)). Notice that another path that consists of using (60) can be followed. Let Λ = Λ T . If one assumes that 0 ∈ W (q), then the projection is contractive because W (q) is convex, so that proj[ W (q); R(Λ + I )v(t -)] R(Λ + I )v(t -) , and one obtains:

T L (t) = 1 2 v T t -ΛDΛ -D v t - + 1 2 proj T W (q); R(Λ + I )v t -R -1 DR -1 proj W (q); R(Λ + I )v t - -v T t -ΛDR -1 proj W (q); R(Λ + I )v t - - 1 2 λ min (D -ΛDΛ) v t -2 + 1 2 λ max R -1 DR -1 R(Λ + I ) 2 + 2 R -1 DΛ 2 R(Λ + I ) 2 ) v t -2 (66)
It appears that Λ plays the role of a restitution matrix that should be small enough. We do not investigate further implicitly defined kinematic restitution operators, because this has been thoroughly analyzed by Moreau, Frémond, and others, while such type of formulations is not the main topic of this paper.

Remark 7 (Necessity of E tt = -I ) The rocking block of Sect. 2.3.2 has been analyzed in [START_REF] Brogliato | Analysis of a generalized kinematic impact law for multibodymulticontact systems, with application to the planar rocking block and chains of balls[END_REF] with the above generalized kinematic law. It has been shown that Coulomb friction is unable to model some rocking motions, which necessitate the introduction of a quasitangential restitution coefficient that makes qtan jump. The obtained impact law is in turn equivalent to the angular velocity restitution that is widely used in earthquake engineering.

Coulomb's friction at contact points

A condition for the tangent quasi-velocities qtan,i to be able to correctly represent the local tangential effects is that their number, which is always np , is larger or equal to the number of independent local tangent velocities, at impact times. To illustrate this, we may consider the disk in an angle and the rocking block examples in Sect. 2.3. For the disk system, n = 3, p = 2, there are two tangential velocities v t,1 and v t,2 , which may be zero (stick) or non-zero (slip) independently. For the block system, n = 3, p = 2, one always has v t,1 = v t,2 when the two corners are in contact. One expects that restitution applied to qtan will be a crude approximation of Coulomb's friction for the disk, whereas it may be a better model for the block. Let us consider the disk example of Sect. 2.3.5 with n = 3, p = 1. As shown by several authors a correct model for the tangential velocity restitution is to apply v t (t + ) = -e t v t (t -), where the tangential restitution coefficient is related to Coulomb's friction in a certain manner (see [16, pp. 133-134] for a short review and references). Let us choose E tt = diag(e t ) and E = diag(E nn , E tt ). Then the mechanical meaning of the restitution law (41) is clear. Notice from (36) that one has n T q H T = 0, which means that the normal and the tangential quasi-dynamics are decoupled. This is not the case in neither [START_REF] Génot | New results on Painlevé paradoxes[END_REF] nor [START_REF] Glocker | Multiple impacts with friction rigid multi-body systems[END_REF].

We may assume in a rather classical way that at each contact point one has

p t,i ∈ -μ i p n,i sgn v t,i t + + e t,i v t,i t - (67)
for some coefficient e t,i , which is a way to write the Coulomb's model, at the impulse level.

Adopting the normal restitution rule qnorm (t + ) = -E nn qnorm (t -), we obtain the following inclusions:

-

(I + E nn ) qnorm t -∈ G(q) pn -n T q H T (q)[ μ][ pn ] Sgn v t t + + E tt v t t - qtan t + -qtan t -∈ -t T q H T (q)[ μ][ pn ] Sgn v t t + + E tt v t t - (68)
with:

[ μ] = diag μ i ∇h i (q) M -1 ∈ R p ×p , [ pn ] = diag( pn,i ), E tt = diag(e t,i ), v t = H T T (q)Ξ T (q)v (v is in (6)), Sgn v t t + + E tt v t t -= sgn v t,1 t + + e t,1 v t,1 t -, . . . , sgn v t,p t + + e t,p v t,p t -T
The unknowns of the generalized equation ( 68) are the p impulses pn,i , and the np quasi-velocities qtan (t + ). The generalized equation in ( 68) is quite complex, in part because it is non-associated (i.e. the multi-valued part is not the sub-differential of a convex nonsmooth potential function). It may be analyzed as follows. Suppose that all components v t,i (t + ) + e t,i v t,i (t -) > 0, i.e. the contact points are all in the same sliding status. In order to verify if such a mode is the mode of the system, one constructs the following system:

(a) -(I + E nn ) qnorm t -= G(q) -n T q H T (q)[ μ] pn (b) qtan t + -qtan t -= -t T q H T (q)[ μ] pn (c) H T T (q)M -1 (q)Ξ -1 (q) -E nn qnorm (t -) qtan (t + ) +E tt H T T (q)M -1 (q)Ξ -1 (q) qnorm (t -) qtan (t -) > 0 (d) pn 0. (69)
If the system in (69) has a unique solution for pn and qtan (t + ), then the impact is in this "allsliding" mode. Suppose that one has pn = -(G(q)n T q H T (q)[ μ]) -1 (I + E nn ) qnorm (t -) 0. If this is the case, then (69) boils down to finding qtan (t + ) such that:

qtan t + = qtan t -+ t T q H T (q)[ μ] G(q) -n T q H T (q)[ μ] -1 (I + E nn ) qnorm (t -) H T T (q)M -1 (q)Ξ -1 (q) -E nn qnorm (t -) qtan (t + ) + E tt H T T (q)M -1 (q)Ξ -1 (q) qnorm (t -) qtan (t -) > 0. (70)
The second matrix in the second line in ( 72) is skew-symmetric. It is checked that the first matrix is positive definite if and only if 0 μ < 4 L , which implies that the whole matrix is also positive definite. The inverse matrix is given by

G(q) + n T q H T (q)μ -1 = 4f (0) 1 16L 2 4 + L 2 -μlL -4 + L 2 -μlL -4 + L 2 + μlL 4 + L 2 + μlL ( 73
)
This matrix is non-negative if and only if μ 4-L 2 lL and μ L 2 -4 lL , which implies that L = 2 and μ = 0: the frictionless case. Thus, the criterion (iii) of Proposition 10 does not hold. Let us, however, focus on the rocking motion, which implies that qnorm,2 (t -) = 0 (the corner A is in contact at the collision instant at B). Suppose also that we model plastic impacts so that E nn = 0. Then qnorm (t -) = qnorm,1 (t -) 0 . Thus, only the first column of (G(q) + n T q H T (q)μ) -1 plays a role and the non-negativity condition yields 4 + L 2 -μlL 0 and L 2 + μlL 0, that is 4-L 2 lL μ 4+L 2 lL . So conditions (i), (ii), and (iii) of the proposition can be fulfilled for some values of the friction coefficient 0 μ < min( 4l , 4+L 2 lL ).

Another source of conservativeness of the conditions of Proposition 10 comes from the fact that (iii) may not be satisfied, however, there exists some pre-impact velocities such that pn 0. The next example illustrates this point.

Example 4 Consider now the system in Sect. 2.3.6, with μ 1 = μ 2 = μ. One calculates that

G(q) + n T q H T (q)μ -1 = √ m 1 -μ 2 1 -μ -μ 1 (74)
so that the non-negativity condition (iii) cannot hold if μ > 0. This means that the analyzed mode does not exist for arbitrary pre-impact condition. However, one has

pn = m 1 -μ 2 (1 + e n,1 ) ẏ(t -) -μ(1 + e n,2 ) ẋ(t -) -μ(1 + e n,1 ) ẏ(t -) + (1 + e n,2 ) ẋ(t -) (75)
The kinetic consistency holds if μ max(

(1+e n,1 ) ẏ(t -) (1+e n,2 ) ẋ(t -) , (1+e n,2 ) ẋ(t -) (1+e n,1 ) ẏ(t -) ).
Motivated by these examples we can thus reformulate the problem of existence of this mode by stating that it exists only if this condition on the friction and the pre-impact velocity holds true.

Proposition 11 Let S -denote the largest set of pre-impact velocities such that pn = -(G(q)n T q H T (q)[ μ]) † (I + E nn ) qnorm (t -) is componentwise non-negative and solves (69) (a)-(d). Then the mode corresponding to v t,i (t + ) + e t,i v t,i (t -) > 0 at each impacting point exists for pre-impact velocities in S -if and only if the inequality:

H T T (q)M -1 (q)Ξ -1 (q) -E nn qnorm (t -) qtan (t -) + t T q H T (q)[ μ](G(q) -n T q H T (q)[ μ]) † (I + E nn ) qnorm (t -) > -E tt H T T (q)M -1 (q)Ξ -1 (q) qnorm (t -) qtan (t -) ( 76
)
is satisfied.

T L (t) = - 1 2 qT norm t -(I -E nn ) T M(q, μ)(I + E nn ) qnorm t - + 1 2 qT norm t -(I + E nn ) T E T (q)[ μ]H T T (q)t q D t (q)t T q H T (q)[ μ] × E(q)(I + E nn ) qnorm t - + qT tan t -D t (q)t T q H T (q)[ μ]E(q)(I + E nn ) qnorm t - (78)
The matrix in the second term is always positive semi-definite, so this term is always 0. Notice that one cannot transform the first quadratic term as was done to pass from [START_REF] Khulief | Modeling of impact in multibody systems: an overview[END_REF] to [START_REF] Kovecses | Finite and impulsive motion of constrained mechanical systems via Jourdain's principle: discrete and hybrid parameter models[END_REF] because M(q, μ) may not be symmetric.

Lemma 2 Suppose that the conditions of Lemma 1 all hold. Then:

E nn 2 1 + 2 E nn 2 < 1 M(q, μ) 2 1 ( M(q,μ)+M T (q,μ) 2 ) -1 2 ⇒ (I -E nn ) T M(q, μ)(I + E nn ) > 0 ( 79
)
Proof One has (I -E nn ) T M(q, μ)(I +E nn ) = M(q, μ)+H (q, μ, e n,i ), with H (q, μ, e n,i ) = -E T nn M(q, μ)E nn -E T nn M(q, μ) + M(q, μ)E nn . Consider Theorem 3, with M = M(q, μ) and A = M(q, μ) + H (q, μ, e n,i ). Using Proposition 9.3.5 in [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF] and the triangular inequality of norms one finds H (q, μ, e n,i ) 2

E nn 2 2 M(q, μ) 2 + 2 E nn 2 M(q, μ) 2 . Thus, it suffices that M(q, μ) + M T (q, μ) 2 -1 2 E nn 2 2 M(q, μ) 2 + 2 E nn 2 M(q, μ) 2 < 1
and the result follows.

Assume that the conditions of Lemmas 1 and 2 hold, and let us denote M(q, μ, e n,i ) = (I -E nn ) T M(q, μ)(I + E nn ) > 0 but not necessarily symmetric, and

K(q, μ, e n,i ) = (I + E nn ) T E T (q)[ μ]H T T (q)t q D t (q)t T q H T (q)[ μ]E(q)(I + E nn ) (80)
Notice that E nn = I implies that M(q, μ, e n,i ) = 0. One has:

T L (t) - 1 2 qT norm t -M(q, μ, e n,i ) qnorm t - + 1 2 qT norm t -K(q, μ i , e n,i ) qnorm t - + qT tan t -D t (q)t T q H T (q)[ μ]E(q)(I + E nn ) qnorm t - (81)
Theorem 1 Assume that the conditions of Lemmas 1 and 2 hold. If qnorm (t -) for some > 0, then there exist

δ 1 > 0, δ 2 > 0, δ 3 > 0, δ 4 > 0 such that if I + E nn 2 δ 1 , qtan (t -) δ 2 , D t (q)t T q H T (q) 2 δ 3 , [ μ] 2 δ 4 , then T L (t) 0.
Proof From (81), one has

T L (t) - 1 2 λ min M(q, μ, e n,i ) qnorm t -2 + 1 2 λ max K(q, μ i , e n,i ) qnorm t -2 + D t (q)t T q H T (q) 2 [ μ] 2 I + E nn 2 qnorm t - qtan t - (82)
Eigenvalues are continuous functions of the matrix entries. Choosing [μ] = 0 yields K(q, 0, e n,i ) = 0. Therefore, λ max (K(q, μ, e n,i )) can be made as small as desired by decreasing the coefficients μ i . On the other hand in view of the definitions of both E(q) and M(q, μ), λ min ( M(q, μ, e n,i )) can be made as close as possible to λ min (G -1 (q)) > 0 that does not depend on μ i . The result follows.

Roughly speaking, Theorem 1 says that if the frictional effects are small enough, then the impact under the considered mode is energetically consistent provided that the quasi-normal pre-impact velocity is sufficiently large while the quasi-tangential pre-impact velocity is sufficiently small. In other words, the impact, viewed in a generalized space, should be close to a "normal" collision. It is noteworthy that this is true even if the normal/tangential coupling term n T q H T (q) vanishes. The interest of Theorem 1 is that the various bounds can be calculated for a system with a given Delassus' matrix. This result is, however, conservative because it does not take into account the sign of the last term in [START_REF] Rockafellar | Convex Analysis[END_REF], but only its norm. It is also noteworthy that even if the coupling term n T q H T (q) vanishes, the energetic consistency is not guaranteed because of the quadratic term with K(q, μ i , e n,i ) 0. This is nevertheless the case even in the simplest case of a particle hitting a rough plane; see Example 5. Another way to formulate Theorem 1 is as follows:

Theorem 2 Assume that the conditions of Lemmas 1 and 2 hold. Then provided that

(i) λ min M(q, μ, e n,i ) > λ max K(q, μ i , e n,i ) (ii) qnorm (t -) qtan (t -) 2 D t (q)t T q H T (q) 2 [ μ] 2 I + E nn 2 λ min ( M(q, μ, e n,i )) -λ max (K(q, μ i , e n,i )) (83) one has T L (t) 0.
Proof Follows directly from [START_REF] Seghete | Conditions for uniqueness in simultaneous impact with application to mechanical design[END_REF].

Remark 8 It is clear that one can replace [ μ] in (68) by the matrix [ μ] = diag(μ i ξ i ∇h i (q) M -1 )
, where ξ i ∈ sgn(v t,i (t + ) + e t,i v t,i (t -)). The analysis of Theorem 1 remains valid, since

[ μ] 2

[ μ] 2 . This means that the energetic consistency criterion is valid for all the modes (whereas existence still has to be checked by a generalized equation like ( 68)), and this explains its conservativeness.

Example 5 Let us illustrate the above developments on the simplest case of a planar particle hitting a line. The horizontal position is x, the vertical one (normal to the line) is y. One

has qnorm = √ m ẏ, qtan = √ m ẋ, pn = 1 √ m p n , E(q) = 1, M(q, μ) = 1, M(q, μ, e n ) = 1 - e 2 n , K(q, μ i , e n,i ) = μ 2 (1 + e n ) 2 , D t = 1, t T q D t t q = 1 m , t T q H T (q) = 1 √ m , G(q) = 1, pn =
step of the analysis made in [START_REF] Génot | New results on Painlevé paradoxes[END_REF] about the characterization of the LCP matrix. Contrarily to some previous results which state the well-posedness of the contact LCP showing only the existence of small enough friction [START_REF] Matrosov | Right-hand solutions of the differential equations of dynamics for mechanical systems with sliding friction[END_REF][START_REF] Trinkle | On dynamic multi-rigid-body contact problems with Coulomb friction[END_REF] or for the all-sticking mode only [START_REF] Pang | Complementarity formulations and existence of solutions of dynamic multirigid-body contact problems with Coulomb friction[END_REF], here criteria with explicit bounds are given for general tangential behaviors in the multi-contact case. To the best of the author's knowledge, only [START_REF] Ivanov | Singularities in the dynamics of systems with non-ideal constraints[END_REF] has proposed explicit criteria to test the existence and uniqueness of solutions to the contact LCP when p = 2, with detailed conditions that guarantee that the contact LCP matrix is a P-matrix.

Single-contact systems

The analysis starts from the properties of the LCP in [START_REF] Duindam | Singularity-free dynamic equations of open-chain mechanisms with general holonomic and nonholonomic joints[END_REF] to determine the normal multiplier. It essentially relies on the properties of the matrix (here a scalar)

1n T q H T (q) μ ∇h(q) M -1 sgn H T T (q) q , [START_REF] Yen | An efficient Newton-type iteration for the numerical solution of highly oscillatory constrained multbody dynamic systems[END_REF] because G(q) = 1 when p = 1. For instance, in the case of the rod of Sect. 2.3.3, this scalar is equal to [1 + 3 1+3 cos 2 (θ) μ cos(θ ) sin(θ ) sgn(H T T (q) q)]. The different modes of the LCP for the sliding rod are easy to obtain in the scalar case, however, the subsequent analysis of the dynamics is rather intricate [START_REF] Génot | New results on Painlevé paradoxes[END_REF]. Indeed even if the rod is initialized in a well-posed mode, it may evolve towards regions of the (θ, θ) plane where singularities (non-uniqueness or unboundedness of the contact force) occur. This is in contrast with the analysis of impact conditions where q is considered as being constant. A central question in such analysis is: Does there exist a critical value μ c > 0 such that no singularity occurs for all 0 < μ < μ c ? This is the case for the system of Sect. 2.3.3 where μ c (θ ) = 1+3 cos 2 (θ) 3 sin(θ) cos(θ) 4 3 [START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF]. But this is not the case for the system of Sect. with minimum value at the configuration θ = arctan( √ m 1 /m 2 ) given by 2 m 1 m 2 , where m 1 and m 2 are the lumped masses of a rigid arm connected to a mass-spring-damper, respectively. It is thus possible to design these systems in a way such that singularities occur for arbitrarily small μ > 0 by properly choosing the masses m 1 and m 2 . In the scalar case, one nevertheless has the following result: Proposition 12 Let p = 1. Suppose that n T q H T (q)

1 ∇h(q) M -1 is a bounded quantity and that G(q) G min > 0 for any q in the admissible domain of C. Then there exists μ c > 0 such that for all 0 μ < μ c the LCP in [START_REF] Duindam | Singularity-free dynamic equations of open-chain mechanisms with general holonomic and nonholonomic joints[END_REF] has a unique bounded solution for any tangential relative velocity v t = H T T (q) q.

Proposition 12 applies to the classical Painlevé example of a rod sliding on a rough surface. Then μ c = 4 3 [START_REF] Génot | New results on Painlevé paradoxes[END_REF]. Let us examine the system in Sect. 2.3.4. Then we obtain n T q H T (q)

1 ∇h(q) M -1 = l l sin(ϕ)+r cos(ϕ) k 2 +(l sin(ϕ)+r cos(ϕ)) 2 . This quantity is bounded for any k 0 as long as ϕ ∈ [0, π 2 ]. One finds when k = 0 that (86) vanishes for μ c = sin(ϕ) + r l cos(ϕ) if sgn(H T T (q) q) = +1 (in [54, Appendix A] the value r l is found because the system is analyzed at ϕ = 0), and that μ < μ c guarantees the positivity of (86).

Multi-contact systems

The analysis is more intricate in the multi-contact case (i.e. p 2) for the simple reason that the LCP modes are much less easy to enumerate (we let apart the subsequent complex anal-ysis of the system's trajectories). Roughly speaking, the term -n T q H T (q) diag( μ i ξ i ∇h i (q) M -1 ) in ( 26) should be a "nice" perturbation of R p×p G(q) > 0, so that the LCP matrix is a P-matrix (the LCP (26) has a unique solution for any b(q, q, t)), or a strictly co-positive matrix (the LCP has at least a solution for any b(q, q, t) [START_REF] Cottle | In: The Linear Complementarity Problem[END_REF]Theorems 3.3.7 and 3.8.5]).

Example 6 Consider the system of Sect. 2.3.6. One has -n T q H T (q) diag μ i ξ i ∇h i (q) M -1 = -0 μ 2 sgn(v t,2 ) μ 1 sgn(v t,1 ) 0 and G(q) = I.

Notice that the perturbation needs not be symmetric, except if μ 1 = μ 2 and v t,1 = v t,2 , which occurs if ẋ = ẏ. A characterization of a P-matrix is that all its principal minors are positive [START_REF] Cottle | In: The Linear Complementarity Problem[END_REF]Definition 3.3.1]. The matrix of interest here is 1 -μ 2 sgn(v t,2 ) -μ 1 sgn(v t,1 ) 1

Its first principal minor is 1; the second one is 1μ 1 μ 2 sgn(v t,2 ) sgn(v t,1 ). So if the two tangential velocities have reversed signs, the matrix is always a P-matrix whatever the friction.

If the tangential velocities have the same signs, then it is a P-matrix if and only if μ 1 μ 2 < 1 (the strict inequality is important since otherwise one obtains a P 0 -matrix). This allows for a high friction at one contact, provided the other contact has low friction.

Apart from the basic result that if the LCP matrix is a P-matrix then the LCP always has a unique solution, some refined results may be stated, taking into account that G(q) = n T q M(q)n q is symmetric positive definite. Let us denote the perturbation matrix as P (q, μ, v t ) = -n T q H T (q) diag μ i ξ i ∇h i (q) M -1

, ξ i ∈ sgn(v t,i ) (

Proposition 13 Suppose that G(q) > 0 and P (q, μ, v t ) = A(q) with A(q) symmetric, real. Then there exists * such that for all | | < * the LCP matrix D(q, μ, v t ) is positive definite, and the LCP in [START_REF] Duindam | Singularity-free dynamic equations of open-chain mechanisms with general holonomic and nonholonomic joints[END_REF] has a unique solution for any b(q, q, t).

Proof As for Proposition 4, the result follows from [START_REF] Lancaster | The Theory of Matrices[END_REF]Exercise 8,p. 218].

The smallness of may be a consequence of either small friction, or of small couplings between normal and tangential directions, that is the matrix n T q H T (q) may have small entries. One difficulty with the analysis of this problem is that in general, the LCP matrix is not symmetric. Perturbations that preserve the positive definiteness without requiring symmetry can also be characterized as follows.

Proposition 14 Let G(q) > 0. Suppose that n T q H T (q) 2 ( p i=1

μ 2 i ξ 2 i ∇h i (q) 2 M -1 ) 1 2 < 1 G -1 (q) 2 .
Then the LCP matrix D(q, μ, v t ) = G(q) + P (q, μ, v t ) is positive definite, and the LCP in [START_REF] Duindam | Singularity-free dynamic equations of open-chain mechanisms with general holonomic and nonholonomic joints[END_REF] has a unique solution for any b(q, q, t).

Proof Let M = G(q), A = G(q)+P (q, μ, v t ) = D(q, μ, v t ), so that M -A = -P (q, μ, v t ). From [START_REF] Chen | Perturbation bounds of P-matrix linear complementarity problems[END_REF]Theorem 2.11] (see Theorem 3 in Appendix A), the condition G -1 (q) 2 × P (q, μ, v t ) 2 < 1 guarantees that D(q, μ, v t ) > 0. Now from Proposition 9.3.5 in [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF] one has P (q, μ, v t ) 2 n T q H T (q) 2 diag( μ i ξ i ∇h i (q) M -1 ) 2 ; the result follows.

Again the interest of Proposition 14 is that it provides explicit upper-bounds on both the friction and normal/tangential couplings, as a function of the kinetic angles. Another result using the co-positiveness is as follows:

Proposition 15 Let G(q) + P (q, μ, v t ) be co-positive. If b(q, q, t) ∈ Q * G+P , then the LCP in [START_REF] Duindam | Singularity-free dynamic equations of open-chain mechanisms with general holonomic and nonholonomic joints[END_REF] has a solution.

Proof Directly from [START_REF] Cottle | In: The Linear Complementarity Problem[END_REF]Theorem 3.8.6].

For small dimensions of the contact LCP (say, no more than 3 or 4 contact points), the co-positivity may be checked [START_REF] Hiriart Urruty | A variational approach to copositive matrices[END_REF]. Further necessary conditions may be derived from the fact that G(q) has diagonal entries equal to unity, but are not discussed here for the sake of briefness of the paper. The criteria derived in Propositions 6, 7, 9, 14, Lemmas 1, 2, Theorem 2, can be refined when dealing with particular systems.

Conclusions

In this paper, it is shown that a specific state-space transformation of the Lagrange dynamics with unilateral constraints and impacts, yields a suitable framework for the analysis of kinematic impact laws, and of the contact linear complementarity problem (closely related to the occurrence of the Painlevé paradoxes). The transformed velocity may also be seen as a quasi-velocity, constructed from projections in the kinetic metric. Kinematic impact laws with and without friction are analyzed, in particular the energetical consistency is carefully studied in all cases. Explicit, computable bounds on the friction coefficients and on the inertial couplings that guarantee the consistency are provided. Preliminary results on the well-posedness of the contact LCP in the multi-contact case are also given. Examples illustrate the developments.

1 2 .

 2 The co-positivity conditions of Proposition 2 for 2 × 2 matrices and diagonal E nn yield8 |e n,1 | 1, |e n,2 | 1, and8 Necessary and sufficient conditions for co-positivity of A = A T are a 11 0, a 22 0, a 12 + √ a 11 a 22 0[START_REF] Hiriart Urruty | A variational approach to copositive matrices[END_REF].

  2.3.4, as shown in [54, Appendix A]. The frictional impact oscillator shown in [54, Fig. 4.1] and analyzed therein, also has μ c = tan(θ ) + m 1

whose diagonal entries are equal to 1. It is positive definite if and only if n q has full rank m + p. Notice that we can split qnorm as qnorm = qb norm qu

	with qb norm ∈ R m corresponds to bilateral constraints, norm ∈ R p corresponds to unilateral constraints. Similarly, one has norm and qu

  + 2a 2 e n,1 e n,2 + 1a 2 0. Suppose that a = 0. Clearly, the second criterion holds for any |e n,1 | 1, |e n,2 | 1 and is less conservative than the first one which implies that

		2 n,1 -1)(e 2 n,2 -1) 0, which is equivalent to (1 -a 2 )e 2 n,1 e 2 n,2 -e 2 n,1 -
	e 2 n,2 e 2 n,1 + e 2 n,2	1 4 , which in turn implies that e n,1	1 2 and e n,2	1

Clearly, if the vectors t q,i are chosen mutually orthogonal then D t (q) = I .

This is equivalently stated as q(t -) ∈ -T Φ u (q(t))[START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF].

x T E T nn G -1 (q)x = x T (E T nn G -1 (q)) T x = x T G -1 (q)E nn x for any vector x ∈ R p .

equalities from (iii) and 2p 2 inequalities from (i) and (ii). Since p 2 -p (p +1) 2 > 0 for all p 2, there are more unknowns than equalities: the problem of calculating the coefficients of restitution (the entries of E nn ) from Corollary 1 seems to be under-determined, consequently the uniqueness of the coefficients

Recall that q is a constant during the shock, so that all the matrices that depend on q only are also constant. Thus, P and H may depend on q as well.

The idea of doing such a change of coordinates was first introduced in[START_REF] Brogliato | Absolute stability and the Lagrange-Dirichlet theorem with monotone multivalued mappings[END_REF] in the context of maximal monotone differential inclusions.

This leads us to the following result. Proposition 10 Assume that (i) E nn is non-negative, (ii) G(q)n T q H T (q)[ μ] is full-rank, (iii) (G(q)n T q H T (q)[ μ]) -1 is non-negative. Then the mode corresponding to v t,i (t + ) + e t,i v t,i (t -) > 0 at each impacting point exists if and only if the inequality:

-E nn qnorm (t -) qtan (t -) + t T q H T (q)[ μ](G(q)n T q H T (q)[ μ]) -1 (I + E nn ) qnorm (t -)

> -E tt H T T (q)M -1 (q)Ξ -1 (q) qnorm (t -) qtan (t -) [START_REF] Neimark | Dynamics of Non-holonomic Systems[END_REF] is satisfied.

If e t,i = 0 for all i, then the mode corresponds to a sliding mode. Let us recall that all functions are computed at q which is the impact position. Suppose that μ i = μ for all 1 i p , G(q) > 0 and that n T q H T (q) is a symmetric matrix. Then for μ small enough, the matrix G(q)n T q H T (q)[ μ] is symmetric positive definite (see the proof of Proposition 4) so (ii) is satisfied. If in addition it is a Z and a K-matrix, 11 then (iii) of the proposition is satisfied by [START_REF] Cottle | In: The Linear Complementarity Problem[END_REF]Theorem 3.11.10]. Such properties usually yield quite stringent conditions on the mechanical system under study. In fact, they may be over-conservative: examples show that solutions may exist even with relatively large μ. Another point is that not all the elements qnorm,i (t -) may be < 0, some of them may be equal to zero. Therefore, the nonnegativity condition (iii) for kinetic consistency may also be too conservative. This is the case for the rocking block as Example 3 shows.

A similar system as ( 70) can be constructed and examined for each mode. In practice, such an enumerative procedure is quite time-consuming and almost impossible for more than few impact points. For a single contact like a falling rod in Sect. 2.3.3, one can derive the conditions of existence of each mode [START_REF] Payr | Oblique frictional impact of a bar: analysis and comparison of different impact laws[END_REF]. What is more reasonable is therefore to prove that the generalized equation in [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques. Frictionless unilateral constraints and inelastic shocks[END_REF] possesses at least one solution, and then to use a numerical method to calculate it. In the simplest case of a planar particle hitting a rough ground with m = 1 kg, coordinates as in Fig. 1(d), and using (67) (with

, where proj denotes the orthogonal projection on [-1, 1]. One may extend this result to the inclusion [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques. Frictionless unilateral constraints and inelastic shocks[END_REF] to obtain a generalized equation with unknown qtan (t + ), and then deduce general conditions for its well-posedness. Propositions 10 and 11 are a preliminary analysis.

Example 3 Let us consider the system in Sect. 2.3.2, with μ 1 = μ 2 = μ. Starting from [START_REF] Génot | New results on Painlevé paradoxes[END_REF] one can compute that at a double impact with θ = 0 (noticing that ∇h i (q) q = √ f (0) for i = 1, 2):

11 Also called an M-matrix [25, p. 222].

Energetic consistency (Coulomb's friction)

Besides the existence of solution to systems like in [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques. Frictionless unilateral constraints and inelastic shocks[END_REF] or [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], it is necessary to examine the energetic consistency of the law that couples normal impact restitution with Coulomb's friction at the impulse level. It is known since the celebrated Kane's counter-example that such models may violate the dissipativity [START_REF] Glocker | Energetic consistency conditions for standard impacts[END_REF]. It is therefore of utmost importance to provide some conditions on the parameters that guarantee the dissipativity. The kinetic energy loss is expressed as follows at a time t of impact:

Let us focus on the mode of Proposition 11. Let us denote

Lemma 1 Assume that G(q) is positive definite. Then:

Proof (i) G(q) is symmetric positive definite. Applying Theorem 3 with M = G(q) and A = G(q)n T q H T (q)[ μ], one finds that the inequality in (i) guarantees that G(q)n T q H T (q)[ μ] is positive definite. Then this matrix has a positive definite inverse which is E(q). (ii) The proof follows from Corollary 3, with M = E(q), B = I -G(q)n T q H T (q)[ μ] and A = M(q, μ). Applying Proposition 9.3.5 in [START_REF] Bernstein | Matrix, Mathematics: Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF] to upper-bound G(q)n T q H T (q)[ μ] 2 by the product of norms, the result follows.

Given the Delassus matrix, the conditions of (i) can be computed as restrictions on both the normal/tangential couplings and the friction coefficients. Once these conditions hold, the conditions of (ii) can be calculated.

We now suppose that (i) of Lemma 1 holds and that pn can be calculated consistently from the first equation in [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques. Frictionless unilateral constraints and inelastic shocks[END_REF], i.e. pn = -E -1 (q)(I + E nn ) qnorm (t -) 0. Therefore, p t = -[ μ] pn in the considered mode. After some manipulations, T L (t) in (3.2.5) can be split into three parts (using in particular the second equation in [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques. Frictionless unilateral constraints and inelastic shocks[END_REF] to obtain qtan (t + ) + qtan (t -)):

with ξ ∈ sgn( ẋ(t + ) + e t ẋ(t -)). The conditions of the Theorems imply that e n 1, while the kinematic admissibility implies that e n 0. The direct application of Theorem 2 gives:

Notice that condition (i) implies that e n < 1. If e n = 0, then μ < 1 and

1 only, the frictionless case is admitted, because as noted above in that case M(q, μ, e n ) = 0, and we have excluded this case from the beginning. For this system, the energy loss is calculated to be:

The collision is dissipative if and only if ((e n -1)

Suppose that e t = 0, ẋ(t -) = 0, and the collision does not reverse the tangential velocity, so that |ξ | = 1 and ξ ẋ(t -) = | ẋ(t -)| 0. We can rewrite the dissipativity condition as

. Case (a) strictly contains [START_REF] Towne | One-dimensional collisions and Chebyshev polynomials[END_REF], because the sign of the third term in [START_REF] Seghete | Variational solutions to simultaneous collisions between multiple rigid bodies[END_REF] has been taken into account.

Glocker [START_REF] Glocker | Energetic consistency conditions for standard impacts[END_REF] analyses the energetic consistency of such kinematic impact laws for various tangential models (isotropic and orthotropic Coulomb friction) in particular cases: orthogonal constraints (all contacts decoupled), restricted coefficients of restitution (small, elastic impacts, plastic impacts). He also formulates friction in a way similar to [START_REF] Morarescu | Trajectory tracking control of multiconstraint complementarity Lagrangian systems[END_REF] in [START_REF] Glocker | An introduction to impacts[END_REF]. Leine and van de Wouw [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF] also derive general criteria, starting from an inclusion formalism for both the normal and the frictional parts. Lubarda [START_REF] Lubarda | The bounds on the coefficients of restitution for the frictional impact of rigid pendulum against a fixed surface[END_REF] derives bounds on the kinematic, kinetic, and energetic restitution coefficients for single collisions between a planar body and a rough fixed surface. Lankarani et al. [START_REF] Lankarani | Treatment of impact with friction in planar multibody mechanical systems[END_REF] use Routh's graphical method to analyze single impacts with Poisson's coefficient of restitution in multibody systems. Both articles use Coulomb's friction at the impulse level, as in ( 67) with e t,i = 0.

Analysis of the contact LCP (26)

In this section, it is assumed that p 1 contacts are lasting and are either sticking or sliding. Painlevé paradoxes occur in constrained systems with Coulomb's friction [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering[END_REF]. They represent some kind of singularities in the dynamics, which yield phenomena like nonuniqueness, non-existence, or unboundedness of the contact force [START_REF] Génot | New results on Painlevé paradoxes[END_REF]. The central tool for the analysis of the Painlevé paradoxes is the complementarity problem [START_REF] Duindam | Singularity-free dynamic equations of open-chain mechanisms with general holonomic and nonholonomic joints[END_REF]. One notices that the matrix of the LCP in [START_REF] Duindam | Singularity-free dynamic equations of open-chain mechanisms with general holonomic and nonholonomic joints[END_REF] is the same as the matrix that appears in the first line of [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques. Frictionless unilateral constraints and inelastic shocks[END_REF]. The main properties for the analysis of the LCP are the P-property, and co-positivity. Another major difference between Painlevé issues and the impact issues is that q varies in the former, so that the LCP matrix may have changing properties along the system's trajectory. In this section, we just point out some preliminary results that may be used to extend the first Appendix A: Theorem 2.11 in [START_REF] Chen | Perturbation bounds of P-matrix linear complementarity problems[END_REF] We give here just an excerpt of this result, and two easy corollaries of it.

Theorem 3 Let M ∈ R n×n be a positive definite matrix. Then every matrix Theorem 4 Let M ∈ R n×n be a co-positive matrix, and let q ∈ R n be given. If the implication v 0, Mv 0, v T Mv = 0 ⇒ v T q 0 is valid, then the LCP: 0 λ ⊥ Mλ + q 0 has a solution.

The implication is equivalently rewritten as q ∈ Q * M . Let us start from [START_REF] Liu | Another form of equations of motion for constrained multibody systems[END_REF]. Since R has full rank and using P Ḡ(q) = Γ T (⇔ Γ R -1 = ḠT R), this is equivalent to

Appendix

Let y = z + RΛv(t -). In view of the definition of W (q), and since N W (q) (•) = ∂ψ W (q) (•), one has N W (q) (y) = ∂ψ W (q) (y) = (Γ R -1 ) T ∂ψ W (q) (Γ R -1 y), where the chain rule of convex analysis has been used. The expression in ( 60) is then a consequence of the fact that xy ∈ -N K (x) ⇔ x = proj[K; y] for any vectors x, y, and convex set K of appropriate dimensions.