Bernard Brogliato

Inertial couplings between unilateral and bilateral holonomic constraints in frictionless Lagrangian systems

Keywords: Solid mechanics, Multibody systems, Unilateral constraints, Bilateral constraints, Linear complementarity problem, Lagrangian system, Delassus' matrix, Lagrange multiplier, Reduced coordinates, Kinetic angle, Gauss' principle, Impact, Singularity

In this paper, the following problem is analyzed: Given a frictionless Lagrangian system subject to complementarity relations (due to a set of unilateral constraints) that define a linear complementarity problem whose matrix is the so-called Delassus' matrix, study the influence of a set of bilateral constraints added to the dynamics on the Delassus' matrix. Two main paths are followed: the Lagrange multipliers method and the reduced coordinates method. The link with optimization (the Gauss' principle of mechanics) and the case of impacts, are also examined. The kinetic angles between the bilateral and the unilateral constraints are used to study the definiteness of the Delassus' matrix.

Introduction

The analysis, modeling, control, and simulation of mechanical systems subject to unilateral constraints have witnessed an intense research activity in the past two decades; see, for instance, [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Brogliato | Nonsmooth Mechanics[END_REF][START_REF] Frémond | Non-smooth Thermomechanics[END_REF][START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-smooth Systems[END_REF][START_REF] Leine | Dynamics and Bifurcations of Non-smooth Mechanical Systems[END_REF][START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF][START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF][START_REF] Stronge | Impact Mechanics[END_REF][START_REF] Studer | Numerics of Unilateral Contacts and Friction: Modeling and Numerical Time Integration in Non-smooth Dynamics[END_REF]. Historically, the issues related to the existence and uniqueness of contact forces, i.e., of the well-posedness of the unilateral contact problem stated with complementarity conditions, have been first analyzed rigorously in particular cases by Etienne Delassus in [START_REF] Delassus | Mémoire sur la théorie des liaisons finies unilatérales[END_REF]. Later Jean Jacques Moreau used convex analysis tools to study the unilateral contact problem in [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF], and proved that under a "compatibility" condition of the constraints (i.e., the constraints are linearly independent) the contact force multipliers and the generalized acceleration exist and are unique. Later the problem has been analyzed by Lötstedt in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF] who relaxed some assumptions and proved the local Bipop team-project, ZIRST Montbonnot, INRIA, 655 avenue de l'Europe, 38334 Saint-Ismier cedex, France e-mail: bernard.brogliato@inria.fr well-posedness of the dynamical system with persistent contact, using tools from complementarity theory (that were not available to Moreau at the time of writing of [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF]). See [START_REF] Brogliato | Nonsmooth Mechanics[END_REF]Sect. 5.1] for some details on this issue and the Delassus' example. Those works consider only what happens during phases of persistent contact with the constraints. Issues related to the existence and uniqueness of solutions to the dynamical system made of the Lagrangian dynamics plus the complementarity conditions and the impact dynamics have been tackled in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF][START_REF] Dzonou | A sweeping process approach to inelastic contact problems with general inertia operators[END_REF][START_REF] Dzonou | A convergence result for a vibro-impact problem with a general inertia operator[END_REF][START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF][START_REF] Monteiro Marques | Differential Inclusions in Non-smooth Mechanical Problems: Shocks and Dry Friction[END_REF][START_REF] Paoli | Time-stepping approximation of rigid-body dynamics with perfect unilateral constraints. I: The inelastic impact case. II: The partially elastic impact case[END_REF][START_REF] Paoli | A proximal-like method for a class of second order measure-differential inclusions describing vibro-impact problems[END_REF][START_REF] Paoli | Mouvement à nombre fini de degrés de liberté avec contraintes unilatérales: cas avec perte d'énergie[END_REF][START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF]. Most often in these works it is considered a unique differentiable unilateral constraint, exceptions being [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF][START_REF] Paoli | Time-stepping approximation of rigid-body dynamics with perfect unilateral constraints. I: The inelastic impact case. II: The partially elastic impact case[END_REF][START_REF] Paoli | A proximal-like method for a class of second order measure-differential inclusions describing vibro-impact problems[END_REF] (see Remark 2).

A central tool in the study of unilaterally constrained mechanical system is the linear complementarity problem for the contact forces, whose matrix is the so-called Delassus' matrix. The Delassus' matrix, whose accurate definition is recalled later in this paper, rules most of the properties of such nonsmooth dynamical systems: existence and uniqueness of solutions, existence and uniqueness of contact forces multipliers, stability issues [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF][START_REF] Leine | Stability properties of equilibrium sets of non-linear mechanical systems with dry friction and impact[END_REF], and numerical issues [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Studer | Numerics of Unilateral Contacts and Friction: Modeling and Numerical Time Integration in Non-smooth Dynamics[END_REF]. It is a fact that when the Delassus' matrix is not positive definite (due for instance to the presence of Coulomb's friction) the dynamics becomes much more complex [START_REF] Génot | New results on Painlevé paradoxes[END_REF][START_REF] Or | Investigation of Painlevé's paradox and dynamic jamming during mechanism sliding motion[END_REF]. In this paper, the following issues are tackled that concern the dynamics of the system when it evolves on the constraints boundaries:

• The case with only unilateral constraints is briefly examined, revisiting, and extending some results of Moreau [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF] and of Lötstedt [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. • The case when both unilateral and bilateral constraints act on the system is deeply examined, mainly using the Lagrange multipliers method. Various criteria are proposed to study the constrained Delassus' matrix, and the existence and uniqueness of the contact force multipliers. • The link with optimization is made, and various quadratic programs for the contact force multipliers and for the generalized acceleration are derived, thereby extending the results of Moreau [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF]. • In all the first three items, the solvability of the contact force problem is studied, in the case when the Delassus' matrix is positive semidefinite only. • Concerning impacts: the existence and uniqueness of the percussion vector is studied, quadratic programs for the percussion and the post-impact velocity are constructed, and the link with Moreau's impact law is made.

Singularities of bilaterally constrained systems are widely studied; see, e.g., [START_REF] Arponen | Regularization of constraint singularities in multibody systems[END_REF][START_REF] Blajer | Augmented Lagrangian formulation: geometrical interpretation and application to systems with singularities and redundancy[END_REF][START_REF] Müller | A Lie-group formulation of kinematics and dynamics of constrained MBS and its application to analytical mechanics[END_REF][START_REF] Müller | A conservative elimination procedure for permanently redundant closure constraints in MBSmodels with relative coordinates[END_REF][START_REF] Neto | Stabilization methods for the integration of DAE in the presence of redundant constraints[END_REF]. This work may also be seen as a study of the singularities that bilateral constraints may add in a well-posed unilaterally constrained Lagrangian system. The paper is organized as follows. In Sect. 2, the dynamics are introduced, and some results that concern the unilaterally constrained case are recalled or revisited. In Sect. 3, the notion of kinetic angles is recalled. Section 4 is dedicated to study the problem with the Lagrange multipliers method. The new linear complementarity problem for the unilateral/bilateral case is constructed, and its well-posedness is studied. The link with optimization and the Gauss' principle of mechanics is made, and various quadratic programs for the acceleration and the contact force multipliers are derived and studied. In Sect. 5, the issue is examined in a reduced coordinate framework. The issue of collisions is studied in Sect. [START_REF] Brogliato | Nonsmooth Mechanics[END_REF], with a generalized kinematic (Newton's like) impact law. Various optimization problems are constructed, and the link with Moreau's impact law is made. Conclusions end the paper in Sect. [START_REF] Cottle | On a problem in linear inequalities[END_REF]. In the Appendices, some theorems are recalled, which are used at several places in the paper.

Notations and definitions For f : R n → R p , ∇f (•) is its Euclidean gradient and ∇f (q) ∈ R n×p for all q ∈ R n . Its Jacobian is the transpose of its gradient: ∂f ∂q (q) = ∇f T (q). When needed, we shall sometimes denote ∇ q f (q). A > 0 means that A is positive definite, A ≥ 0 means that it is positive semidefinite. For two n × n positive definite matrices A and B, A > B means that A -B is positive definite. The product of two vectors x and y in the metric defined by a matrix A = A T > 0 is denoted as x, y A = x T Ay. Let x ∈ R n be a vector, x ≥ 0 (> 0) means that all entries x i ≥ 0 (> 0). N (A) denotes the kernel of the matrix A, R(A) denotes its image. I is the identity matrix with proper dimension. Let K ⊂ R n be a nonempty closed convex cone, N K (x) denotes its normal cone at x, T K (x) is its tangent cone at x. The projection of a vector x on K in the metric defined by a matrix A = A T > 0 on K is denoted as proj A [K; x]. A square matrix A is copositive if x T Ax ≥ 0 for all x ∈ R n + . An LCP(M, q) is a problem of the form 0 ≤ λ ⊥ w = Mλ + q ≥ 0, for a matrix M and a vector q of suitable dimensions. An LCP is said solvable if it has a solution. M is called a P-matrix if and only if the LCP(M, q) has a unique solution for any vector q. A positive definite matrix is a P-matrix. A Mixed LCP (MLCP) is a problem of the form 0 ≤ w ⊥ b + Dλ + Bw ≥ 0, a + Aλ + Cw = 0. Let S ⊂ R n be a nonempty set. Its dual set is S * = {z ∈ R n |z T y ≥ 0 for all y ∈ S}, and is always a convex cone. A matrix A is sometimes denoted as [a ij ]. The matrix diag(a ii ) denotes the diagonal matrix with ith entry a ii . A † is the Moore-Penrose generalized inverse of A.

Nomenclature Constraints acceleration terms: H 0 (q, q, t) in (3), H 1 (q, q) in ( 10), H 2 (q, q, t) in ( 12), H 3 (q, q, t) in ( 15), H 4 (q, q, t) in ( 21); constraints matrices: A(q) in Assumption 1, A c (q) in ( 16) and ( 18), M -1 c (q) in ( 17), G(q) in ( 18); n: number of degrees of freedom of the unconstrained system, p: number of bilateral holonomic constraints, m: number of unilateral constraints.

The dynamical system

Unilaterally constrained system

Let us consider a complementarity frictionless Lagrangian system with a configuration space C and a generalized coordinate vector q ∈ C ⊂ R n . Let us for the moment disregard the impacts that may occur in the system. In other words, it is assumed that the velocities are of class C 1 (I ; R n ), where I ⊂ R + is an interval, possibly open (in which case q(•) is supposed to be right continuous). Its dynamics are given by (the time argument is dropped for convenience):

⎧ ⎨ ⎩ M(q) q + F (q, q, t) = ∇h(q)λ u ,

0 ≤ λ u ⊥ h(q) ≥ 0, (1)
where h : R n → R m is a set of m unilateral constraints, λ u ∈ R m is a vector of Lagrange multipliers, M(q) ∈ R n×n is a positive definite inertia matrix, and F (q, q, t) contains nonlinear torques (Coriolis, centrifugal forces), forces that derive from a potential (gravity, elasticity), and possible external actions (like control inputs, or disturbances). The unilateral constraints define a finitely represented admissible domain Φ ⊂ C: Φ = {q ∈ C|h(q) ≥ 0}, and we assume that h(•) is of class C 1 (R n ; R m ) and Φ = ∅. The vector of generalized coordinates q is a minimal set of coordinates for the system in the interior of Φ, i.e., when h(q) > 0. Let us state the following:

Assumption 1 The so-called Delassus' matrix A(q) = ∇h T (q)M -1 (q)∇h(q) ∈ R m×m is positive definite for all q ∈ C.

Obviously, A(q) is symmetric, and Assumption 1 implies that m ≤ n. In many applications, this is not satisfied, however, our objective here is to study how bilateral constraints influence the Delassus' matrix, so assuming that A(q) > 0 is sufficient for us. The positive definiteness of A(q) holds if and only if the constraints h i (q) are linearly independent, equivalently ∇h(q) has rank m. Indeed N (A(q)) = N (M -1 2 (q)∇h(q)) (see Theorem 2.4.3 iii) [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF] in Appendix A.5), which in turn is equal to N (∇h(q)) since M -1 2 (q) has full rank n. Thus N (A(q)) = {0} is equivalent to N (∇h(q)) = {0}, that is, ∇h(q) has full rank m (which implies that m ≤ n).

Proposition 1

The following linear complementarity problem (LCP):

0 ≤ λ u ⊥ w = A(q)λ u + H 0 (q, q, t) ≥ 0 ( 2 )
with

H 0 (q, q, t) = -∇h T (q)M -1 (q)F (q, q, t) + d dt ∇h T (q) q, ( 3
)
has a unique solution whatever H 0 (q, q, t), if and only if Assumption 1 holds.

Proof The condition of Assumption 1 is necessary and sufficient for A(q) to be a P-matrix, since it is symmetric. Since the vector d dt (∇h T (q)) q may take arbitrary directions, the result follows from the fundamental result on existence and uniqueness of solutions of LCPs [START_REF] Cottle | The Linear Complementarity Problem[END_REF]Theorem 3.3.7].

The contact LCP in (2), whose unknown is the multiplier λ u , appears in the dynamics when the trajectories evolve smoothly on the boundary of the admissible domain bd(Φ), so that the complementarity conditions in (1) imply

0 ≤ λ u (t) ⊥ w(t) = d 2 dt 2 h q(t) = ∇h T q(t) q(t) + d dt ∇h T q(t) q(t) ≥ 0. (4)
Using [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF] and Assumption 1, it follows that the LCP in ( 2) is in turn equivalent to

0 ≤ ḧ(q) ⊥ y(t) = A -1 (q) ḧ(q) -A -1 (q)H 0 (q, q, t) ≥ 0, (5)
where ḧ(q) denotes d 2 dt 2 (h(q(t))), and we notice that w in (2) is equal to ḧ(q), while y in ( 5) is equal to λ u . Thus, the LCP in ( 5) is just the "inverse" of the LCP in [START_REF] Arponen | Regularization of constraint singularities in multibody systems[END_REF]. Theorem 1 in [START_REF] Cottle | On a problem in linear inequalities[END_REF] applies to [START_REF] Blajer | Augmented Lagrangian formulation: geometrical interpretation and application to systems with singularities and redundancy[END_REF] and we infer that both ḧ(q) and y(t) exist and are unique. Since λ u solution of ( 2) is a Lipschitz continuous function of H 0 (q, q, t), it follows that (1) is an ordinary differential equation and classical existence and uniqueness results may be used to study the local well-posedness of (1). In case Assumption 1 does not hold, A(q) may be only positive semidefinite. In such a case, one cannot construct the LCP in [START_REF] Blajer | Augmented Lagrangian formulation: geometrical interpretation and application to systems with singularities and redundancy[END_REF], and one has to work with (2), i.e., at the multiplier level. Then the following holds.

Proposition 2 Let m and n be arbitrary. If λ u,1 and λ u,2 are any two solutions of the LCP in [START_REF] Arponen | Regularization of constraint singularities in multibody systems[END_REF], then (i) ∇h T (q)λ u,1 = ∇h T (q)λ u,2 and (ii) H T 0 (q, q, t)λ u,1 = H T 0 (q, q, t)λ u,2 . If the LCP in (2) has at least one solution at t , and given unique (q(t), q(t)), then q(t) is unique.

Proof A(q) is positive semidefinite and symmetric, therefore, from Theorem 3.1.7(d) in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] (see Sect. A.1) one has A(q)λ u,1 = A(q)λ u,2 for any two solutions λ u,1 and λ u,2 . Let us denote M -1 2 (q) the symmetric square root of M -1 (q), and B(q) = M -1 2 (q)∇h(q). We have B T (q)B(q)(λ u,1λ u,2 ) = 0. If λ u,1λ u,2 ∈ N (B(q)), the result (i) follows. If not then necessarily B(q)(λ u,1λ u,2 ) ∈ N (B T (q)). Now from Proposition 1, page 199 in [START_REF] Lancaster | The Theory of Matrices[END_REF] one has N (B T (q)) = (R(B(q))) ⊥ . Since B(q)(λ u,1λ u,2 ) ∈ R(B(q)) and N (B T (q)) ∩ R(B(q)) = {0}, the inclusion holds only if B(q)(λ u,1λ u,2 ) = 0. Noting that M -1 2 (q) has full rank n (i) follows. The result (ii) follows from Theorem 3.1.7(a) in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] noting that A(q) is symmetric. Uniqueness of the acceleration follows from the dynamics. Proposition 2 therefore says that ∇h(q)λ u is unique, even if the Delassus' matrix is not full rank. Let us now state the following.

Proposition 3

Let m and n be arbitrary.

(i) The LCP in (2) is solvable if for any z ∈ R n in the set of solutions of the homogeneous LCP 0 ≤ z ⊥ A(q)z ≥ 0, one has z T H 0 (q, q, t) ≥ 0. (ii) Suppose that d dt (∇h(q(t))) = 0. Then the LCP in (2) is solvable.
Proof since A(q) is positive semidefinite, it is also copositive [8, p. 176]. (i) The solvability follows directly from Theorem 3.8.6 in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] (recalled in Sect. A.2). (ii) Let λ T A(q)λ = 0, then ∇h(q)λ = 0 and from the proposition assumption in item (ii) one has λ T H 0 (q, q, t) = -λ T ∇h T (q)M -1 (q)F (q, q, t) = 0. This shows that all the solutions of the homogeneous LCP 0 ≤ λ ⊥ A(q)λ ≥ 0, which in particular all satisfy the orthogonality condition λ T A(q)λ = 0, also satisfy λ T H 0 (q, q, t) ≥ 0. Thus, Theorem 3.8.6 in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] applies and one infers that the contact LCP in (2) has a solution.

It is interesting to note that the fact that ∇h(q(t)) is constant along the solution curve q(t) implies the solvability of the contact LCP in [START_REF] Arponen | Regularization of constraint singularities in multibody systems[END_REF]. In such a case, it follows from Proposition 2 that the acceleration q(t) at time t is unique, provided (q(t), q(t)) is unique. The homogeneous LCP 0 ≤ λ ⊥ A(q)λ ≥ 0 corresponds to the contact LCP of a fictitious system with the same mass matrix as the original system, where the inertial and exogeneous torques F (q, q, t) are orthogonal in the kinetic metric to the unilateral constraints (the kinetic metric is defined in Sect. 3), and ∇h(q(t)) is constant along the solution curve q(t). It is noteworthy that when Assumption 1 is not satisfied so that A(q) is positive semidefinite only, there is no reason that λ u exists or, in case existence holds, that it is unique. The LCP in (2) has been studied previously in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. Proposition 1 is proved in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF] and is recalled here only for convenience. In [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF], the constraints h i (q) ≥ 0 are assumed to be compatible, which we understand as functionally independent, and the uniqueness of λ u is proved from a quadratic program for the acceleration; see Sect. 4.4. Propositions 2 and 3 are, to the best of our knowledge, original results. The result of Proposition 2 is indirectly stated in the proof of Theorem 5.4 in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF], which proves the local existence and uniqueness of an analytic solution (q(•), q(•)) when the set of active constraints remains unchanged (see nevertheless Remark 1 below). The time functions in (4) are assumed to be right continuous, and are therefore to be understood as their right limits at each t . This allows us to study the system on an interval I = [t 1 , t 2 ), and to determine the values of the functions on the right of t 2 from the LCP. This LCP plays a fundamental role in numerical methods as well (see, e.g., [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]). Since an LCP can be rewritten as a quadratic problem, the above developments may be given the interpretation of the extension of Gauss' principle to systems with unilateral constraints, as done in [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF] (see Sect. [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]

.4 of this paper).

Example 1 To fix the ideas, let us consider a very simple example. Let a planar block whose mass center has coordinates q T = (xy) be subject to three unilateral constraints representing Fig. 1 The unilaterally constrained block the signed distance from three aligned points of its base, to a perfectly flat ground; see Fig. 1.

x is the vertical coordinate and y the horizontal one. The contact is frictionless. We assume that the block cannot rotate (it has two degrees of freedom in the plane). If the block's height is 2l, then this corresponds to three constraints h 1 (q) = h 2 (q) = h 3 (q) = xl ≥ 0. One has M(q) = diag(m) ∈ R 2×2 , F (q, q, t) = (mg 0) T + f (t) for some external action f (t) ∈ R 2 , λ u = (λ u,1 λ u,2 λ u,3 ) T , and

A(q) = 1 m ⎛ ⎝ 1 1 1 1 1 1 1 1 1 ⎞ ⎠ .
Suppose for instance that the block is initialized with the three points at contact and with f (•) = 0. Then h 1 (q) = h 2 (q) = h 3 (q) = 0 on some time interval. At the equilibrium, one has that ∇h(q)λ u = (mg 0) T + f (t) so that λ u,1 + λ u,2 + λ u,3 = mg on the considered time interval. At this stage, the problem is a statical problem and there is an infinity of contact multipliers that solve the static equilibrium: this is a simple example of a hyperstatic system. Consider now that f (t) becomes nonzero and takes values such that the block may detach from the ground. Then we have to consider a dynamical problem that is ruled by the LCP (2

): 0 ≤ λ u ⊥ Aλ u -(g + f 1 (t) m )(1 1 1) T ≥ 0.
One has rank(A) = 1, and Propositions 2 and 3 apply. The block has a unique acceleration, despite the contact multipliers cannot be determined in a unique way. See Remark 7 for more details on hyperstatic systems.

Unilaterally/bilaterally constrained system

Let us now suppose that in addition to the unilateral constraints, the system is subjected to a set of p bilateral holonomic constraints g(q) = 0, with g : R n → R p , and we assume that

g(•) is of class C 2 (R n ; R p ). The dynamics becomes ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩
M(q) q + F (q, q, t) = ∇h(q)λ u + ∇g(q)λ b ,

0 ≤ λ u ⊥ h(q) ≥ 0, g(q) = 0, (6)
where λ b ∈ R p is a vector of multipliers associated with the holonomic constraints. The system is therefore supposed to evolve in the intersection between the manifold

M = q ∈ C|g(q) = ∇g T (q) q = 0 (7)
and the admissible domain Φ. One sees that (4) is then rewritten as:

0 ≤ λ u ⊥ w(t) = ∇h T (q)M -1 (q)∇h(q)λ u + ∇h T (q)M -1 (q)∇g(q)λ b + H 0 (q, q, t) ≥ 0. ( 8
)

The kinetic angles

The kinetic angle between two unilateral constraints is a quantity which reflects the couplings that exist between the normal directions to the constraints boundaries and the inertial properties of the system. Kinetic angles are known to play an important role in frictionless multiple impacts [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF][START_REF] Brogliato | Nonsmooth Mechanics[END_REF][START_REF] Glocker | Concepts for modeling impacts without friction[END_REF][START_REF] Glocker | The geometry of Newton's cradle[END_REF][START_REF] Paoli | Continuous dependence on data for vibro-impact problems[END_REF]. In particular, it is known that the value π 2 is a critical value of the kinetic angle. Roughly speaking, trajectories are continuous with respect to initial data for kinetic angles ≤ π 2 , and discontinuous for kinetic angles > π 2 , depending on the restitution coefficients [START_REF] Paoli | Continuous dependence on data for vibro-impact problems[END_REF]. Kinetic angles are angles calculated in the kinetic metric, that is, the metric defined from the inertia matrix. The kinetic angle θ 12 (q) between two manifolds finitely represented as f 1 (q) = 0 and f 2 (q) = 0 is given by

θ 12 (q) = π -arccos ∇f 1 (q) T M -1 (q)∇f 2 (q) ∇f 1 (q) T M -1 (q)∇f 1 (q) ∇f 2 (q) T M -1 (q)∇f 2 (q) . ( 9
)
When two unilateral constraints are collided at the same time, their kinetic angle can be computed this way. One may also compute the kinetic angle between two bilateral constraints, or one unilateral constraints that is hit and a bilateral constraint. One checks from (9) that if ∇f 1 (q) T M -1 (q)∇f 2 (q) = 0 then θ 12 (q) = π 2 (the two constraints are orthogonal in the kinetic metric at the point q ∈ C). When the two constraints are linearly dependent, i.e., ∇f 1 (q) = α∇f 2 (q) for some nonzero α, then θ 12 (q) = π . The rationale behind [START_REF] Delassus | Mémoire sur la théorie des liaisons finies unilatérales[END_REF] is to start from the definition of the normal to a manifold in a Riemannian metric. Let f : R n → R be differentiable. At a point q ∈ C, the normal to a manifold {q ∈ R n |f (q) = 0} in the kinetic metric is the vector n(q) = M -1 (q)∇f (q) √ ∇f T (q)M -1 (q)∇f (q) . The scalar product between two such normal vectors is n 1 (q), n 2 (q

) M = n T 1 (q)M(q)n 2 (q) = ∇f T 1 (q)M -1 (q)∇f 2 (q) √ ∇f 1 (q) T M -1 (q)∇f 1 (q) √ ∇f 2 (q) T M -1 (q)∇f 2 (q) = ∇f 1 (q),∇f 2 (q) M -1 √ ∇f 1 (q) T M -1 (q)∇f 1 (q) √ ∇f 2 (q) T M -1 (q)∇f 2 (q)
. The π in (9) stems from the fact that the normal vectors point outwards the admissible domain Φ. The Delassus' matrix A(q) has a clear relation with the kinetic angles between the unilateral constraints. Indeed its entries are equal to ∇h T i (q)M -1 (q)∇h j (q), i, j ∈ {1, . . . , m}. If the unilateral constraints are mutually orthogonal then A(q) = diag(∇h T i (q)M -1 (q)∇h i (q)), i ∈ {1, . . . , m}. The constraints are then dynamically decoupled at the impact times.

The Lagrange multipliers method

4.1 The unilateral/bilateral LCP Since the system evolves in M, it follows that d 2 dt 2 (g(q(t))) = 0 for all t ≥ 0. One has

d dt g q(t) = ∇g T q(t) q(t) and d dt ∇g T q(t) q(t) = ∇g T q(t) q(t) + d dt ∇g T q(t) q(t).
For simplicity, we shall denote d 2 dt 2 (g(q(t))) = ∇g T (q(t)) q(t) + H 1 (q(t), q(t)), with

H 1 (q, q) = d dt ∇g T (q) q. ( 10
)
Since d 2 dt 2 (g(q(t))) = 0 for all t ≥ 0, it follows using ( 6) that

∇g T (q)M -1 (q)∇g(q)λ b + ∇g T (q)M -1 (q)∇h(q)λ u + H 2 (q, q, t) = 0, (11)
where

H 2 (q, q, t) = -∇g T (q)M -1 (q)F (q, q, t) + H 1 (q, q). (12)
Let us make the following:

Assumption 2

The p × p matrix ∇g T (q)M -1 (q)∇g(q) has full rank p for all q ∈ C, that is, it is positive definite.

Similarly to Assumption 1, Assumption 2 means that the functions g i (•) are mutually independent, it implies in particular that p ≤ n and that rank(∇g(q)) = p for all q ∈ C. Thus,

λ b = -∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 (q)∇h(q)λ u -∇g T (q)M -1 (q)∇g(q) -1 H 2 (q, q, t). ( 13
)
Inserting ( 13) into ( 6), we obtain the constrained dynamics:

M(q) q + F (q, q, t) = ∇h(q) -∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 (q) × ∇h(q) λ u + H 3 (q, q, t), ( 14
)
where

H 3 (q, q, t) = -∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 H 2 (q, q, t). ( 15
)
The constrained Delassus' matrix is therefore given by A c (q) = ∇h T (q)M -1 (q)∇h(q) -∇h T (q)M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q)

-1 × ∇g T (q)M -1 (q)∇h(q). ( 16
)
Thus, this corresponds to changing M -1 (q) to a new "constrained" inverse mass matrix in the Delassus' matrix:

M -1 c (q) = M -1 (q) -M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 (q) (17)
so that the new Delassus' matrix is A c (q) = ∇h T (q)M -1 c (q)∇h(q). The notation M -1 c (q) does not mean that this is the inverse of some matrix but is chosen just to mimic the unconstrained case. Let us denote G(q) = ∇g(q) ∇g T (q)M -1 (q)∇g(q)

-1 ∇g T (q), [START_REF] Glocker | The geometry of Newton's cradle[END_REF] so that

A c (q) = ∇h T (q)M -1 (q) I -G(q)M -1 (q) ∇h(q). (19)
This will be used later.

Remark 1 Notice that ( 8) and ( 11) form a Mixed LCP [START_REF] Cottle | The Linear Complementarity Problem[END_REF] that is transformed into an LCP using Assumption 2. The set of equations in [START_REF] Dzonou | A convergence result for a vibro-impact problem with a general inertia operator[END_REF] and w(t) = 0 with w(t) in ( 8), does not fit in general with Eqs. (5.2a) and (5.2b) in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF], because of the presence of the terms d dt (∇g T (q)) q and d dt (∇h T (q)) q. Thus, Lemma 5.2 in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF] cannot be used in general, except if both the gradients of h(q) and g(q) are constant in time along the solution curve q(t) (as in Proposition 4 below). Moreover, an incorrect use of Theorem 1 in [START_REF] Cottle | On a problem in linear inequalities[END_REF] is made in the proof of [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]Lemma 5.2] to study the LCP in Eq. (5.5) of [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. The proof should rather follow the lines of the proofs of Propositions 2 and 3, using Theorem 3.1.7 in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] for LCP with positive semidefinite matrices.

During persistently unilaterally constrained phases of motion, the multiplier λ u is consequently computed as the solution of the LCP:

0 ≤ λ u ⊥ w = A c (q)λ u + H 4 (q, q, t) ≥ 0 ( 20
)
with

H 4 (q, q, t) = ∇h T (q)M -1 (q)H 3 (q, q, t) + H 0 (q, q, t). ( 21
)
The discrepancy between A c (q) and A(q) (and between H 0 (q, q, t) and H 4 (q, q, t)) is directly related to the inertial coupling matrix ∇h T (q)M -1 (q)∇g(q). Using similar calculations as in Sect. 2 to construct the LCP (8), we may deduce from ( 14) and ( 20) an LCP that is equivalent to [START_REF] Lancaster | The Theory of Matrices[END_REF], in case A c (q) has full rank m:

0 ≤ y = A -1 c (q) ḧ(q) -A -1 c (q)H 4 (q, q, t) ⊥ ḧ(q) ≥ 0. ( 22
)
Remark 2 Assumptions 1 and 2 are made in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF], Assumption 1 is made in [START_REF] Paoli | Time-stepping approximation of rigid-body dynamics with perfect unilateral constraints. I: The inelastic impact case. II: The partially elastic impact case[END_REF][START_REF] Paoli | A proximal-like method for a class of second order measure-differential inclusions describing vibro-impact problems[END_REF]. In most of the other works dealing with the well-posedness of (1) (without any bilateral constraints), only the case m = 1 is analyzed, so that Assumption 1 is trivially satisfied. Notice that bilateral constraints and unilateral constraints are analyzed separately in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF].

Analytical study of the unilateral/bilateral LCP

We first focus on the general case, when A c (q) is positive semidefinite only. Then conditions that guarantee A c (q) > 0 are examined.

Positive semi-definite A c (q)

Let us first state a result that can be deduced from [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF], where neither Assumption 1 nor Assumption 2 are supposed to hold. Thus, there are no assumptions on the number of bilateral or unilateral constraints and on the number of degrees of freedom.

Proposition 4 Let m, p and n be arbitrary. Assume that d dt (∇g(q(t))) = 0 and d dt (∇h(q(t))) = 0. (i) If λ u,1 , λ u,2 , λ b,1 , λ b,2 are any solutions of the problem in (8) and (11) at time t , ∇h(q)λ u,1 + ∇g(q)λ b,1 = ∇h(q)λ u,2 + ∇g(q)λ b,2 . Consequently, the acceleration q(t) at time t is unique, given unique (q(t), q(t)). (ii) Moreover, if the matrix ∇h T (q) ∇g T (q) M -1 (q)(∇h(q) ∇g(q)) ∈ R (m+p)×(m+p) is full rank, then the problem possesses at most one solution λ u and λ b .

Proof (i) Under the conditions of the proposition, the MLCP in ( 8) and ( 11) exactly fits with the MLCP in Eqs. (5.2a), (5.2b), (5.2c) in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. The first result follows from Lemma 5.2 in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF], taking into account the proof modification as pointed out in Remark 1. Taking this into account the fact that the acceleration is unique follows directly from the system's dynamics in [START_REF] Brogliato | Nonsmooth Mechanics[END_REF].

(ii) The uniqueness of the solutions (when they exist) of the problem in ( 8) and ( 11) at time t also follows from Lemma 5.2 in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF] (taking once again into account the proof modification).

From now on, it is no longer supposed that the functions h(q) and g(q) are constant along the solution curve q(t). Relaxing this assumption is important because as pointed out previously this leads to LCP(M, q) with a vector q that may take arbitrary directions. Then the P property of matrices (that is equivalent in our context to the positive definiteness because of the symmetry) is required to guarantee the existence and uniqueness of solutions. As Propositions 1, 2, and 3 show, when the LCP matrix is positive semidefinite only, the mere existence of solutions (without uniqueness) requires some care.

Proposition 5 Let Assumption 2 hold. Suppose that all the bilateral constraints are orthogonal in the kinetic metric to all the unilateral constraints. Then A c (q) = A(q) and H 4 (q, q, t) = H 0 (q, q, t), so that the LCP in [START_REF] Arponen | Regularization of constraint singularities in multibody systems[END_REF] and the LCP in [START_REF] Lancaster | The Theory of Matrices[END_REF] are identical.

Proof The conditions of the proposition imply that ∇g T (q)M -1 (q)∇h(q) = 0; the proof is complete from ( 12), ( 15), [START_REF] Glocker | Concepts for modeling impacts without friction[END_REF], and [START_REF] Leine | Dynamics and Bifurcations of Non-smooth Mechanical Systems[END_REF].

It is clear from [START_REF] Glocker | Concepts for modeling impacts without friction[END_REF] that if the inertial couplings between unilateral and bilateral constraints are small enough, then A c (q) and A(q) should possess the same definiteness (in particular A(q) > 0 ⇒ A c (q) > 0). This will be examined in Proposition 9 and Corollary 4. Let us now characterize precisely the matrix M -1 c (q) in [START_REF] Glocker | An introduction to impacts[END_REF]. It is crucial to note that the constrained Delassus' matrix A c (q) in ( 16) may be positive definite despite M -1 c (q) is not. This is due to the fact that M -1 c (q) ∈ R n×n and A c (q) ∈ R m×m do not have the same dimensions in general, except if m = n. For instance, 1 0 0 0 is positive semidefinite, but

(1 0) 1 0 0 0 1 0 = 1 is positive definite. Lemma 1 Let Assumption 2 hold. The kernel of the matrix M -1 c (q) in (17) is given by N M -1 c (q) = R ∇g(q) . ( 23
)
Its range is given by

R M -1 c = N ∇g T (q) . ( 24
)
Proof From ( 17), it follows that x ∈ N (M -1 c (q)) if and only if y ∈ N I -M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q) and y = M -1 (q)x.

Now since I -M -1 (q)∇g(q)[∇g T (q)M -1 (q)∇g(q)] -1 ∇g T (q) and M -1 (q)∇g(q) [∇g T (q)M -1 (q)∇g(q)] -1 ∇g T (q) are idempotent matrices (and are projectors since they are symmetric [4,

N I -M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q) = R M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q) .
Therefore, the kernel of M -1 c (q) is characterized by the range of M -1 (q)∇g(q) [∇g T (q)M -1 (q)∇g(q)] -1 ∇g T (q). Now we have

R M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q) = z ∈ R n |∃x ∈ R n such that z = M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)x . (25)
Thus,

N M -1 c (q) = M(q)y|y ∈ N I -M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q) = M(q)y|y ∈ R M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q) = M(q)y|y = M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)x, x ∈ R n = z ∈ R n |z = ∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)x, x ∈ R n = R ∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q) = R ∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 2 = R ∇g(q) , ( 26
)
where the last two equalities are obtained using Theorem 2.4.3 in [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF] and [20, Exercise 7, p. 78]. The proof for the range is similar. One has:

R M -1 c (q) = z ∈ R n |∃x ∈ R n such that z = M -1 (q)x -M -1 (q)∇g(q) × ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 (q)x = z ∈ R n |∃y ∈ R n such that z = y -M -1 (q)∇g(q)
× ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)y = R I -M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q)

-1 ∇g T (q)

= N M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)

= N ∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)

= N ∇g T (q)M -1 (q)∇g(q)

-1 2 ∇g T (q) = N ∇g T (q) , ( 27
)
using again Theorem 2.4.3 [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF] (see Appendix A.5) and [START_REF] Lancaster | The Theory of Matrices[END_REF]Exercise 7,p. 78]. Above

[∇g T (q)M -1 (q)∇g(q)] -1 2 is the positive definite square root of [∇g T (q)M -1 (q) ∇g(q)] -1 . Lemma 2 Let Assumption 2 hold. The kernel of M -1 c (q) is a linear subspace of R n of dimension p. Its range is a linear subspace of dimension n -p, and R(M -1 c (q)) = N (M -1 c (q)) ⊥ . Thus, M -1 c (q) is symmetric of rank n -p.
Proof The proof follows from Lemma 1 since R(∇g(q)) = N (∇g T (q)) ⊥ .

Lemma 3

Let Assumption 2 hold. The constrained inverse mass matrix M -1 c (q) is positive semidefinite.

Proof From (17) and 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] (see Appendix A.8), it follows that M -1 c (q) is positive semidefinite if and only if

I -M -1 2 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 2 (q) ≥ 0, ( 28
)
where M -1 2 (q) is the symmetric positive definite square root of M -1 (q). The matrix

M -1 2 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 2 (q)
is idempotent, and it is an orthogonal projector onto R(∇g(q)) [20, Corollary 1, p. 430]. Since R(∇g(q)) is a linear subspace containing {0}, the orthogonal projection onto it is nonexpansive [19, Proposition 3.1.3] and one has for any x ∈ R n :

x T M -1 2 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 2 (q)x ≤ x T x. ( 29
)
The result follows.

The matrix M -1 c (q) is therefore symmetric positive semi definite, of rank np. An easy consequence of Lemma 3 is that A c (q) is at least positive semidefinite under Assumption 2.

Let us now study extensions of the results of Sect. 2. Proposition 6 Let Assumption 2 hold. Let λ u,1 and λ u,2 be two solutions of the LCP in [START_REF] Lancaster | The Theory of Matrices[END_REF].

Then (i) ∇h(q)(λ u,1 -λ u,2 ) ∈ N (M -1 c (q)), (ii) (λ T u,1 -λ T u,2 )H 4 (q, q, t) = 0.
Proof (i) From Theorem 3.1.7 (d) in [START_REF] Cottle | The Linear Complementarity Problem[END_REF], it follows that A c (q)(λ u,1λ u,2 ) = 0. Since A c (q) is symmetric and at least positive semidefinite, it has a square root B(q) [20, p. 181] that is also symmetric positive semidefinite. Redoing a similar reasoning as in the proof of Proposition 2, we find that B(q)(λ u,1λ u,2 ) = 0. Now we can write B(q) as B(q) = P (q)M -1 2 (q)∇h(q), with P (q) = I -M -1 2 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 2 (q), and P (q)P T (q) = P (q) since this is a projector [4, Fact 3.9.1]. Thus, we get that P (q)M -1 2 (q)∇h(q)(λ u,1λ u,2 ) = 0. From the expression of P (q) and since M(q) is full rank, it follows that this is equivalent to M -1 c (q)∇h(q)(λ u,1λ u,2 ) = 0, from which (i) follows.

(ii) From Theorem 3.1.7(a) in [START_REF] Cottle | The Linear Complementarity Problem[END_REF], it follows that λ T u,1 (H 4 (q, q, t) + A c (q)λ u,2 ) = λ T u,2 (H 4 (q, q, t) + A c (q)λ u,1 ). Since A c (q) is symmetric the result follows.

If there are no bilateral constraints or if the conditions of Proposition 5 hold, then the results of Proposition 2 are recovered. A consequence of Proposition 6 is as follows.

Corollary 1 Let R(∇h(q)) ∩ R(∇g(q)) = {0} for all q ∈ C. Let λ u,1 and λ u,2 be two solutions of the LCP in [START_REF] Lancaster | The Theory of Matrices[END_REF]. Then ∇h(q)(λ u,1λ u,2 ) = 0.

Proof From Lemma 1 one has N (M -1 c (q)) = R(∇g(q)) {0}. Since ∇h(q)(λ u,1λ u,2 ) ∈ R(∇h(q)) from Proposition 6(i) the result follows.

Corollary 1 shows that it is sufficient that the unilateral and the bilateral constraints be linearly independent, to recover the result of Proposition 2(i). Let us now state a result that concerns the existence of solutions of the LCP in [START_REF] Lancaster | The Theory of Matrices[END_REF].

Lemma 4 Let Assumption 2 hold. (i)

The LCP in [START_REF] Lancaster | The Theory of Matrices[END_REF] is solvable if for any z ∈ R n in the set of solutions of the homogeneous LCP 0 ≤ z ⊥ A c (q)z ≥ 0 one has z T H 4 (q, q, t) ≥ 0; (ii) if d dt (∇g(q(t))) = 0 and d dt (∇h(q(t))) = 0 and R(∇h(q)) ∩ R(∇g(q)) = {0} for all q ∈ C, the LCP in (20) is solvable. Proof (i) A c (q) is positive semidefinite, so it is copositive. From Theorem 3.8.6 in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] the solvability follows.

(ii) The solutions of the homogeneous LCP 0 ≤ z ⊥ A c (q)z ≥ 0 satisfy ∇h(q)z ∈ N (M -1 c (q)). From Lemma 1, one has N (M -1 c (q)) = R(∇g(q)). Since ∇h(q)z ∈ R(∇h(q)), the conditions of (ii) imply that ∇h(q)z = 0. Now the conditions of (ii) also imply that H 4 (q, q, t) = ∇h T (q)M -1 (q)H 3 (q, q, t) -∇h T (q)M -1 (q)F (q, q, t), see ( 21) and (3). Therefore, z T H 4 (q, q, t) = 0, so the implication of (i) is satisfied and from Theorem 3.8.6 in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] the solvability follows.

Remark 3 Considering constant gradients is a restrictive assumption. Lemma 4(ii) is coherent with Proposition 4. The condition of Lemma 4(ii) can also be stated as H 4 (q, q, t) ∈

Q * Ac(q) , with Q Ac(q) = {z ∈ R m |0 ≤ z ⊥ A c (q)z ≥ 0} and Q * Ac(q)
is its dual set. It is easily seen that Q Ac(q) is a nonempty (since {0} belongs to it) convex cone, since A c (q) is symmetric.

Positive definite A c (q)

The previous results concern positive semi definite A c (q). It is of interest to study the conditions under which this matrix is positive definite, i.e., the addition of bilateral constraints does not destroy the positive definiteness of the LCP matrix.

Proposition 7 Let p = m = 1, then A c (q) is symmetric positive definite, except if the kinetic angle θ u 1 b 1 between the two constraints is θ u 1 b 1 = 0.
Proof In this case, [∇g T (q)M -1 (q)∇g(q)] -1 and ∇h T (q)M -1 ∇h(q) are nonzero scalars, so we can rewrite [START_REF] Glocker | Concepts for modeling impacts without friction[END_REF] as

A c (q) = ∇h T (q)M -1 ∇h(q) 1 - ∇g(q), ∇h(q) 2 M -1
∇h T (q)M -1 ∇h(q)∇g T (q)M -1 (q)∇g(q)

= ∇h T (q)M -1 ∇h(q) 1 -cos 2 (θ u 1 b 1 ) . ( 30
)
Using ( 9) and noting that cos(πθ) =cos(θ ) for any θ the equality in [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF] follows. The proof is complete by noting that the cosine belongs to

[-1, 1].
Obviously, under the conditions of Proposition 7, the matrix A c (q) is a positive scalar. The following is a first generalization of Proposition 7.

Proposition 8 Let m ≥ 1, p ≥ 1 and let Assumption 2 hold. A necessary condition for A c (q) to be positive definite is that there does not exist i ∈ {1, . . . , m} and j ∈ {1, . . . , p} such that the vector ∇h i (q) is colinear to ∇g j (q), for all q ∈ C.

Proof From lemma, if ∇h i (q) = α∇g j (q) for some i and some j , then A c (q) ∈ R m×m has its ith row and column which are both zero. Thus, it is of rank strictly less than m.

The above results (as well as Corollary 1) suggest that when the bilateral and the unilateral constraints are linearly independent, then the addition of bilateral constraints does not modify the properties of the contact force LCP. We can now state generalizations of the above results. The first result characterizes A c (q) from the kinetic angles between the constraints, and extends Propositions 5, 7, and 8.

Theorem 1

Let Assumptions 1 and 2 hold. Suppose that the unilateral constraints h i (•), i ∈ {1, . . . , m} are mutually orthogonal in the kinetic metric, and that the bilateral constraints g i (•), i ∈ {1, . . . , p} are also mutually orthogonal in the kinetic metric. Then A c (q) is positive definite if and only if

I > p i=1
cos θ ji (q) cos θ ki (q) , j,k ∈ {1, . . . , m}, where θ ij (q) is the kinetic angle between the constraint h i (•) and the constraint g j (•) at q ∈ C and I is the m × m identity matrix.

Proof It follows from [START_REF] Glocker | Concepts for modeling impacts without friction[END_REF] and from [4, Proposition 8.1.2(xii), (xiv)] that A c (q) > 0 if and only if I > ∇h T (q)M -1 (q)∇h(q) -1 2 ∇h T (q)M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 2 × ∇g T (q)M -1 (q)∇g(q) -1 2 ∇g T (q)M -1 (q)∇h(q) ∇h T (q)M -1 (q)∇h(q)

-1 2 . ( 31
)
This is obtained from ( 16) that is equivalent (by definition) to ∇h T (q)M -1 (q)∇h(q) > ∇h T (q)M -1 (q)∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 × ∇g T (q)M -1 (q)∇h(q) [START_REF] Müller | A Lie-group formulation of kinematics and dynamics of constrained MBS and its application to analytical mechanics[END_REF] and multiplying on the right and on the left both sides of (32) by the symmetric positive definite square root (∇h T (q)M -1 (q)∇h(q)) -1 2 . From the theorem conditions, one has (∇h T (q)M -1 (q)∇h(q)) -

1 2 = diag( 1 √ ∇h T i (q)M -1 (q)∇h i (q)
), and (∇g T (q)M -1 (q)∇g(q)) -

1 2 = diag( 1 √ ∇g T i (q)M -1 (q)∇g i (q)
). The first three matrices in the right-hand side of ( 31) define an m × p matrix C(q) with entries ∇f i (q) T M -1 (q)∇g j (q)

∇f i (q) T M -1 (q)∇f i (q) ∇g j (q) T M -1 (q)∇g j (q) . (

) 33
The term in ( 33) is equal to cos(πθ ij (q)) =cos(θ ij (q)), see [START_REF] Delassus | Mémoire sur la théorie des liaisons finies unilatérales[END_REF]. The right-hand side in ( 31) is equal to C(q)C T (q) and its entries (row j , column k) are the terms p i=1 cos(θ ji (q)) cos(θ ki (q)), j, k ∈ {1, . . . , m}. The proof is complete.

Clearly if the conditions of Proposition 5 are satisfied then this criterion yields I > 0. If m = p = 1, one sees that if ∇h(q) = α∇g(q) for some nonzero real α, then from (31) one gets I > I which is impossible, hence recovering the result of Proposition 8. In such a case, the orthogonality conditions of the theorem are trivially satisfied. If all the unilateral and the bilateral constraints are linearly dependent, then all the kinetic angles satisfy θ ij = π and the theorem condition states that I > [p], where [p] is the matrix whose all entries are equal to p. Then for all p ≥ 2 one has the first entry of I -[p] that is equal to 1p < 0 and this matrix cannot be positive definite. It is noteworthy that the theorem's conditions imply that p ≤ n and m ≤ n.

Corollary 2 Consider the conditions of Theorem 1. Then A c (q) is positive definite only if for each j ∈ {1, . . . , m} one has

p i=1 cos 2 θ ji (q) < 1. ( 34
)
Proof The j th diagonal element of the m × m matrix [ p i=1 cos(θ ji (q)) cos(θ ki (q))] is given by p i=1 cos 2 (θ ji (q)). The necessary condition thus follows from Theorem 1, since a positive definite matrix necessarily has all its diagonal elements positive.

It easily follows from Corollary 2 (see also Sect. 3) that if there exists one h i (•) that is colinear to one g j (•), then A c (q) is positive semidefinite only, and we recover the result of Proposition 8. From [START_REF] Neto | Stabilization methods for the integration of DAE in the presence of redundant constraints[END_REF], it follows that if p increases, then the kinetic angles between the unilateral and the bilateral constraints should approach π 2 in order to keep the positive definiteness of A c (q). Corollary 3 Consider the conditions of Theorem 1. Then A c (q) is positive definite if [START_REF] Neto | Stabilization methods for the integration of DAE in the presence of redundant constraints[END_REF] is satisfied and for each j ∈ {1, . . . , m} one has

m i=1,i =j p k=1 cos(θ jk ) cos(θ ik ) < 1 - p i=1 cos 2 (θ ji ). ( 35
)
Proof The conditions in [START_REF] Or | Investigation of Painlevé's paradox and dynamic jamming during mechanism sliding motion[END_REF] state the diagonal dominance of the matrix

I - p i=1
cos θ ji (q) cos θ ki (q) , j,k ∈ {1, . . . , m} [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Fact 4.10.15], where we notice that the right-hand side in [START_REF] Or | Investigation of Painlevé's paradox and dynamic jamming during mechanism sliding motion[END_REF] is positive due to [START_REF] Neto | Stabilization methods for the integration of DAE in the presence of redundant constraints[END_REF].

The condition [START_REF] Neto | Stabilization methods for the integration of DAE in the presence of redundant constraints[END_REF] assures that this matrix has all its diagonal elements positive. Hence, this matrix is positive definite [20, p. 373] and so is A c (q) from Theorem 1.

Notice that since we already know that A c (q) is positive semidefinite, only the strict inequalities in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], [START_REF] Neto | Stabilization methods for the integration of DAE in the presence of redundant constraints[END_REF], and ( 35) are of interest. Proposition 5 on one hand, and Theorem 1 and its two Corollaries 2 and 3 on the other hand, are complementary: the first one uses the orthogonality (in the kinetic metric) between the unilateral and the bilateral constraints, whereas the second ones uses the orthogonality between the unilateral constraints themselves and the orthogonality between the bilateral constraints themselves. Let us now relax the orthogonality conditions. To this aim, let us first make an assumption. Assumption 3 Let Assumptions 1 and 2 hold. Then ∇h T (q)M -1 (q)∇h(q)

-1 2 = D -1 2 h (q) + h (q)M h (q)
and ∇g T (q)M -1 (q)∇g(q)

-1 = D -1 g (q) + g (q)M g (q), where D h = diag(∇h T i (q)M -1 (q)∇h i (q)), D g = diag(∇g T i (q)M -1 (q)∇g i (q)), M h (q), and M g (q) have diagonal entries equal to zero and are symmetric matrices, h (q) ≥ 0 and g (q) ≥ 0 are real numbers.

In Assumption 3, the matrices h (q)M h (q) ∈ R m×m and g (q)M g (q) ∈ R p×p stem from the inertial couplings between the constraints, which are the off-diagonal entries of the matrices ∇h T (q)M -1 (q)∇h(q) and ∇g T (q)M -1 (q)∇g(q). In Appendix A.3, it is shown on a simple example that Assumption 3 is plausible. Proposition 9 Let Assumption 3 hold. Then A c (q) > 0 for small enough h (q) and g (q) if

I > D -1 2 h (q) ∇h T (q)M -1 (q)∇g(q) D -1 g (q) ∇g T (q)M -1 (q)∇h(q) D -1 2 h (q). ( 36
)
Proof The inequality in ( 31) is equivalent to A c (q) > 0. Let us denote E hg (q) = ∇h T (q)M -1 (q)∇g(q). After calculations, one finds ∇h T (q)M -1 (q)∇h(q)

-1

2 E hg (q) ∇g T (q)M -1 (q)∇g(q)

-1 × E T hg (q) ∇h T (q)M -1 (q)∇h(q)

-1 2 = D -1 2 h (q)E hg (q)D -1 g (q)E T hg (q)D -1 2 h (q) + O( h , g ), (37)
where

O( h , g ) = D -1 2 h (q)E hg (q) g (q)M g (q)E T hg (q)D -1 2 h (q) + D -1 2 h (q)E hg (q) D -1 g (q) + g (q)M g (q) E T hg (q) h (q)M h (q) + h (q)M h (q)E hg (q) D -1 g (q) + g (q)M g (q) E T hg (q)D -1 2 h (q) + h (q)M h (q)E hg (q) D -1 g (q) + g (q)M g (q) E T hg (q) h (q)M h (q). ( 38
)
We notice that O( h , g ) ∈ R m×m is symmetric. Hence, from [20, Exercise 8, p. 218] (see Appendix A.9), it follows that for small enough h (q) and g (q) the matrix I -

D -1 2 h (q)E hg (q)D -1 g (q)E T hg (q)D -1 2 h (q)+O( h , g ) is positive definite provided I -D -1 2 h (q) × E hg (q)D -1 g (q)E T hg (q)D -1 2 h (q) is positive definite. The result follows.
This proposition says that if the inertial couplings between the unilateral constraints h i (q), i ∈ {1, . . . , m}, and the inertial couplings between the bilateral constraints g i (q), i ∈ {1, . . . , p}, are small enough, then the test for the positive definiteness of A c (q) is the same as for Theorem 1 and Corollaries 2 and 3. The next result completes Proposition 5 and comes after Theorem 1. It concerns the couplings between the unilateral and the bilateral constraints.

Corollary 4 Let Assumption 3 hold. Suppose that E hg (q) = hg (q)I for some hg (q) ≥ 0. Then for small enough hg (q) one has A c (q) > 0.

Proof We rewrite [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] as

I > 2 hg (q) ∇h T (q)M -1 (q)∇h(q) -1 2 ∇g T (q)M -1 (q)∇g(q) -1 × ∇h T (q)M -1 (q)∇h(q) -1 2 .
Using again [20, Exercise 8, p. 218] (see Appendix A.9), the result follows.

The next result provides necessary and sufficient conditions for positive definiteness, from a different point of view than Theorem 1.

Proposition 10 Let Assumption 2 hold and let M

-1 2 c (q) be the unique symmetric positive semidefinite square root of M -1 c (q). The matrix A c (q) is positive definite if and only if the vectors M -1 2 c (q)∇h i (q), i ∈ {1, . . . , m}, are independent vectors of R n .

Proof From Lemma 3, M -1 c (q) ≥ 0. We have A c (q) = D T (q)D(q) where D(q) = M -1 2 c (q)∇h(q) ∈ R n×m . Thus, A c (q) > 0 ⇔ rank(D(q)) = m [20, p. 180], i.e., the m vectors M -1 2 c (q)∇h i (q) of R n are independent. Proposition 10 is general but it is not very constructive because the square root of M -1 c (q) has to be known, and it does not explicitly relate the bilateral and the unilateral constraints. Still another necessary and sufficient characterization of the positive definiteness of A c (q) from ( 31) is that the maximum eigenvalue of the matrix in the right-hand side of (31) be < 1 [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Lemma 8.4.1]. This, however, requires the calculation of the maximum eigenvalue of this matrix, and does not provide us with a convenient way to check that A c (q) > 0 when designing a mechanism.

Corollary 5 Let Assumption 2 hold. Then

A c (q) > 0 only if m + p ≤ n.
Proof Let D(q) be as in Proposition 10. We have that rank(D(q)) ≤ min[rank(M -1 2 c (q)), rank(∇h(q))] (see [20, p. 97]). Since rank(M -1 2 c (q)) = np (see Lemma 2 and Theorem 1, p. 181 [START_REF] Lancaster | The Theory of Matrices[END_REF]), if np < m then necessarily rank(D(q)) < m which contradicts Proposition 10. So necessarily np ≥ m.

Remark 4

In practice, one often has m > n, so the conditions of Corollary 5 are not met. This result is coherent with Proposition 5 in the sense that the conditions of both Proposition 5 and of Corollary 5 imply that np ≥ m.

Example

Let us consider a planar homogeneous slender rod of length 2l > 0, mass m > 0, with three degrees of freedom q T = (x y θ ), x and y being the center of mass coordinates and θ the orientation angle. The moment of inertia at the gravity center is I = ml 2 3 . The unilateral constraint is given by h(q) = xl cos(θ ) ≥ 0 and corresponds to the signed distance between a vertical wall and the slender's tip, where we assume that the x-axis points outwards the wall. The bilateral constraint is given by g(q) = yl sin(θ ) -H = 0 for some H > 0 and corresponds to constraining the other tip to a horizontal motion. The system is depicted in Fig. 2. After some calculations, one obtains

M -1 c (q) = ⎛ ⎜ ⎝ 1 m 0 0 0 3 cos 2 (θ) m(1+3 cos 2 (θ)) 3 cos(θ) ml(1+3 cos 2 (θ)) 0 3 cos(θ) ml(1+3 cos 2 (θ)) 3 ml 2 (1+3 cos 2 (θ)) ⎞ ⎟ ⎠ . ( 39
)
The second and third columns are colinear vectors of R 3 , and rank(M -1 c (q)) = 2. One also finds A c (q) = 1 m + 3 sin 2 (θ) m(1+3 cos 2 (θ)) > 0. Given the above unilateral constraint, the class of bilateral constraints which are orthogonal to it have the gradient ∇g T (q) = (-3 sin(θ) l γ β γ ), β, γ ∈ R. Such bilateral constraints are of the form g(q) = βy + δ = 0 for some δ ∈ R. This means that the center of gravity is constrained to move horizontally. In this simple case, the Euclidean and the kinetic orthogonality are the same.

Geometrical and optimization interpretations

In this section, we review several equivalent formulations of the LCPs in ( 2) and ( 20), strongly inspired by the results of Moreau [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF] and Lötstedt [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. In particular, it is shown that the bilateral constraints yield a new quadratic program for the acceleration, that is, the direct extension of the quadratic program derived in [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF]. Gauss' principle of Mechanics is therefore extended to systems with mixed unilateral/bilateral constraints. We first review the results for the purely unilaterally constraint case.

The unilaterally constrained case

From a classical result of convex analysis, the LCP in ( 2) is equivalent to the following inclusion:

A(q)λ u + H 0 (q, q, t) ∈ -N K (λ u ) [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF] with K = (R + ) m . When A(q) > 0 this is in turn equivalent to [12, p. 85] (recall that A(q) = A T (q)):

λ u = proj A(q) K; -A -1 (q)H 0 (q, q, t) . ( 41
)
It is also a classical fact [START_REF] Cottle | The Linear Complementarity Problem[END_REF] that the LCP ( 2) is equivalent to the quadratic program:

min λu≥0 1 2 λ T u A(q)λ u + H 0 (q, q, t) T λ u ( 42
)
and it is noteworthy that this equivalence holds also when A(q) is symmetric positive semidefinite only (in the positive definite case the solution is unique so the optimum is global [1, Theorem 12.36(a)]). Equivalence means that if λ * u is a solution of (42) then it is a solution of (2), and vice versa.

Proposition 11

Let m and n be arbitrary. Consider the quadratic program in [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF]. (i) If H 0 (q, q, t) ∈ R(A(q)), or (ii) if d dt (∇h(q(t))) = 0, then it has a solution.

Proof (i) From Exercise 2.10.25(b) in [START_REF] Cottle | The Linear Complementarity Problem[END_REF], it follows that the quadratic function of ( 42) is bounded below on the set {λ u ∈ R m |λ u ≥ 0}. Using Theorem 2.8.1 in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] the result follows.

(ii) One has R(A(q)) = R(∇h T (q)M -1 (q)∇h(q)) = R(∇h T (q)) (using Exercise 6, p. 180 and Exercise 7, p. 78 in [START_REF] Lancaster | The Theory of Matrices[END_REF]). So the proposition's condition guarantees that H 0 (q, q, t) ∈ R(A(q)) and (i) applies.

Proposition 11 is similar to Proposition 3, starting from a quadratic program rather than from an LCP. From the LCP in [START_REF] Blajer | Augmented Lagrangian formulation: geometrical interpretation and application to systems with singularities and redundancy[END_REF], it follows that ḧ(q) = proj A -1 (q) K; H 0 (q, q, t) [START_REF] Stronge | Impact Mechanics[END_REF] that is equivalent to min ḧ(q)≥0 1 2 ḧT (q)A -1 (q) ḧ(q) -ḧT (q)A -1 (q)H 0 (q, q, t) [START_REF] Studer | Numerics of Unilateral Contacts and Friction: Modeling and Numerical Time Integration in Non-smooth Dynamics[END_REF] and to

ḧ(q) = argmin z≥0 1 2 z -H 0 (q, q, t) T A -1 (q) z -H 0 (q, q, t) . ( 45
)
One may in turn formulate quadratic programs for the acceleration as follows. Starting from [START_REF] Studer | Numerics of Unilateral Contacts and Friction: Modeling and Numerical Time Integration in Non-smooth Dynamics[END_REF] and using the expression of ḧ(q) in ( 4), one finds that q is the solution of ⎧ ⎨ ⎩ min 1 2 qT ∇h(q)A -1 (q)∇h T (q) q + qT ∇h(q)A -1 (q)∇h T (q)M -1 (q)F (q, q, t)

subject to: ∇h T (q) q + d dt (∇h T (q)) q ≥ 0.

(

) 46
And starting from [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF], one obtains that q is the solution of ⎧ ⎨ ⎩ min 1 2 qT M(q) q + qT d dt (∇h(q))A -1 (q)∇h T (q) q subject to: A -1 (q)(∇h T (q) q + d dt (∇h T (q)) q) ≥ 0.

(47)

As shown in [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF] and later in [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF], the acceleration is also the solution of ⎧ ⎨ ⎩ min 1 2 ( q + M -1 (q)F (q, q, t)) T M(q)( q + M -1 (q)F (q, q, t))

subject to: ∇h T (q) q + d dt (∇h T (q)) q ≥ 0 (48) which may be obtained from Dorn's dual problem of ( 42) that is given by (see Appendix A.4):

⎧ ⎨ ⎩ min 1 2 λ T u A(q)λ u subject to: A(q)λ u + H 0 (q, q, t) ≥ 0 ( 49
)
and using (1) to replace ∇h(q)λ u . Dorn's duality theorem says that if λ * * u solves (49) then there exists λ * u that solves [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF] with [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF], then it solves (49) and the two extrema are equal. It follows that q is the unique solution of (48) even if A(q) is not full rank. We recover here the results of Propositions 1 and 2. Similarly, starting from (44) and using Dorn's duality, one finds that the acceleration is the solution of ⎧ ⎨ ⎩ min 1 2 [ qT ∇h(q) + qT d dt (∇h(q))]A -1 (q)[∇h T (q) q + d dt (∇h T (q)) q] subject to: A -1 (q)∇h T (q)( q + M -1 (q)F (q, q, t)) ≥ 0.

A(q)(λ * u -λ * * u ) = 0. In case A(q) > 0, then λ * u = λ * * u . Conversely, if λ * u solves
(50)

From a numerical point of view, solving the programs in ( 46), (47), and (50) requires to invert two matrices: M(q) and A(q), whereas the one in (48) involves only the inverse of M(q). Moreover, (48) holds even if A(q) is semidefinite, whereas (46), (47), and (50) require A(q) invertible. The results of this section are summarized as follows:

Solvability condition:

H 0 (q, q, t) ∈ R(A(q)) ⇓ LCP (2) λ * u solution ⇐⇒ QP (42) λ * u solution =⇒ if A(q) > 0 ⇐= ⎧ ⎨ ⎩ QP (49) λ * * u solution A(q)(λ * u -λ * * u ) = 0 ⇑ (if A(q) > 0)
Solvability condition:

H 0 (q, q, t) ∈ Q * A(q) (43) ⇔ (44) ⇔ (46) ⇔ (47) ⇔ (40) ⇔ (41) QP (48) unique solution q (51)
where Q * A(q) is defined in Remark 3. In [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF], Moreau derived the quadratic programs (48) and (49) using convex analysis tools and under Assumption 1, whereas Lötstedt [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]Sect. 6] used Dorn's duality to derive (49) from ( 42), and then (48), with no rank assumption on A(q). The issue of the solvability for the existence of a multiplier λ u was not tackled previously.

Mixed case: unilateral/bilateral constraints

Similarly the unilateral/bilateral case yields in case A c (q) > 0:

λ u = proj Ac(q) K; -A -1 c (q)H 4 (q, q, t) . ( 52
)
The addition of bilateral constraints therefore modifies the metric and the vector that is projected, but does not change the cone on which the projection is made. In the general case where A c (q) ≥ 0, the LCP ( 20) is equivalent to the quadratic program:

min λu≥0 1 2 λ T u A c (q)λ u + H T 4 (q, q, t)λ u . ( 53
)
The extension of Proposition 11 is as follows.

Proposition 12 Let m, p, and n be given.

(i) If H 4 (q, q, t) ∈ R(A c (q)), or (ii) if d dt (∇h(q(t)))
= 0 and all the unilateral constraints are orthogonal in the kinetic metric to all the bilateral constraints, or (iii) if d dt (∇h(q(t))) = 0, d dt (∇g(q(t))) = 0 and F (q, q, t) ∈ R(∇h(q)), then the quadratic program in (53) has a solution.

Proof (i) The proof is the same as for Proposition 11(i).

(ii) Follows from Proposition 11(ii) and using Proposition 5.

(iii) If the gradients are constant, then H 4 (q, q, t) = ∇h T (q)M -1 (q)[H 3 (q, q, t) -F (q, q, t)]. It follows from ( 12) and ( 15) that H 3 (q, q, t) -F (q, q, t) = -I -G(q)M -1 (q) F (q, q, t).

Therefore, under the proposition's condition and using [START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF] one has H 4 (q, q, t) = A c (q)x for some vector x.

Under the conditions of Proposition 12(iii), the LCP in (20) becomes 0 ≤ λ u ⊥ w = A c (q)λ u + A c (q)x ≥ 0 for some x ∈ R m . However, Theorem 1 in [START_REF] Cottle | On a problem in linear inequalities[END_REF] does not apply in the general case where A c (q) ≥ 0 since M -1 c (q) is positive semidefinite. The conditions of Proposition 12(iii) mean that the nonlinear and external forces have to be normal (in the Euclidean metric) to the unilateral constraints boundaries. This is much more restrictive than the conditions of Lemma 4, and of low practical interest for the design of mechanisms. The dual program of (53) is by Dorn's duality:

⎧ ⎨ ⎩ min 1 2 λ T u A c (q)λ u subject to: A c (q)λ u + H 4 (q, q, t) ≥ 0. (54)
In case A c (q) > 0, the two programs are equivalent. From Proposition 5, if all the bilateral constraints are orthogonal in the kinetic metric to all the unilateral ones, then both quadratic programs in (49) and in (54) are the same. From ( 14), one has ∇h T (q) q + ∇h T (q)M -1 (q)F (q, q, t) = A c (q)λ u + ∇h T (q)M -1 (q)H 3 (q, q, t). Let A c (q) be positive definite. Inserting this into (54), one has that the acceleration is the solution of:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ min 1 2 {∇h T (q) q + ∇h T (q)[M -1 (q)F -H 3 ]} T A -1 c (q){∇h T (q) q + ∇h T (q)[M -1 (q)F -H 3 ]}
subject to: ∇h T (q) q + ∇h T (q)M -1 (q)F (q, q, t) -∇h T (q)H 3 (q, q, t) ≥ 0.

(

) 55
When all the bilateral constraints are orthogonal to all the unilateral constraints, (55) reduces to (46). The next result is interesting because it presents the direct extension of (48) to the case with bilateral constraints, and holds for A c (q) ≥ 0.

Proposition 13 Consider the dynamics (6) and let Assumption 2 hold. Suppose that the LCP (20) is solvable. Then the acceleration is the solution of the quadratic program:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ min q 1 2 [ q + M -1 (q)[F (q, q, t) -H 3 (q, q, t)] T M(q)[ q + M -1 (q)[F (q, q, t) -H 3 (q, q, t)]
subject to: ∇h T (q)M -1 c (q)M(q)[ q + M -1 (q)(F (q, q, t) -H 3 (q, q, t))]

+ H 4 (q, q, t) ≥ 0. ( 56
)
Proof From ( 14), ( 16), [START_REF] Glocker | An introduction to impacts[END_REF], and ( 18), one has M(q) q + F (q, q, t) = I -G(q)M -1 (q) ∇h(q)λ u + H 3 (q, q, t).

(57)

The matrix G(q)M -1 (q) is a projector (see [START_REF] Glocker | The geometry of Newton's cradle[END_REF] for the definition of G(q)), and N (I -G(q)M -1 (q)) = R(G(q)M -1 (q)) = R(∇g(q)) (using similar arguments as in the proofs of Lemma 1 and Lemma 3). Thus, (I -G(q)M -1 (q)) ∈ R n×n is not full rank, but from the solvability assumption the existence of ∇h(q)λ u implies that M(q) q + F (q, q, t) -H 3 (q, q, t) ∈ R(I -G(q)M -1 (q)). Thus, (57) can be inverted and one has ∇h(q)λ u = I -G(q)M -1 (q)

† M(q) q + F (q, q, t) -H 3 (q, q, t) , (

where (I -G(q)M -1 (q)) † is the Moore-Penrose inverse of I -G(q)M -1 (q), and one has (I -G(q)M -1 (q)) † = I -G(q)M -1 (q) [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Fact 6.3.13]. Inserting (58) into (54) one obtains the program (the arguments of H 3 (q, q, t) and F (q, q, t) are dropped):

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ min 1 2 (M(q) q + F -H 3 ) T (I -G(q)M -1 (q)) T M -1 (q)(I -G(q)M -1 (q))(M(q) q + F -H 3 )
subject to: ∇h T (q)M -1 c (q)(I -G(q)M -1 (q))(M(q) q + F -H 3 ) + H 4 (q, q, t) ≥ 0.

(59) Now we have the following:

G(q)M -1 (q) is a projector so (G(q)M -1 (q)) † = G(q)M -1 (q) and N (G(q)M -1 (q)) = N (M -1 (q)G T (q)) = N (G(q)) = N (∇g T (q)) = R(∇g(q)) ⊥ . Further, M(q) q +F -H 3 ∈ R(I -G(q)M -1 (q)) ⊂ N (∇g T (q)) = R(∇g(q)) ⊥ .
Therefore, M(q) q +F -H 3 ∈ N (G(q)M -1 (q)). We infer that (I -G(q)M -1 (q))(M(q) q + F -H 3 ) = M(q) q + F -H 3 , and (59) simplifies to (56). When all the bilateral constraints are orthogonal to all the unilateral ones, one has ∇h T (q)M -1 c (q)M(q) = ∇h T (q), ∇h T (q)M -1 c (q)H 3 (q, q, t) = 0 (see ( 15) and ( 17)), and H 4 (q, q, t) = H 0 (q, q, t) (see [START_REF] Leine | Dynamics and Bifurcations of Non-smooth Mechanical Systems[END_REF], and Proposition 1 for the definition of H 0 (q, q, t)). It can be checked that in this case the quadratic programs in (56) and in (48) are identical. Let us now examine the feasible set of the quadratic program (56).

Proposition 14 Let us denote

F = z ∈ R n |∇h T (q)M -1 c (q)M(q)z + ∇h T (q)M -1 c (q) F (q, q, t) -H 3 (q, q, t) + H 4 (q, q, t) ≥ 0 (60)
the feasible set of the quadratic program in (56). (i) Then F = ∅ for any d dt (∇h T (q)) q and any d dt (∇g T (q)) q, only if there is no i ∈ {1, . . . , m} such that ∇h i (q) ∈ R(∇g(q)). (ii) If F is not empty, then the quadratic program in (56) has a unique solution.

Proof It is verified that M -1 c (q)M(q) is a projector, so that using [20, Theorem 1, p. 194] R(M -1 c (q)M(q)) = N (M -1 (q)∇g(q)[∇g T (q)M -1 (q)∇g(q)] -1 ∇g T (q)), and from Lemma 1, R(M -1 c (q)M(q)) = N (∇g T (q)). Similarly, N (M -1 c (q)M(q)) = R(∇g(q)). So M -1 c (q)M(q) is a projector onto N (∇g T (q)) = R(∇g(q)) ⊥ parallel to R(∇g(q)). Consequently, ∇h T (q)M -1 c (q)M(q)z = ∇h T (q)x with x ∈ R(∇g(q)) ⊥ . Suppose that there exists an index i such that ∇h i (q) ∈ R(∇g(q)). Then the ith entry of the m-vector ∇h T (q)M -1 c (q)M(q)z is zero. The ith line of the inequality in (56) thus boils down to H 4,i (q, q, t) ≥ 0. In view of ( 21), ( 15), [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF], and the definition of H 0 (q, q, t) in Proposition 1, one can choose a vector d dt (∇h T (q)) q such that H 4,i (q, q, t) < 0. So (i) is proved. (ii) Follows from the fact that M(q) = M T (q) > 0 for any q ∈ C, which means that the quadratic program in (56) is strictly convex. This result is coherent with Propositions 7 and 8.

Remark 5

The program in (55) has simpler constraints than the one in (56), but a more complex quadratic term which implies the inverse of a matrix. Moreover, it holds only in case A c (q) is full rank, which is not the case of (56) that holds independently of A c (q) being or not full rank. This should not come as a surprise since we know that the multiplier may not be unique while the acceleration is; see Proposition 2 and Corollary 1.

Let us now relate the feasibility of the program (56) with the solvability conditions of Lemma 4, where we use the notation of Remark 3.

Lemma 5 Suppose that the conditions of Lemma 4(i) hold, i.e., H 4 (q, q, t) ∈ Q * Ac(q) . Then the quadratic program (56) is feasible, i.e., the set F in (60) is non empty and the acceleration is unique.

Proof Consider a z as in Lemma 4(i). Following the proof of Lemma 4, one sees that ∇h(q)z ∈

N (M -1 c (q)) = R(∇g(q)). Let S = {z ∈ R m |∇h(q)z ∈ R(∇g(q))}. S is nonempty, since {0} ∈ R(∇h(q)) ∩ R(∇g(q)). Thus, z ∈ Q Ac(q) is equivalent to z ∈ S ∩ R m + . We deduce that H 4 (q, q, t) ∈ Q * Ac(q) = (S ∩ R m + ) * = S * + R m + because
S is an affine and convex cone, and [41, Corollary 16.4.2] is used. Thus, we can write H 4 (q, q, t) = x + y with x ∈ S * and y ∈ R m + . The set S * is given by S * = {y ∈ R m |y = ∇h T (q)x, x ∈ R(∇g(q)) * }. Thus H 4 (q, q, t) = x + y with x = ∇h T (q)α, α ∈ R(∇g(q)) * , and y ∈ R m + . Now notice that ∇h T (q)M -1 c (q)M(q)z + ∇h T (q)M -1 c (q)(F (q, q, t) -H 3 (q, q, t)) ∈ R(∇h T (q)). Therefore, the inequality in the definition of F in (60) is equivalent to ∇h T (q)(x + α) + y ≥ 0 with x ∈ R n , α ∈ R(∇g(q)) * and y ≥ 0. Whatever the value of α, one can always pick an x such that the inequality is satisfied. The proof is complete using Proposition 14.

Obviously Lemma 5 can also be proved indirectly by a sequence of equivalences and implications of the above LCPs and quadratic programs. The direct proof is interesting per se, however. It is worth reminding that Assumption 2 is supposed to hold for Lemma 5, because the quadratic program in (56) is constructed from (59) which needs that ∇g T (q)M -1 (q)∇g(q) be full rank p. Let us now briefly come back on the case without bilateral constraints.

Corollary 6

Suppose that the conditions of Proposition 3 hold, i.e., H 0 (q, q, t) ∈ Q * A(q) . Then the quadratic program in (48) is feasible, i.e., the set {z ∈ R n |∇h T (q)z + d dt (∇h T (q)) q ≥ 0} is nonempty. Then the acceleration is unique.

Proof The proof follows from Lemma 5 noting that the program in (56) reduces to the program (48) when there are no bilateral constraints (one may take g(q) = 0 so that the orthogonality conditions of Proposition 5 are trivially satisfied). The uniqueness of q follows since M(q) > 0.

The results of this section are summarized as follows: projection (52) QP (56) unique solution q

(if A c (q) > 0) LCP (20) λ * u solution ⇐⇒ QP (53) λ * u solution =⇒ if A c (q) > 0 ⇐= ⎧ ⎨ ⎩ QP (54) λ * * u solution A c (q)(λ * u -λ * * u ) = 0 ⇑ (if A c (q) > 0)
Solvability condition:

H 4 (q, q, t) ∈ Q * Ac(q)
QP (55) unique solution q

(61)

The reduced coordinates method

An alternative method to analyze a system with bilateral holonomic constraints is to reduce its generalized coordinate vector q, exploiting the fact that the system lives on the manifold M in [START_REF] Cottle | On a problem in linear inequalities[END_REF]. Hereafter, we shall use the McClamroch-Wang framework for coordination reduction [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]. Let us start from the Lagrangian dynamics in [START_REF] Brogliato | Nonsmooth Mechanics[END_REF].

Assumption 4

There is an open set Θ ⊂ R n-p and a function Ω : Θ → R p such that g(Ω(q 2 ), q 2 ) = 0 for all q 2 ∈ Θ.

This holds if in the neighborhood of any q such that g( q) = 0 the constraints are linearly independent 1 (possibly after some suitable coordinates reordering is performed).

Assumption 5 Assumption 4 holds with

Θ = R n-p .
Let us now describe the coordinate transformation. Let us partition q = q 1 q 2 , with q 1 ∈ R p and q 2 ∈ R n-p . The idea is that on the manifold M one has q 1 = Ω(q 2 ), so that the new bilateral constraints are expressed as

f (q) = q 1 -Ω(q 2 ) = 0. (62)
Define the nonlinear transformation X :

R n → R n by x = X(q) = q 1 -Ω(q 2 ) q 2
, that is differentiable with a differentiable inverse transformation Q : R n → R n , given by: q

= Q(x) = x 1 +Ω(x 2 )
x 2

. The vector partition x =

x 1

x 2 , with x 1 ∈ R p and x 2 ∈ R n-p is used. The Jacobian matrix of Q(•) at x is given by

T (x) = ∂Q ∂x = I p ∂Ω ∂x 2 (x 2 ) 0 I n-p ,
and it is full rank. It is easily verified that

∂X ∂q (q) = I p -∂Ω ∂q 2 0 I n-p = ∇f T (q) 0I n-p = T -1 (x).
The n × n identity matrix is partitioned as

I n = [E T 1 E T 2 ], E 1 ∈ R p×n , E 2 ∈ R (n-p)×n
. Then the dynamics in ( 6) is rewritten in the new coordinates as:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ E 1 M(x 2 )E T 2 ẍ2 + E 1 F (x 2 , ẋ2 , t) = E 1 T T (x 2 )∇ q h(x 2 )λ u + E 1 T T (x 2 )∇ q g(x 2 )λ b , E 2 M(x 2 )E T 2 ẍ2 + E 2 F (x 2 , ẋ2 , t) = E 2 T T (x 2 )∇ q h(x 2 )λ u , x 1 = 0, 0 ≤ λ u ⊥ h(x 2 ) ≥ 0.
(63) 1 One also says functionally independent.

Example 2 To illustrate the method let us recall the example given in [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]Sect. 7] of a planar Cartesian manipulator with an elliptic bilateral constraint, and we suppose in addition that a unilateral constraint acts on the manipulator. One has q T = (q 1 q 2 ), g(q) = 4(q 1 ) 2 + (q 2 ) 2 -1 and ∇g(q)λ b = (8q

1 λ b 2q 2 λ b ) T . Furthermore, x = X(q) = q 1 -0.5[1-(q 2 ) 2 ] 1 2 q 2
, and

T (x) = 1 -0.5[1-(q 2 ) 2 ] 1 2 0 1
, E 1 = (1 0), E 2 = (0 1). The transformed dynamics is given by

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -0.5x 2 (1 -x 2 2 ) -1 2 ẍ2 -0.5(1 -x 2 2 ) -3 2 ẋ2 2 = Λ u,1 + Λ b,1 , [1 + 0.25(1 -x 2 2 ) -1 x 2 2 ] ẍ2 + 0.25(1 -x 2 2 ) -2 x 2 ẋ2 2 = Λ u,2 , x 1 = 0, 0 ≤ λ u ⊥ h(x 2 ) ≥ 0, (64)
where Λ u = (Λ u,1 Λ u,2 ) T and Λ b = (Λ b,1 0) T are the generalized forces associated to the unilateral and the bilateral constraints, respectively, so that

Λ u,1 = E 1 T T (x 2 )∇ q h(x 2 )λ u , Λ u,2 = E 2 T T (x 2 )∇ q h(x 2 )λ u , Λ b,1 = E 1 T T (x 2 )∇ q g(x 2 )λ b .
The second equation in (63) represents the dynamics of the system in the manifold M (see [START_REF] Cottle | On a problem in linear inequalities[END_REF]), where λ b does not appear since all the constraints are perfect (no tangential forces since there is no friction). The first equation is an algebraic relation between the acceleration and the two multipliers. One has

M(x) = T T (x)M(Q(x))T (x), F (x, ẋ, t) = T T (x){F (Q(x), T (x) ẋ) + M(Q(x)) Ṫ (x) ẋ}. The notation f (x 2
) is for f (x) evaluated at x 1 = 0. Notice that using the above assumptions it follows that the reduced (np) × (np) mass matrix E 2 M(x 2 )E T 2 is positive definite. It is noteworthy that (63) is equivalent to (6), using that x 1 (t) = 0 for all t ≥ 0.

Notice that h(x 2 ) in (63) stands for the function h

(x) = h • Q • E 2 (x)
, and its gradient is given by ∇ h(x) = E 2 T T (x)∇ q h(q), where ∇ q h(q) denotes the transpose of the Jacobian ∂h ∂q at q with q = Q(x). According to our previous notation, we denote it as ∇ q h(x 2 ) in (63). Notice that E 2 T T (x)∇ q g(q) = ∇Ω(x 2 )I n-p ∇g(q) = ∇Ω(x 2 )∇ q 1 g(q) + ∇ q 2 g(q) = 0 (65) from Assumption 4. Also, E 1 T T (x 2 )∇ q g(x 2 ) = ∇ q 1 g(q), and ∇ q 1 g(q) is a p × p square matrix which is full rank if the constraints are linearly independent (possibly after reordering of the coordinates). Let us denote the multiplier associated with the constraints f (q) = 0 as λ f,b . Then one has

T T (x)∇ q f (q)λ f,b = T T (x) I p -∇Ω(q 2 ) λ f,b = I p 0 λ f,b = E 1 T T (x)∇ q g(q) E 2 T T (x)∇ q g(q) λ b = T T (x)∇ q g(q)λ b ( 66
)
from which we deduce that

λ f,b = ∇ q 1 g(q)λ b , ( 67
)
where λ b is associated with g(q) = 0 as in [START_REF] Brogliato | Nonsmooth Mechanics[END_REF]. It follows that in order to obtain the Lagrangian dynamics (63) one can also start from (62), and then apply the coordinate change q = Q(x).

Proceeding as above, we can construct the LCP:

0 ≤ λ u ⊥ w = A r (x 2 )λ u + H 5 (x 2 , ẋ2 , t) ≥ 0 (68) with A r (x 2 ) = ∇h T (q)T (x)E T 2 E 2 M(x 2 )E T 2 -1 E 2 T T (x)∇h(q) ≥ 0, (69)
and

H 5 (x 2 , ẋ2 , t) = d dt ∇h T (q)T (x)E T 2 ẋ2 -∇h T (q)T (x)E T 2 E 2 M(x 2 )E T 2 -1 E 2 F (x 2 , ẋ2 , t).
Since it is assumed that M(x 2 ) is full rank, the rank properties of the reduced order Delassus' matrix A r (x 2 ) depend only on the rank of the (np) × m matrix E 2 T T (x)∇h(q). Then exactly the same developments as in Sects. 2 and 4 can be redone to analyze the LCP in (68).

One sees that a necessary condition for A r (x 2 ) to be full rank is that m ≤ n-p ⇒ m+p ≤ n.

Using the first equation in (63), one can calculate λ b :

λ b = ∇ q 1 g(x 2 ) -1 E 1 M(x 2 )E T 2 E 2 M(x 2 )E T 2 -1 E 2 T T (x 2 )∇ q h(x 2 ) -∇ q 1 h(x 2 ) λ u + ∇ q 1 g(x 2 ) -1 E 1 F (x 2 , x2 , t) -E 1 M(x 2 )E T 2 E 2 M(x 2 )E T 2 -1 E 2 F (x 2 , x2 , t) , (70)
where we used that E 1 T T (x 2 )∇ q h(q) = ∇ q 1 h(x 2 ). Equation (70) corresponds to [START_REF] Frémond | Non-smooth Thermomechanics[END_REF], and (68) corresponds to [START_REF] Cottle | The Linear Complementarity Problem[END_REF]. Let ∇h p (q) ∈ R p×m denote the first p lines of ∇h(q), and ∇h n-p (q) ∈ R (n-p)×m its last np lines (i.e., ∇h p (q) = ∇ q 1 h(q) and ∇h n-p (q) = ∇ q 2 h(q)).

Proposition 15 A r (x 2 ) is positive definite if and only if ∇h n-p (x 2 ) + ∇Ω(x 2 )∇h p (x 2 ) has full rank np.

Proof One has

E 2 T T (x)∇h(q) = [0 I n-p ] I p 0 ∇Ω(x 2 ) I n-p ∇h(q) = ∇Ω(x 2 )∇ q h p (q) + ∇ q h n-p (q). One has A r (q) = (∇Ω(x 2 )∇ q h p (x 2 ) + ∇ q h n-p (x 2 )) T (E 2 M(x 2 )E T 2 ) -1 (∇Ω(x 2 )∇ q h p (x 2 ) + ∇ q h n-p (x 2 )). One infers that due to the positive definiteness of (E 2 M(x 2 )E T 2 ) -1 , A r (q) has rank n -p if and only if N (∇h n-p (x 2 ) + ∇Ω(x 2 )∇h p (x 2 )) = {0}.
Corollary 7 A r (x 2 ) is positive definite only if there exist no i ∈ {1, . . . , m} such that ∇h i (q) = α∇g j (q) for some j ∈ {1, . . . , p} and some α = 0.

Proof If the conditions of the corollary are not fulfilled, then from (65) there exists one index i ∈ {1, . . . , m} such that A r (q) has its ith diagonal element equal to zero.

Corollary 7 is coherent with Propositions 7, 8, and 14. The reduced coordinate method is a particular generalized coordinate transformation which allows one to write the gradient of the bilateral constraint in a specific form. It does not bring additional information on the constraints orthogonality, however. But the condition of Proposition 15 is easy to check. 6 Analysis at the impact times

The impact dynamics

All the above developments concern the dynamics outside the velocity discontinuities. At an impact time t , the Lagrangian dynamics becomes [START_REF] Brogliato | Nonsmooth Mechanics[END_REF]: M q(t) q t + -q t -= ∇h q(t) p u (t) + ∇g q(t) p b (t),

where p u (t) and p b (t) are the densities of the measures λ u and λ b , respectively, with respect to the Dirac measure with atom at t . Velocities are assumed to be of bounded variation so they have right and left limits at every t ≥ 0, and positions are continuous, hence q(t) is a constant at the impact. We denote q(t + ) = q+ = lim τ →t,τ >t q(τ ), and the same for the left limit. In particular, one has p u (t) ≥ 0, as we require that the contact force, considered as a measure, be positive. Since the bilateral constraints are satisfied at all times, one has ∇g T (q(t)) q(t + ) = ∇g T (q(t)) q(t -) = 0 for all t ≥ 0. Premultiplying Eq. ( 71) by ∇g T (q) and using Assumption 2, this allows one to compute that (the time argument is dropped):

p b = -∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 (q)∇h(q)p u . ( 72
)
Therefore, using (71) and ( 18), M(q) q+ -q-= ∇h(q) -∇g(q) ∇g T (q)M -1 (q)∇g(q) -1 ∇g T (q)M -1 (q)∇h(q) p u = I -G(q)M -1 (q) ∇h(q)p u .

(73)

Calculation of the percussion vector p u

We have the following proposition.

Proposition 16 If all the bilateral constraints are orthogonal in the kinetic metric to all the unilateral constraints, then at an impact time t one has p b (t) = 0.

Proof In this case, one has ∇g T (q)M -1 (q)∇h(q) = 0, and from (73) the result follows.

Example 3

In the example of Sect. 4.3, the constraints g(q) = βy + δ, β, δ ∈ R, undergo no percussion when the rod's tip hits the vertical wall; see Fig. 3. Let us provide some insight on this. The impact dynamics (71) for this rod is at an impact time t given by

⎧ ⎨ ⎩ m( ẋ(t + ) -ẋ(t -)) = p u (t), m( ẏ(t + ) -ẏ(t -)) = βp b (t), I ( θ(t + ) -θ(t -)) = l sin(θ (t))p u (t). ( 74
)
Since g(q) = 0, at all times it follows that y ≡ -δ β at all times and, therefore, from the second line in (74), one obtains that p b (t) = 0 as expected. Suppose now that the bilateral constraint is g(q) = αx + βy + γ for some constants α, β, γ (see the dashed prismatic joint on Fig. 3). One has ∇g T (q)M -1 (q)∇h(q) = α m , and the impact dynamics is given by ⎧

⎨ ⎩ m( ẋ(t + ) -ẋ(t -)) = p u (t) + αp b (t), m( ẏ(t + ) -ẏ(t -)) = βp b (t), I ( θ(t + ) -θ(t -)) = l sin(θ (t))p u (t).
(75) Remark 6 Notice from (72) that even if the system does not collide with all the unilateral constraints but only some of them, inertial couplings may induce a vector p b with no zero components. Proposition 16 shows that one may reduce the impacts effects inside the bilateral constraints, by designing mechanisms such that the bilateral constraints and the unilateral constraints have kinetic angles close to π 2 .

Let us denote qN = ∇h T (q) q, i.e. qN is the m-vector with entries ∇h T i (q) q. We may now adopt a collision rule as follows:

q+ N = -E q- N with: q- N ≤ 0, q+ N ≥ 0, h(q) = 0, ( 76
)
where E is the m × m restitution matrix (see, e.g., [6, pp. 299-304]). This is a generalized kinematic (Newton's like) impact rule, and its relationship with Moreau's impact law [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] is clarified later. Here, we assume that all the surfaces h i (q) = 0 are collided at the same time t , i.e., the m constraints are active at t . Inserting (76) into (73), one obtains

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ A c (q)p u = -(I + E) q- N , p u ≥ 0, E(I + E) -1 A c (q)p u ≥ 0. (77)
The constrained Delassus' matrix is the matrix that rules the calculation of p u , in addition to λ u in [START_REF] Lancaster | The Theory of Matrices[END_REF]. We may use some results of Sect. 4.2 to analyze (77). The last inequality comes from the kinematic constraints q+ N ≥ 0. The second one is a kinetic constraint stating the nonnegativity of the percussion vector. One has to add energetic constraints on the impact law to solve the impact problem, which is to be solved under three sets of constraints: energetic, kinetic, and kinematic [6, p. 302], [START_REF] Glocker | An introduction to impacts[END_REF]Eq. (5.9)]. Suppose that E = diag(e), e ≥ 0, then the equality in (77) becomes A c (q)p u = -(1 + e) q-N ≥ 0, the second inequality is e 1+e A c (q)p u = -e q-N ≥ 0 and is satisfied. If m = 1, the first inequality is satisfied as well since A c (q) is a positive scalar (in this case the two inequalities are identical).

Proposition 17

The equation in (77) has a solution p * u :

Proof In the case that p u,1 and p u,2 are solutions of (77), they satisfy p u,i = -A † c (q)(I + E) q-N + (I -A † c (q)A c (q))x, x ∈ R m , and (I -A † c (q)A c (q)) is the orthogonal projector onto N (A c (q)). Therefore, A c (q)p u,1 = A c (q)p u,2 . The result follows from the proofs of Proposition 6 and Corollary 1.

Proposition 18 is the counterpart of Propositions 2(i), Proposition 6 and Corollary 1, for the percussion vector. It is clear that the results of Sect. 4.2.2 may be used to guarantee that A c (q) > 0 so that p u is unique.

Remark 7 (Hyperstatic systems) The case of hyperstatic systems has been alluded to in Example 1. As its name indicates, this refers to statics, i.e., systems at an equilibrium (q, q) = (q * , 0) for some constant q * . One characterizes a static equilibrium as ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∇h(q * )λ * u = F (q * , 0, t), h(q * ) = 0,

λ * u ≥ 0, ( 78
)
where it is assumed that the system is fixed by the means of stops that can only push and not pull, hence the nonnegativity condition on the multiplier. The hyperstaticity exists if the problem (78) does not possess a unique solution λ * u , but several solutions (possibly an infinity), while ∇h(q * )λ u is unique. Provided that F (q * , 0, t) ∈ R(∇h(q * )) the solutions of (78) are given by λ * u = (∇h(q * )) † F (q * , 0, t) + [I -(∇h(q * )) † ∇h(q * )]x, for some x ∈ R dim(N (∇h(q)) , and with λ * u ≥ 0. One sees that the impact problem in (77) with p u ≥ 0 is a statics problem, similar to (78), though involving different terms. Hyperstaticity implies that m ≥ n (there must be at least as many constraints as degrees of freedom to guarantee that the equilibrium exists and persists). The complementarity problem that is studied in this paper is on the contrary a dynamical problem, whose solvability conditions differ from those of (78), as Propositions 2 and 3 show. However, if q(t + ) = 0 then the contact LCP in (2) becomes equal to 0 ≤ λ u t + ⊥ ∇h T q(t) M -1 q(t) ∇h q(t) λ u t + -F q(t), 0, t + ≥ 0.

(79)

Then the solutions of (78) solve also the contact LCP (79) at time t . Depending on F (q(t), 0, t + ) the LCP (79) may nevertheless possess other solutions for which the first equality in (78) is not satisfied. Finally, from Proposition 2, we know that the uniqueness of ∇h(q(t))λ u holds so that the dynamical system can be further integrated safely despite the contact force may not be uniquely determined. From the dynamics point of view, hyperstaticity is therefore not an issue.

Remark 8 (Impact geometry distortion) If all the unilateral constraints are orthogonal to all the bilateral ones, then the right-hand-side of (73) belongs to the normal cone -N Φ (q). In the general case of G(q) = 0 with G(q) in ( 18), it belongs to -[I -G(q)M -1 (q)]N φ (q) = {z ∈ R n |z = M(q)x -G(q)x, x ∈ -Nφ (q)}, where NΦ (q) is the normal cone to Φ at q in the kinetic metric. The bilateral constraints therefore introduce a distortion in the geometry of the problem, that may be measured by M(q) -G(q) since the right-hand side belongs to the set [G(q) -M(q)] Nφ (q). Using the same arguments as in the proof of Lemma 3, it follows that M(q) -G(q) ≥ 0, but it is not positive definite. The extension of the geometrical analysis of multiple impacts done in [START_REF] Glocker | Concepts for modeling impacts without friction[END_REF][START_REF] Glocker | An introduction to impacts[END_REF][START_REF] Glocker | The geometry of Newton's cradle[END_REF] is not tackled in this paper.

constraints. Then the post-impact velocity is the unique solution of:

(a)

⎧ ⎨ ⎩ min 1 2 ( q+ -q-) T M(q)( q+ -q-) subject to: ∇h T (q) q+ = q+ N ≥ 0

⇔ (b)

⎧ ⎨ ⎩ M(q)( q+ -q-) = ∇h(q)p u 0 ≤ p u ⊥ q+ N ≥ 0 (84) which is Moreau's impact rule.

Proof Due to the hypotheses, one has A c (q) = A(q) (see Proposition 5). Since we are dealing with an impact with p u > 0 and q-N < 0, and in view of the hypotheses, the program in (80) is equivalent to the program: ⎧ ⎨ ⎩ min 1 2 p T u A(q)p u + p T u q-N subject to: p u ≥ 0 (85) since they are both equivalent to their Karush-Kuhn-Tucker problem, which both possess the same unique solution. Its dual is given from Dorn's theorem by ⎧ ⎨ ⎩ min 1 2 p T u A(q)p u subject to: A(q)p u + q-N ≥ 0.

(86)

Using the dynamics (73), which reduces to M(q)( q+ -q-) = ∇h(q)p u , one obtains (84)(a). Notice that q+ N ≥ 0 is equivalent to q+ ∈ T Φ (q), i.e., the right velocity is admissible and points inward the admissible domain Φ. Therefore, the quadratic program in (84)(a) is exactly Moreau's impact rule; see [6, Eq. (5.66)] or [START_REF] Glocker | An introduction to impacts[END_REF]Eq. (5.20)]. The equivalence with (84)(b) follows from the Karush-Kuhn-Tucker's conditions applied to (84)(a) which are necessary and sufficient due to the problem's convexity.

An impact as in Proposition 19 is called an m-impact. It is clear that the result is not altered if less surfaces are attained simultaneously, where one simply has to adapt the matrices dimensions. The quadratic program in (84) may be seen as the impact counterpart of (48).

Conclusions

The problem considered in this paper is: given a Lagrangian system subject to a set of unilateral constraints (i.e., a complementarity Lagrangian system), how does the addition of a set of bilateral constraints alter the calculation of the Lagrange multipliers associated to the unilateral constraints? In other words, we study the properties of the extended Delassus' matrix when bilateral constraints also act on the system. Various criteria are proposed, and the kinetic angles between the constraints are extensively used. The link with Gauss' principle of Mechanics is made. An extended quadratic program for the acceleration is proposed. This work therefore extends the previous results by Moreau [START_REF] Moreau | Les liaisons unilatérales et le principe de Gauss[END_REF][START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF] and Lötstedt [START_REF] Lötstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. The case of impacts is also examined carefully. The major conclusion to be drawn from this study is that mixing unilateral and bilateral constraints does not destroy the problem's convexity, in the sense that the Delassus' matrix always remains at least positive semidefinite. This is therefore quite different from Coulomb's friction constraints that do destroy the convexity and yield other types of (much) more complex dynamics; see, e.g., [START_REF] Génot | New results on Painlevé paradoxes[END_REF][START_REF] Or | Investigation of Painlevé's paradox and dynamic jamming during mechanism sliding motion[END_REF].

Fig. 2

 2 Fig.2The constrained rod

Fig. 3

 3 Fig.3Orthogonal bilateral and unilateral constraints (at an impact time)

  Def. 3.1.1 xxiii)]; see Appendix A.7), one has from [20, Theorem 1, p. 194] (see Appendix A.6):

Acknowledgements This work was performed with the support of the NSFC/ANR project Multiple Impact, ANR-08-BLAN-0321-01.

(i) if I + E is full rank, and (I + E) -1 A c (q) and E(I + E) -1 A c (q) are monotone matrices; (ii) only if (I + E) q-N ∈ Q * Ac(q) ; (iii) only if (I + E) q-N ∈ R(∇h T (q)).

Proof (i) Let A c (q) be full rank. Then p u = -A c (q) -1 (I + E) q-N . Since -q-N ≥ 0 and since a monotone matrix is such that its inverse is nonnegative (it has only nonnegative entries) [20, p. 531], the first inequality in (77) is satisfied, and for the same reason the second inequality is satisfied, also.

(ii) (77) possesses a solution only if the LCP 0 ≤ p u ⊥ A c (q)p u + (I + E) q-N ≥ 0 is feasible (i.e., there exists p * u such that p * u ≥ 0 and A c (q)p * u + (I + E) q-N ≥ 0). Since A c (q) is positive semidefinite, hence copositive plus, this LCP is feasible if and only if [START_REF] Cottle | The Linear Complementarity Problem[END_REF]Corollary 3.8.10].

(iii) Solving (77) boils down to solving a linear equation of the form Ax = b, with additional inequality constraints. It is clear that a necessary condition for the existence of a solution of (77) is that

x, x ∈ R(∇g(q)) * }; see the proof of Lemma 5. Since R(∇g(q)) * = N (∇g T (q)), the necessary condition in Proposition 17(ii) can be checked. In any case when existence holds, the solutions of (77) are given by p u = -A † c (q)(I + E) q-N + (I -A † c (q)A c (q))x, where A † c (q) is the Moore-Penrose inverse of A c (q) and x ∈ R m , and (I -A † c (q)A c (q)) is the orthogonal projector onto N (A c (q)) [START_REF] Bernstein | Matrix, Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF]Prop. 6.1.6(xii)]. If m = 1, this gives a p u that satisfies the two inequalities. Let us investigate further the condition (ii) of Proposition 17.

Lemma 6

Let the bilateral constraints be orthogonal in the kinetic metric to the unilateral constraints. Let R((I + E)∇h T (q)) ⊆ R(∇ T (q)).

Proof In this case, one has A c (q) = A(q) = ∇h T (q)M -1 (q)∇h(q).

Under the lemma's condition, the result follows.

One notices that the lemma's second condition holds when E = diag(e). Also, we noticed that in this case the second inequality in (77) trivially holds, and (I + E) q+ N < 0 whenever q+ N < 0 which holds if an impact exists (the degenerate case q-N,i = 0 for some component i can be ignored since it yields q+ N,i = 0 for the diagonal restitution matrix case). The LCP 0

We infer that when q+ N < 0 and E = diag(e), e ≥ 0, these two LCPs are equivalent to (77) and the solution satisfies p * u > 0. Under these conditions, we recover Moreau's impact rule, according to the formulation of this impact rule in [START_REF] Glocker | An introduction to impacts[END_REF]Proposition 5.6]. It is worth recalling that the fact that the restitution matrix is equal to diag(e) = eI is crucial in these developments, for otherwise the analysis of (77) complicates. The conclusion of this discussion is that the impact law in (76) generalizes Moreau's law, at the price of being less tractable.

From (77), it is obvious that A c (q)p u is unique, that is if p u,1 and p u,2 are two solutions of (77) then A c (q)(p u,1p u,2 ) = 0. We can state the following. Proposition 18 Let p u,1 and p u,2 be two solutions of (77), equivalently of (80). Then ∇h(q)(p u,1p u,2 ) ∈ N (M -1 c (q)). If R(∇h(q)) ∩ R(∇g(q)) = {0}, then p u,1p u,2 ∈ N (∇h(q)).

Optimization framework

Let us assume that the impact law is well-posed, so that (77) has a solution p * u . Mimicking what has been done in Sect. 4.4, we may formulate the impact problem as follows. The percussion vector p * u is also a solution of the quadratic program:

subject to: p u ≥ 0,

Since A c (q), and hence A 2 c (q) are positive semidefinite, the quadratic program in (80) is equivalent to the LCP [8, p. 4]:

in the sense that if there exists a vector y and p * u that solve the LCP (81) then p * u solves (80), and vice versa. One sees that M is positive semi definite. One sees that the program in (80) differs from the one in (53) because there are additional inequality constraints. It is therefore not possible to deduce a program for q+ just by replacing q by q+ -qand letting terms like H i (q, q, t) vanish, in (56), without additional assumptions (see Proposition 19 below). From Dorn's duality and (80), one obtains:

and using the dynamics q+

subject to:

We recover here that despite the program (82) may not be strictly convex because A c (q) ≥ 0, the program (83) is strictly convex. It is also feasible since w = 0 and q+ N + E q-N = 0 satisfies the inequality constraints. The next result makes the link with Moreau's impact law.

Proposition 19

Suppose that A c (q) is positive definite, and let E = 0 (plastic impact). Assume that all the surfaces h i (q) = 0, 1 ≤ i ≤ m, are collided at the same time with q-N < 0, and that all the bilateral constraints are orthogonal in the kinetic metric to all the unilateral

Appendix: Useful mathematical results

A.1 Theorem 3.1.7 in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] (excerpts) Theorem 2 Let M ∈ R n×n be positive semidefinite, and let q ∈ R n be arbitrary. The following hold:

(a) If z 1 and z 2 are two solutions of the LCP(M, q) then (z 1 ) T (q + Mz 2 ) = (z 2 ) T (q + Mz 1 ). (d) If M is symmetric (as well as positive semidefinite) then Mz 1 = Mz 2 for any two solutions z 1 and z 2 .

A.2 Theorem 3.8.6 in [START_REF] Cottle | The Linear Complementarity Problem[END_REF] Theorem 3 Let M ∈ R n×n be copositive and let q ∈ R n be given. If the implication

Notice that the implication can also be written as q ∈ Q * M , with the notation of Remark 3.

A.3 Square root calculation

Consider the matrix A = 1 1 whose off-diagonal terms mimic the inertial couplings between the constraints. Since we can rewrite A = 1 0 0 1 + 0 0 , calculating the square root of A boils down to calculating the square root of the sum of two matrices: a diagonal matrix that mimics the matrices D h and D g (which do no depend on the inertial coupling parameter ), and a matrix with the off-diagonal terms and zero diagonal. Then lengthy but straightforward calculations yield

and:

where O( n ) denotes terms of degree higher or equal to n, and the approximations are valid if | | < 1. This simple example demonstrates that the approximations in Assumption 3 are not unrealistic.

A.4 Dorn's duality and converse duality theorems A.9 Exercise 8 in [START_REF] Lancaster | The Theory of Matrices[END_REF] p. 218

This exercise is reformulated here as a lemma for convenience.

Lemma 7

Let A be symmetric. Then the matrix I + A is positive definite for sufficiently small real numbers .