Keywords: Systèmes Lagrangiens, Probleme de complémentarité, Impacts, Stabilité, Poursuite des trajectoires, Commande basé sur la passivité, Systèmes non-réguliers

Commande passive pour la poursuite de trajectoires dans les systèmes Lagrangiens non-réguliers multi-contraintes

Résumé : Dans cette étude, on considère le problème de la commande pour assurer la poursuite des trajectoires pour une classe de systèmes Lagrangiens non-réguliers soumis à des contraintes unilatérales sans frottement. Une commande basée sur la passivité garantit certaines propriétés de stabilité pour le système en boucle fermée. Une attention particulière est accordée aux phases de transition avec impacts et aux phases de détachement. Ce travail étend précédents travaux sur le sujet car on considère des systèmes avec plusieurs contraintes et ndegré de liberté.

Introduction

The tracking control problem under consideration was studied in [START_REF] Brogliato | On the control of finite-dimensional mechanical systems with unilateral constraints[END_REF] mainly in the 1dof (degree-of-freedom) case and in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF] in the n-dof case. Both of these paper consider systems with only one unilateral frictionless constraint. Here we not only consider the multiconstraint case but the results in Section 7 relax some very hard to verify condition imposed in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF]. We note that in the case of a single nonsmooth impact the exponential stability and bounded-input bounded state (BIBS) stability was studied in [START_REF] Menini | A: Exponential and BIBS stabilisation of one degree of freedom mechanical system subject to single non-smooth impacts[END_REF] using a state feedback control law. A study for a multiple degrees-of-freedom linear systems subject to nonsmooth impacts can be found in [START_REF] Menini | A: Dynamic position feedback stabilisation of multidegrees-of-freedom linear mechanical system subject to nonsmooth impacts[END_REF]. That approach proposes a proportional-derivative control law in order to study BIBS stability via Lyapunov techniques. Another approach for the tracking control of nonsmooth mechanical systems (so-called billiards) can be found in [START_REF] Galeani | Trajectory tracking for a particle in elliptical billiards[END_REF][START_REF] Menini | Asymptotic tracking of periodic trajectories for a simple mechanical system subject to nonsmooth impacts[END_REF].

This paper focuses on the problem of tracking control of complementarity Lagrangian systems [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] subject to frictionless unilateral constraints whose dynamics may be expressed as:

   M (X) Ẍ + C(X, Ẋ) Ẋ + G(X) = U + ∇F (X)λ X 0 ≤ λ X ⊥ F (X) ≥ 0, Collision rule (1)
where X(t) ∈ R n is the vector of generalized coordinates, M (X) = M T (X) ∈ R n×n is the positive definite inertia matrix, F (X) ∈ R m represents the distance to the constraints, C(X, Ẋ) is the matrix containing Coriolis and centripetal forces, G(X) contains conservative forces, λ X ∈ R m is the vector of the Lagrangian multipliers associated to the constraints and U ∈ R n is the vector of generalized torque inputs. For the sake of completeness we precise that ∇ denotes the Euclidean gradient ∇F (X) = (∇F 1 (X), . . . , ∇F m (X)) ∈ R n×m where ∇F i (X) ∈ R n represents the vector of partial derivatives of F i (•) w.r.t. the components of X. We assume that the functions F i (•) are continuously differentiable and that ∇F i (X) = 0 for all X with F i (X) = 0. It is worth to precise here that for a given function f (•) its derivative w.r.t. the time t will be denoted by ḟ (•). For any function f (•) the limit to the right at the instant t will be denoted by f (t + ) and the limit to the left will be denoted by f (t -). A simple jump of the function f (•) at the moment t = t ℓ is denoted σ f (t ℓ ) = f (t + ℓ )f (t - ℓ ).

Definition 1 A Linear Complementarity Problem (LCP) is a system given by:

   λ ≥ 0 Aλ + b ≥ 0 λ T (Aλ + b) = 0 (2)
which is compactly re-written as

0 ≤ λ ⊥ Aλ + b ≥ 0 (3)
Such an LCP has a unique solution for all b if and only if A is a P-matrix [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF].

The admissible domain associated to the system (1) is the closed set Φ where the system can evolve and it is described as follows:

Φ = {X | F (X) ≥ 0} = 1≤i≤m Φ i ,
where Φ i = {X | F i (X) ≥ 0} considering that a vector is non-negative if and only if all its components are non-negative. In order to have a well-posed problem with a physical meaning we consider that Φ contains at least a closed ball of positive radius.

Definition 2 A singularity of the boundary ∂Φ of Φ is the intersection of two or more codimension one surfaces

Σ i = {X | F i (X) = 0}.
The presence of ∂Φ may induce some impacts that must be included in the dynamics of the system. It is obvious that m > 1 allows both simple impacts (when one constraint is involved) and multiple impacts (when singularities or surfaces of codimension larger than 1 are involved). In order to simplify the presentation we introduce the following notion of p ǫ -impact.

Definition 3 Let ǫ ≥ 0 be a fixed real number. We say that a p ǫ -impact occurs at the instant t if ||F I (X(t))|| ≤ ǫ, If ǫ = 0 the p surfaces Σ i , i ∈ I are stroked simultaneously. When ǫ > 0 the system collides ∂Φ in a neighborhood of the intersection i∈I Σ i .

INRIA

Definition 4 [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF] The tangent cone to Φ = {X | F i (X) ≥ 0, ∀i = 1, . . . , n} at q ∈ R n is defined as:

T Φ (q) = {z ∈ R n | z T ∇F i (q) ≥ 0, ∀i = J(q)}
where J(q) {i ∈ {1, . . . , n} | F i (q) ≤ 0}. When q ∈ Φ \ ∂Φ one has J(q) = ∅ and T Φ (q) = R n . The normal cone to Φ at q is defined as the polar cone to T Φ (•):

N Φ (q) = {y ∈ R n | ∀z ∈ T Φ (q), y T z ≤ 0}
The collision (or restitution) rule in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF], is a relation between the post-impact velocity and the pre-impact velocity. Among the various models of collision rules, Moreau's rule is an extension of Newton's law which is energetically consistent [START_REF] Glocker | Concepts for Modeling Impacts without Friction[END_REF] and is numerically tractable [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]. For these reasons throughout this paper the collision rule will be defined by Moreau's relation [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]:

Ẋ(t + ℓ ) = (1 + e) arg min z∈TΦ(X(t ℓ )) 1 2 [z -Ẋ(t - ℓ )] T × M (X(t ℓ ))[z -Ẋ(t - ℓ )] -e Ẋ(t - ℓ ) (4)
where Ẋ(t + ℓ ) is the post-impact velocity, Ẋ(t - ℓ ) is the pre-impact velocity and e ∈ [0, 1] is the restitution coefficient. Denoting by T the kinetic energy of the system, we can compute the kinetic energy loss at the impact t ℓ as [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF]:

T L (t ℓ ) = - 1 -e 2(1 + e) Ẋ(t + ℓ ) -Ẋ(t - ℓ )] T M (X(t ℓ )) × [ Ẋ(t + ℓ ) -Ẋ(t - ℓ ) ≤ 0 (5)
The collision rule can be rewritten considering the vector of generalized velocities as an element of the tangent space to the configuration space of the system, equipped with the kinetic energy metric. Doing so (see [START_REF] Brogliato | Nonsmooth Mechanics[END_REF] §6.2), the discontinuous velocity components Ẋnorm and the continuous ones Ẋtang are identified. Precisely,

Ẋnorm Ẋtang = M Ẋ, M = n T t T M (X)
where n ∈ R m represents the m unitary normal vec-

tors n i = M -1 (X)∇Fi(X) √ ∇Fi(X) T M -1 (X)∇Fi(X)
, i = 1, . . . , m and t represents nm mutually independent unitary vectors t i such that t T i M (X)n j = 0, ∀i, j. In this case the collision rule (4) at the impact time t ℓ becomes the generalized Newton's rule

Ẋnorm (t + ℓ ) Ẋtang (t + ℓ ) = -η Ẋnorm (t - ℓ ) Ẋtang (t - ℓ )
, η = diag(e 1 , . . . , e m , 0, . . . , 0) where e i is the restitution coefficient w.r.t. the surface Σ i . For the sake of simplicity we consider in this paper that all the restitution coefficients are equal, i.e. e 1 = . . . = e m e.

Remark 1 1) If

X ∈ Σ 1 Σ 2 and the angle ∠(Σ 1 , Σ 2 ) ≤ π then in the neighbor- hood of X one has Φ ≃ T Φ (X).
2) The case e = 0 is called a plastic impact and the case e = 1 is called an elastic impact. In the first case the normal component of the velocity becomes zero and in the second case the normal component of the velocity changes only its direction and preserves its magnitude. As we can easily see from [START_REF] Brogliato | Nonsmooth Mechanics[END_REF] in the second case there is no loss of kinetic energy at the impact moment.

3) One recalls that we deal with frictionless unilateral constraints. Some frictional contact laws that fit within the nonsmooth mechanic framework (1) can be found in [START_REF] Leine | Stability properties of equilibrium sets of nonlinear mechanical systems with dry friction and impact[END_REF].

The structure of the paper is as follows: in Section 2 one presents some basic concepts and prerequisites necessary for the further developments. Section 3 is devoted to the controller design. In Section 4 one defines the desired (or "exogenous") trajectories entering the dynamics. The desired contact-force that must occur on the phases where the motion is constrained, is explicitly defined in Section 5. Section 6 focuses on the strategy for take-off at the end of the constraint phases. The main results related to the closed-loop stability analysis are presented in Section 7. One example and concluding remarks end the paper.

The following standard notations will be adopted: || • || is the Euclidean norm, b p ∈ R p and b n-p ∈ R n-p are the vectors formed with the first p and the last np components of b ∈ R n , respectively. N Φ (X p = 0) is the normal cone N Φ (X) to Φ at X [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF] when X satisfies X p = 0, λ min (•) and λ max (•) represent the smallest and the largest eigenvalues, respectively.

Basic concepts

Typical task

In the case m = 1 (only one unilateral constraint) the time axis can be split into intervals Ω k and I k corresponding to specific phases of motion [START_REF] Brogliato | On the control of finite-dimensional mechanical systems with unilateral constraints[END_REF]. Precisely, Ω 2k corresponds to free-motion phases (F (X) > 0) and Ω 2k+1 corresponds to constrainedmotion phases (F (X) = 0). Between free and constrained phases the dynamical system always passes into a transition phase I k containing some impacts. Since the dynamics of the system does not change during the transition between constrained and free-motion phases, in the time domain one gets the following typical task representation:

R + = Ω 0 ∪ I 0 ∪ Ω 1 ∪ Ω 2 ∪ I 1 ∪ . . . ∪ Ω 2k ∪ I k ∪ Ω 2k+1 ∪ . . . (6)
In the case m ≥ 2 (multiple constraints) things complicate since the number of typical phases increases due to the singularities that must be taken into account. Explicitly, the constrained-motion phases need to be decomposed in sub-phases where some specific constraints are active. Between two such sub-phases a transition phase occurs when the number of active constraints increases. Nevertheless, a typical task can be represented in the time domain as:

R + = k≥0 Ω J k 2k ∪ I J k k ∪ m k i=1 Ω J k,i 2k+1 J k ⊂ J k,1 ; J k+1 ⊂ J k,m k ⊂ J k,m k -1 ⊂ . . . J k,1 (7)
where the superscript J k represents the set of active constraints (J k = {i ∈ {1, . . . , m} | F i (X) = 0}) during the corresponding motion phase, and I J k k denotes the transient between two Ω k phases when the number of active constraints increases. When the number of active constraints decreases there is no impact, thus no other transition phases are needed. We note that J k = ∅ corresponds to free-motion.

INRIA

For the sake of simplicity and without any loss of generality we replace

m k i=1 Ω J k,i 2k+1 by Ω J ′ k 2k+1 where J k ⊂ J ′ k and J k+1 ⊂ J ′ k .
Therefore the typical task simplifies as:

R + = k≥0 Ω J k 2k ∪ I J k k ∪ Ω J ′ k 2k+1 J k ⊂ J ′ k , J k+1 ⊂ J ′ k (8)
Since the tracking control problem involves no difficulty during the Ω k phases, the central issue is the study of the passages between them (the design of transition phases I k and detachment conditions), and the stability of the trajectories evolving along [START_REF] Dzonou | A sweeping process approach to inelastic contact problems with general inertia operators[END_REF] (i.e. an infinity of cycles). Throughout the paper, the sequence

Ω J k 2k ∪ I J k k ∪ Ω J ′ k 2k+1
will be referred to as the cycle k of the system's evolution. For robustness reasons during transition phases I k we impose a closed-loop dynamics (containing impacts) that mimics somehow the bouncing-ball dynamics (see e.g. [START_REF] Brogliato | Nonsmooth Mechanics[END_REF]).

Stability analysis criteria

The system (1) is a complex nonsmooth and nonlinear dynamical system which involves continuous and discrete time phases. A stability framework for this type of systems has been proposed in [START_REF] Brogliato | On the control of finite-dimensional mechanical systems with unilateral constraints[END_REF] and extended in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF]. This is an extension of the Lyapunov second method adapted to closed-loop mechanical systems with unilateral constraints. Since we use this criterion in the following tracking control strategy it is worth to clarify the framework and to introduce some definitions.

Let us introduce the trajectories playing a role in the dynamics and the design of the controller:

• X nc (•) denotes the desired trajectory of the unconstrained system (i.e. the trajectory that the system should track if there were no constraints). We suppose that F (X nc (t)) < 0 for some t, otherwise the problem reduces to the tracking control of a system with no constraints.

• X * d (•) denotes the signal entering the control input and playing the role of the desired trajectory during some parts of the motion.

• X d (•) represents the signal entering the Lyapunov function. This signal is set on the boundary ∂Φ after the first impact of each cycle.

These signals may coincide on some time intervals as we shall see later. In order to clarify the differences and the usefulness of these signals we present a simple example concerning a one-degree of freedom model.

Example 1 Let us consider the following one degree-of-freedom dynamical system:

   ( Ẍ -Ẍ * d ) + γ 2 ( Ẋ -Ẋ * d ) + γ 1 (X -X * d ) = λ 0 ≤ X ⊥ λ ≥ 0 Ẋ(t + k ) = -e n Ẋ(t - k ) (9)
where X * d (•) is a twice differentiable signal and γ 1 , γ 2 are two gains. The real number λ represents the Lagrange multiplier associated to the constraint X ≥ 0.

The system [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF] represents the dynamics of a unitary mass lying in the right half plane (X ≥ 0) and moving on the Ox axis (see figure 1). Thus, the system dynamics is 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000

= Ẍ * d -γ 2 ( Ẋ -Ẋ * d ) -γ 1 (X -X * d ).
One considers that the Lyapunov function V (X, Ẋ, t) of the closed-loop system is given by a quadratic function of the tracking error.

In order to stabilize the system on the constraint surface ∂Φ there are two solutions.

[1st solution] Setting X * d (t) = 0, Ẋ * d (t) = 0 and Ẍ * d (t) = 0 at the time t of contact with ∂Φ. Therefore when the system reaches its equilibrium point on ∂Φ one obtains the contact force λ = 0. It is noteworthy that, since we are interested to study a dynamics containing a constraint movement phase, we need a non-zero contact force when the system is stabilized on ∂Φ. Therefore this type of solution is not convenient for the study developed in this work. Moreover the stabilization on ∂Φ with a non zero tracking error on [tǫ, t] is not obvious, especially in higher dimensions, so that the conditions for a perfect tangential approach are never met in practice.

[2nd solution] Setting X * d (t) = -a, a > 0, Ẋ * d (t) = 0 and Ẍ * d (t) = 0 on some interval of time when approaching ∂Φ. Since the tracking error X = X -X * d = -X * d at the equilibrium point we cannot use X in the definition of the Lyapunov function.

Obviously the system reaches the desired position when X = 0. Thus we shall use the signal X d (t) = 0 in order to express the desired trajectory in the Lyapunov function and the corresponding tracking error will be given by X = X -X d . Concluding, when the equilibrium point is reached one obtains X = 0, X = -a, V (t, X = 0) = 0 and the corresponding contact force λ = -γ 1 X * d = γ 1 a > 0.

Next, let us define Ω as the complement of I = k≥0

I J k k and
Consider I J k k ∆ = [τ k 0 , t k f ]
and V (•) such that there exists strictly increasing functions α(•) and β(•) satisfying the conditions: α(0) = 0, β(0) = 0 and α(||x||) ≤ V (x, t) ≤

INRIA β(||x||).

In the sequel, we consider that for each cycle k the sequence of impact instants t k ℓ has an accumulation point t k ∞ . We note that a finite accumulation period (i.e. t k ∞ < +∞) implies that e < 1 (in [START_REF] Ballard | Formulation and well-posedness of the dynamics of rigid body systems with perfect unilateral constraints[END_REF] it is shown that e = 1 implies that t k ∞ = +∞). The following criterion is inspired from [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF], and will be used for studying the stability of system (1).

Proposition 1 (Weak Stability) Assume that the task admits the representation [START_REF] Dzonou | A sweeping process approach to inelastic contact problems with general inertia operators[END_REF] and that

a) λ[I J k k ] < +∞, ∀k ∈ N, b) outside the impact accumulation phases [t k 0 , t k ∞ ] one has V (x(t), t) ≤ -γV (x(t), t) for some constant γ > 0, c) ℓ≥0 V (t k- ℓ+1 ) -V (t k+ ℓ ) ≤ K 1 V p1 (τ k 0 ), ∀ℓ ≥ 0 for some p 1 ≥ 0, K 1 ≥ 0, d) the system is initialized on Ω 0 such that V (τ 0 0 ) ≤ 1, e) ℓ≥0 σ V (t k ℓ ) ≤ K 2 V p2 (τ k 0 ) + ξ for some p 2 ≥ 0, K 2 ≥ 0 and ξ ≥ 0. If p = min{p 1 , p 2 } < 1 then V (τ k 0 ) ≤ δ(γ, ξ)
, where δ(γ, ξ) is a function that can be made arbitrarily small by increasing the value of γ. The system is practically weakly stable with R = α -1 (δ(γ, ξ)).

Proof:

From assumption (b) one has

V (t k f ) ≤ V (t k ∞ )e -γ(t k f -t k ∞ )
It is clear that condition (c) combined with (e) leads to

V (t k ∞ ) ≤ V (τ k 0 ) + K 1 V p1 (τ k 0 ) + K 2 V p2 (τ k 0 ) + ξ Considering p < 1, the assumption (d) guarantees that max{V (τ k 0 ), V p1 (τ k 0 ), V p2 (τ k 0 )} ≤ V p (τ k 0 ) ≤ 1 and we get V (t k f ) ≤ e -γ(t k f -t k ∞ ) [1 + K 1 + K 2 + ξ] V p (τ k 0 ) ≤ e -γ(t k f -t k ∞ ) [1 + K 1 + K 2 + ξ] δ(γ, ξ) From assumption (b) one has V (τ k+1 0 ) ≤ V (t k f )
thus the first part of the statement holds. The term δ(γ, ξ) can be made as small as desired increasing either γ or the length of the interval

[t k ∞ , t k f ]. The proof is completed by the relation α(||x||) ≤ V (x, t), ∀x, t
Remark 2 Since the Lyapunov function is exponentially decreasing on the Ω k phases, assumption (d) in Proposition 1 means that the system is initialized on Ω 0 sufficiently far from the moment when the trajectory X nc (•) leaves the admissible domain.

Precisely, the weak stability is characterized by an "almost decreasing" Lyapunov function V (x(•).•) as illustrated in Figure 2.

Remark 3 It is worth to point out the local character of the stability criterion proposed by Proposition 1. This character is firstly given by condition d) of the statement and

secondly by the synchronization constraints of the control law and the motion phase of the system (see [START_REF] Galeani | Trajectory tracking for a particle in elliptical billiards[END_REF] and ( 8) below).

The practical stability is very useful because attaining asymptotic stability is not an easy task for the unilaterally constrained systems described by (1) especially when n ≥ 2 and M (q) is not a diagonal matrix (i.e. there are inertial couplings, which is the general case).

Controller design

In order to overcome some difficulties that can appear in the controller definition, the dynamical equations (1) will be expressed in the generalized coordinates introduced by McClamroch & Wang [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]. We suppose that the generalized coordinates transformation holds globally in Φ, which may obviously not be the case in general. However, the study of the singularities that might be generated by the coordinates transformation is out of the scope of this paper. Let us consider

D = [I m . . . O] ∈ R m×n , I m ∈ R m×m the identity matrix. The new coordinates will be q = Q(X) ∈ R n , with q = q 1 q 2 , q 1 =    q 1 1 . . . q m 1    such that Φ = {q | Dq ≥ 0} 1 . The tangent cone T Φ (q 1 = 0) = {v | Dv ≥ 0}
is the space of admissible velocities on the boundary of Φ.

The controller used here consists of different low-level control laws for each phase of the system. More precisely, the switching controller can be expressed as

T (q)U =    U nc for t ∈ Ω ∅ 2k U J t for t ∈ I J k U J c for t ∈ Ω J k (10)
where T (q) = T 1 (q) T 2 (q) ∈ R n×n is full-rank under some basic assumptions (see [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]). The dynamics becomes: [START_REF] Glocker | Set Value Force Laws: Dynamics of Non-Smooth Systems Lecture Notes in Applied Mechanics[END_REF] where the set of complementary relations can be written more compactly as 0 ≤ λ ⊥ Dq ≥ 0.

       M 11 (q)q 1 + M 12 (q)q 2 + C 1 (q, q) q + g 1 (q) = T 1 (q)U + λ M 21 (q)q 1 + M 22 (q)q 2 + C 2 (q, q) q + g 2 (q) = T 2 (q)U q i 1 ≥ 0, q i 1 λ i = 0, λ i ≥ 0, 1 ≤ i ≤ m Collision rule
In the sequel U nc coincides with the fixed-parameter controller proposed in [START_REF] Sadegh | Stability analysis of an adaptive controller for robotic manipulators[END_REF][START_REF] Slotine | Adaptive manipulator control: A case study[END_REF] and the closed-loop stability analysis of the system is based on Proposition 1. First, let us introduce some notations: q = qq d , q = qq * d , s = q + γ 2 q, s = q + γ 2 q, qe = qdγ 2 q where γ 2 > 0 is a scalar gain and q d , q * d represent the desired trajectories defined in the previous section. Using the above notations the controller is given by

T (q)U        U nc = M (q)q e + C(q, q) qe + G(q) -γ 1 s U J t = U J nc , t ≤ t k 0 U J t = M (q)q e + C(q, q) qe + G(q) -γ 1 s, , t > t k 0 U J c = U nc -P d + K f (P q -P d ) (12)
where γ 1 > 0 is a scalar gain, K f > 0, P q = D T λ and P d = D T λ d is the desired contact force during persistently constrained motion. It is clear that during Ω J k not all the constraints are active and, therefore, some components of λ and λ d are zero.

In order to prove the stability of the closed-loop system ( 11) -( 12) we will use the following positive definite function:

V (t, s, q) = 1 2 s T M (q)s + γ 1 γ 2 qT q (13)

Tracking control framework

In this paper we treat the tracking control problem for the closed-loop dynamical system ( 10)-( 12) with the complete desired path a priori taking into account the complementarity conditions and the impacts. In order to define the desired trajectory let us consider the motion of a virtual and unconstrained particle perfectly following a trajectory (represented by X nc (•) on Figure 3) with an orbit that leaves the admissible domain for a given period. Therefore, the orbit of the virtual particle can be split into two parts, one of them belonging to the admissible domain (inner part) and the other one outside the admissible domain (outer part). In the sequel we deal with the tracking control strategy when the desired trajectory is constructed such that:

(i) when no activated constraints, it coincides with the trajectory of the virtual particle (the desired path and velocity are defined by the path and velocity of the virtual particle, respectively), (ii) when p ≤ m constraints are active, its orbit coincides with the projection of the outer part of the virtual particle's orbit on the surface of codimension p defined by the activated constraints (X d between A ′′ and C in Figure 3), (iii) the desired detachment moment and the moment when the virtual particle reenters the admissible domain (with respect to p ≤ m constraints) are synchronized.

Therefore we have not only to track a desired path but also to impose a desired velocity allowing the motion synchronization on the admissible domain.

The main difficulties here consist of:

• stabilizing the system on ∂Φ during the transition phases I J k k and incorporating the velocity jumps in the overall stability analysis;

• deactivating some constraints at the moment when the unconstrained trajectory re-enters the admissible domain with respect to them;

• maintaining a constraint movement between the moment when the system was stabilized on ∂Φ and the detachment moment.

Remark 4

The problem can be relaxed considering that we want to track only a desired path like X nc (•) (without imposing a desired velocity on the inner part of the desired trajectory and/or a given period to complete a cycle). In this way the synchronization problem (iii) disappears and we can assume there exists a twice differentiable desired trajectory outside [t k 0 , t k f ] that assures the detachment when the force control is dropped. In other words, in this case we have to design the desired trajectory only during I J k k phases.

Design of the desired trajectories

Throughout the paper we consider

I J k k = [τ k 0 , t k f ],
where τ k 0 is chosen by the designer as the start of the transition phase

I J k k and t k f is the end of I J k k .
We note that all superscripts (•) k will refer to the cycle k of the system motion. We also use the following notations:

• t k 0 is the first impact during the cycle k,

• t k ∞ is the accumulation point of the sequence {t k ℓ } ℓ≥0 of the impact instants during the cycle

k (t k f ≥ t k ∞ ),
• τ k 1 will be explicitly defined later and represents the instant when the desired signal X * d reaches a given value chosen by the designer in order to impose a closed-loop dynamics with impacts during transition phases,

• t k d is the desired detachment instant, therefore the phases Ω

J ′ k 2k+1 can be expressed as [t k f , t k d ].
It is noteworthy that t k 0 , t k ∞ , t k d are state-dependent whereas τ k 1 and τ k 0 are exogenous and imposed by the designer. To better understand the definition of these specific instants, in the Figure 3 we simplify the system's motion as follows:

• during transition phases we take into account only the constraints that must be activated

J ′ k \ J k .
• at the end of Ω 2k+1 phases we take into account only the constraints that must be deactivated

J ′ k \ J k+1 .
The points A, A ′ , A ′′ and C in Figure 3 correspond to the moments τ k 0 , t k 0 , t k f and t k d respectively. We have seen that the choice of τ k 0 plays an important role in the stability criterion given by Proposition 1. On the other hand in Figure 3 we see that

INRIA A A'' B A' C Φ ∂Φ X * d (t) = X d (t) X d (t) X * d (t) X nc (t) = X * d (t) = X d (t) X nc (t)
X d (•) = X * d (•) is deformed compared to X nc (•).
In order to reduce this deformation τ k 0 and implicitly the point A must be close to ∂Φ (see also Figure 6). Further details on the choice of τ k 0 will be given later. Taking into account just the constraints J ′ k \ J k+1 we can identify t k d with the moment when X d (•) and X nc (•) rejoin at C.

Design of q *

d (•) and q d (•) on the phases

I J k k
During the transition phases the system must be stabilized on ∂Φ. Obviously, this does not mean that all the constraints have to be activated (i.e. q i 1 (t) = 0, ∀i = 1, . . . , m). Let us consider that only the first p constraints (eventually reordering the coordinates) define the border of Φ where the system must be stabilized. On [τ k 0 , t k 0 ) we define q * d (•) as a twice differentiable signal such that q * d (•) approaches a given point in the normal cone N Φ (q p = 0) on [τ k 0 , τ k 1 ]. Precisely, we define q * d (•) such as: • during a small period δ > 0 chosen by the designer the desired velocity becomes zero preserving the twice differentiability of q * d (•). For instance we can use the following definition:

q * d (t) = q nc τ k 0 + (t -τ k 0 -δ) 2 (t -τ k 0 ) δ 2 , t ∈ [τ k 0 , τ k 0 + δ] which means q * d (τ k 0 + δ) = q * d (τ k 0 ) = q nc (τ k 0 ), q * d (τ k 0 + δ) = 0 and q * d (τ k 0 ) = qnc (τ k 0 )
• choosing ϕ > 0 and denoting

t ′ = t-(τ k 0 +δ) τ k 1 -(τ k 0 +δ) , the components q i d * , i = 1, .
. . , p of (q * d ) p are defined as:

q i d * (t) = a 3 (t ′ ) 3 + a 2 (t ′ ) 2 + a 0 , t ∈ [τ k 0 + δ, min{τ k 1 ; t k 0 }] -ϕV 1/3 (τ k 0 ), t ∈ (min{τ k 1 ; t k 0 }, t k f ] (14)
where V (•) is defined in [START_REF] Leine | Stability properties of equilibrium sets of nonlinear mechanical systems with dry friction and impact[END_REF] and with the coefficients given by:

a 3 = 2[ q i nc (τ k 0 ) + ϕV 1/3 (τ k 0 )] a 2 = -3[ q i nc (τ k 0 ) + ϕV 1/3 (τ k 0 )] a 0 = q i nc (τ k 0 ) (15)
RR n°6483

• all the other components of q * d (•) are frozen:

(q * d ) n-p (t) = q nc n-p (τ k 0 ), t ∈ (τ k 0 + δ, t k f ] (16)
The rationale behind the choice of q * d (•) is on one hand to assure a robust stabilization on ∂Φ, mimicking the bouncing-ball dynamics; on the other hand to enable one to compute suitable upper-bounds that will help using Proposition 1 (hence V 1/3 (•) terms in [START_REF] Mabrouk | A unified variational model for the dynamics of perfect unilateral constraints[END_REF] with V (•) in ( 13)).

Remark 5 1) Straightforward computations show that q * d (•) satisfies the following relations.

q i d * (τ k 1 ) = -ϕV 1/3 (τ k 0 ), qi d * (τ k 1 ) = 0, i = 1, . . . , p 2
) Two different situations are possible. The first is given by t k 0 > τ k 1 (see Figure 4) and we shall prove that in this situation all the jumps of the Lyapunov function in [START_REF] Leine | Stability properties of equilibrium sets of nonlinear mechanical systems with dry friction and impact[END_REF] are negative. The second situation was pointed out in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF] and is given by t k 0 < τ k 1 . In this situation the first jump at t k 0 in the Lyapunov function is positive, therefore the system can be only weakly stable. It is noteworthy that q * d (•) will then have a jump at the time t k 0 since (q i d ) * (t k+ 0 ) = -ϕV 1/3 (τ k 0 ), ∀i = 1, . . . , p (see ( 14)).

In order to limit the deformation of the desired trajectory q * d (•) w.r.t. the unconstrained trajectory q nc (•) during the I k phases (see Figure 3 and4), we impose in the sequel

||q nc p (τ k 0 )|| ≤ ψ (17)
where ψ > 0 is chosen by the designer. It is obvious that a smaller ψ leads to smaller deformation of the desired trajectory and to smaller deformation of the real trajectory as we shall see in Section 8. Nevertheless, due to the tracking error, ψ cannot be chosen zero. We also note that ||q nc p (τ k 0 )|| ≤ ψ is a practical way to choose

τ k 0 . O B C A t t k 0 t k 1 t k d (q * d ) i (t) q i 1 q i 1 (t) τ k 0 t k 0 > τ k 1 Ω 2k I k Ω 2k+1 Ω 2k+2 A ′ t k f τ k 1 (q * d ) i (t) (q * d ) i (t) -ϕV 1/3 (τ k 0 )
Figure 4: The design of q * 1d on the transition phases I k

During the transition phases I k we define (q d ) n-p (t) = (q * d ) n-p (t). Assuming a finite accumulation period, the impact process can be considered in some way equivalent to a plastic impact. Therefore, (q d ) p (•) and ( qd ) p (•) are set to zero on the right of t k 0 .

INRIA 5 Design of the desired contact force during constraint phases

For the sake of simplicity we consider the case of the constraint phase Ω J k , J = ∅ with J = {1, . . . , p}. Obviously a sufficiently large desired contact force P d assures a constrained movement on Ω J k . Nevertheless at the end of the Ω J 2k+1 phases a detachment from some surfaces Σ i has to take place. It is clear that a take-off implies not only a well-defined desired trajectory but also some small values of the corresponding contact force components. On the other hand, if the components of the desired contact force decrease too much a detachment can take place before the end of the Ω J k phases which can generate other impacts. Therefore we need a lower bound of the desired force which assures the contact during the Ω J k phases. Dropping the time argument, the dynamics of the system on Ω J k can be written as

M (q)q + F (q, q) = U c + D T p λ p 0 ≤ q p ⊥ λ p ≥ 0 ( 18
)
where F (q, q) = C(q, q) q + G(q) and D p = [I p . . . O] ∈ R p×n . On Ω J k the system is permanently constrained which implies q p (•) = 0 and qp (•) = 0. In order to assure these conditions it is sufficient to have λ p > 0.

In the following let us denote

M -1 (q) = [M -1 (q)] p,p [M -1 (q)] p,n-p [M -1 (q)] n-p,p [M -1 (q)] n-p,n-p
and C(q, q) = C(q, q) p,p C(q, q) p,n-p C(q, q) n-p,p C(q, q) n-p,n-p where the meaning of each component is obvious.

Proposition 2

On Ω J k the constraint motion of the closed-loop system [START_REF] Menini | A: Dynamic position feedback stabilisation of multidegrees-of-freedom linear mechanical system subject to nonsmooth impacts[END_REF], [START_REF] Galeani | Trajectory tracking for a particle in elliptical billiards[END_REF],( 12) is assured if the desired contact force is defined by

(λ d ) p β - Mp,p (q) 1 + K f [M -1 (q)] p,p C p,n-p (q, q)+ [M -1 (q)] p,n-p C n-p,n-p (q, q) + γ 1 [M -1 (q)] p,n-p s n-p (19)
where Mp,p (q) = [M -1 (q)] p,p -1 = D p M -1 (q)D T p -1 is the inverse of the Delassus' matrix and β ∈ R p , β > 0.

Proof: First, we notice that the second relation in [START_REF] Menini | A: Dynamic position feedback stabilisation of multidegrees-of-freedom linear mechanical system subject to nonsmooth impacts[END_REF] implies on Ω J k (see [START_REF] Glocker | Set Value Force Laws: Dynamics of Non-Smooth Systems Lecture Notes in Applied Mechanics[END_REF])

0 ≤ qp ⊥ λ p ≥ 0 ⇔ 0 ≤ D p q ⊥ λ p ≥ 0. (20)
From ( 18) and ( 12) one easily gets:

q = M -1 (q) -F (q, q) + U nc + (1 + K f )D T p (λ -λ d ) p
Combining the last two equations we obtain the following LCP with unknown λ p :

0 ≤ D p M -1 (q) -F (q, q) + U nc -(1 + K f )D T p (λ d ) p +(1 + K f )D p M -1 (q)D T p λ p ⊥ λ p ≥ 0 (21)
Since (1 + K f )D p M -1 (q)D T p > 0 and hence is a P-matrix, the LCP (21) has a unique solution and one deduces that λ p > 0 if and only if

Mp,p (q) 1 + K f D p M -1 (q) U nc -F (q, q) -(1 + K f )D T p (λ d ) p < 0 ⇔ (λ d ) p > Mp,p (q) 1 + K f D p M -1 (q) U nc -F (q, q) ⇔ (λ d ) p = β + Mp,p (q) 1 + K f D p M -1 (q) U nc -F (q, q) (22)
with β ∈ R p , β > 0. Since U nc -F (q, q) = M (q)q e -C(q, q)sγ 1 s, (q e ) p = 0 and s p = 0, (22) rewrites as [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] and the proof is finished. It is noteworthy that

λ p = - Mp,p (q) 1 + K f D p M -1 (q) U nc -F (q, q) -(1 + K f )D T p (λ d ) p = (λ d ) p - Mp,p (q) 1 + K f D p M -1 (q) U nc -F (q, q) = β Remark 6
The control law used in this paper with the design of λ d described above leads to the following closed-loop dynamics on Ω 2k+1 .

   M p,n-p (q) ṡn-p + C p,n-p (q, q)s n-p = (1 + K f )(λ -λ d ) p M n-p,n-p (q) ṡn-p + C n-p,n-p (q, q)s n-p + γ 1 s n-p = 0 q p = 0, λ p = β
It is noteworthy that the closed-loop dynamics is nonlinear and therefore, we do not use the feedback stabilization proposed in [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF].

Strategy for take-off at the end of constraint phases Ω J 2k+1

We have discussed in the previous sections the necessity of a trajectory with impacts in order to assure the robust stabilization on ∂Φ in finite time and, the design of the desired trajectory to stabilize the system on ∂Φ. Now, we are interested in finding the conditions on the control signal U J c that assure the take-off at the end of constraint phases Ω J 2k+1 . As we have already seen before, the phase

Ω J 2k+1 corresponds to the time interval [t k f , t k d ). The dynamics on [t k f , t k d
) is given by (18) and the system is permanently constrained, which implies q p (•) = 0 and qp (•) = 0. Let us also consider that the first r constraints (r < p) have to be deactivated. Thus, the detachment takes place at t k d if qr (t k+ d ) > 0 which requires λ r (t k- d ) = 0. The last pr constraints remain active which means λ p-r (t k- d ) > 0. To simplify the notation we drop the time argument in many equations of this section. We decompose the LCP matrix (which is the Delassus' matrix multiplied by 1 + K f ) as:

(1 + K f )D p M -1 (q)D T p = A 1 (q) A 2 (q) A 2 (q) T A 3 (q) with A 1 ∈ R r×r , A 2 ∈ R r×(p-r) and A 3 ∈ R (p-r)×(p-r)

INRIA

Proposition 3 For the closed-loop system [START_REF] Menini | A: Dynamic position feedback stabilisation of multidegrees-of-freedom linear mechanical system subject to nonsmooth impacts[END_REF], [START_REF] Galeani | Trajectory tracking for a particle in elliptical billiards[END_REF], [START_REF] Glocker | Concepts for Modeling Impacts without Friction[END_REF] the passage when the number of active constraints decreases from p to r < p, is possible if

(λ d ) r (t k d ) (λ d ) p-r (t k d ) = A 1 -A 2 A -1 3 A T 2 -1 b r -A 2 A -1 3 b p-r -C 1 C 2 + A -1 3 b p-r -A T 2 (λ d ) r ( 23
)
where b p b(q, q, U nc ) D p M -1 (q)[U nc -F (q, q)] ≥ 0 and C 1 ∈ R r , C 2 ∈ R p-r such that C 1 ≥ 0, C 2 > 0.
Proof: From ( 12) and ( 18) one gets

qp (t) = b p + (1 + K f )D p M -1 (q)D T p (λ -λ d )
Therefore the LCP (20) rewrites as:

0 ≤ λ r λ p-r ⊥ b r + A 1 (λ -λ d ) r + A 2 (λ -λ d ) p-r b p-r + A T 2 (λ -λ d ) r + A 3 (λ -λ d ) p-r ≥ 0 (24)
Under the conditions λ r = 0 and λ p-r > 0 one has

0 ≤ λ p-r ⊥ b p-r -A T 2 (λ d ) r + A 3 (λ -λ d ) p-r ≥ 0
with the solution

λ p-r = -A -1 3 b p-r -A T 2 (λ d ) r -A 3 (λ d ) p-r (25)
Thus λ p-r > 0 is equivalent to

(λ d ) p-r > A -1 3 b p-r -A T 2 (λ d

Proposition 4

The closed-loop system ( 18),( 10),( 12) is permanently constrained on

[t k f , t k d ) and a smooth detachment is guaranteed on [t k d , t k d + ǫ) (ǫ is a small positive real number chosen by the designer) if (i) (λ d ) p (•) is defined on [t k f , t k d ) by (23) where C 1 is replaced by C 1 (t -t k d ). (ii) On [t k d , t k d + ǫ) q * d (t) = q d (t) = q * r (t) q nc n-r (t)
, where q * r (•) is a twice differentiable function such that

q * r (t k d ) = 0, q * r (t k d + ǫ) = q nc r (t k d + ǫ), q * r (t k d ) = 0, q * r (t k d + ǫ) = qnc r (t k d + ǫ) (26) and q * r (t k+ d ) = a > max 0, -A 1 (q) (λ d ) r (t k- d ) .
Proof: (i) The uniqueness of solution of the LCP [START_REF] Sadegh | Stability analysis of an adaptive controller for robotic manipulators[END_REF] guarantees that [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] 26) is imposed in order to assure the twice differentiability of the desired trajectory. Finally, straightforward computations show that

σ qr (t k d ) = q * r (t k+ d ) + A 1 (q) (λ d ) r (t k- d )
which means that the detachment is guaranteed and no other impacts occur when the desired acceleration satisfies q * r (t k+ d ) > max 0, -A 1 (q) (λ d ) r (t k- d ) .

Closed-loop stability analysis

In the case Φ = R n , the function V (t, s, q) in ( 13) can be used in order to prove the closed-loop stability of the system ( 11), ( 12) (see for instance [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF]). In the case studied here (Φ ⊂ R n ) the analysis becomes more complicated.

To simplify the notation V (t, s(t), q(t)) is denoted as V (t). In order to introduce the main result of this paper we make the next assumption, which is verified in practice for dissipative systems.

Assumption 1

The controller U t in [START_REF] Glocker | Concepts for Modeling Impacts without Friction[END_REF] assures that the sequence {t k ℓ } ℓ≥0 of the impact times possesses a finite accumulation point t k ∞ i.e. lim ℓ→∞

t k ℓ = t k ∞ < +∞, ∀k.
It is worth to precise that during the stabilization on the intersection of p surfaces Σ i we do not know which one and how many surfaces are stroked. However we assume that all the impacts are p ǫ -impacts in the sense of Definition 3.

Lemma 1 Consider the closed-loop system ( 10)-( 12) with (q * d ) p (•) defined on the interval [τ k 0 , t k 0 ] as in ( 14)- [START_REF] Menini | Asymptotic tracking of periodic trajectories for a simple mechanical system subject to nonsmooth impacts[END_REF]. Let us also suppose that condition b) of Proposition 1 is satisfied. The following inequalities hold:

||q(t k- 0 )|| ≤ V (τ k 0 ) γ 1 γ 2 , ||s(t k- 0 )|| ≤ 2V (τ k 0 ) λ min (M (q)) || q(t k- 0 )|| ≤ 2 λ min (M (q)) + γ 2 γ 1 V 1/2 (τ k 0 ) (27) Furthermore, if t k 0 ≤ τ k 1 one has ||(q d ) p (t k- 0 )|| ≤ ǫ + V (τ k 0 ) γ 1 γ 2 ||( qd ) p (t k- 0 )|| ≤ K + K ′ V 1/6 (τ k 0 ) ( 28
)
where ǫ is the real constant fixed in Definition 3 and K, K ′ > 0 are some constant real numbers that will be defined in the proof.

Proof: From (13) we can deduce on one hand that

V (t k- 0 ) ≥ γ 1 γ 2 ||q(t k- 0 )|| 2

INRIA

and on the other hand

V (t k- 0 ) ≥ 1 2 s(t k- 0 ) T M (q(t k- 0 ))s(t k- 0 ) Since condition b) of Proposition 1 is satisfied one has V (τ k 0 ) ≥ V (t k- 0 )
and the first two inequalities in (27) become trivial. Let us recall that s(t) = q(t) + γ 2 q(t)

which implies || q(t k- 0 )|| ≤ ||s(t k- 0 )|| + γ 2 ||q(t k- 0 )||.
Combining this with the first two inequalities in (27) we derive the third inequality in (27). For the rest of the proof we assume that

t k 0 ≤ τ k 1 . Therefore (q d ) p (t k- 0 ) = (q * d ) p (t k 0 ). It is clear that ||(q d ) p (t k- 0 )|| ≤ ||q p (t k- 0 )|| + ||q p (t k 0 )|| Taking into account that t k 0 is a p ǫ -impact (which means ||q p (t k 0 )|| ≤ ǫ), the first in- equality in (28) becomes obvious. Let us denote t ′ k = t k 0 -τ k 0 -δ τ k 1 -τ k 0 -δ ∈ [0, 1].
We recall here that τ k 0 was chosen such that ||q nc p (τ k 0 )|| ≤ ψ. From ( 14), ( 15) and the first inequality in (28), for i = 1, . . . , p one has

q i d (t k- 0 ) = (q i ) nc (τ k 0 ) + ϕV 1/3 (τ k 0 ) 2(t ′ k ) 3 -3(t ′ k ) 2 +(q i ) nc (τ k 0 ) ≤ ǫ+ V (τ k 0 ) γ 1 γ 2 It follows that 3(t ′ k ) 2 -2(t ′ k ) 3 ≥ (q i ) nc (τ k 0 ) -ǫ - V (τ k 0 ) γ1γ2 (q i ) nc (τ k 0 ) + ϕV 1/3 (τ k 0 ) For t > 0 one has 2t -t 2 ≥ 3t 2 -2t 3 , therefore 2t ′ k -(t ′ k ) 2 ≥ (q i ) nc (τ k 0 ) -ǫ - V (τ k 0 ) γ1γ2 (q i ) nc (τ k 0 ) + ϕV 1/3 (τ k 0 ) which means that (1 -t ′ k ) 2 ≤ V (τ k 0 ) γ1γ2 + ϕV 1/3 (τ k 0 ) + ǫ (q i ) nc (τ k 0 ) + ϕV 1/3 (τ k 0 ) Straightforward computations lead to | qi d (t k- 0 )| = 6((q i ) nc (τ k 0 ) + ϕV 1/3 (τ k 0 )) τ k 1 -τ k 0 -δ t ′ k -(t ′ k ) 2 Since t ′ k -(t ′ k ) 2 ≤ 1 -t ′ k and (q i ) nc (τ k 0 ) ≤ ψ, one arrives at (see Appendix 10.1 for details) | qi d (t k- 0 )| ≤ 6 √ ψǫ τ k 1 -τ k 0 -δ + 6 1 √ γ1γ2 + ϕ (ϕ + ψ) + ǫϕ τ k 1 -τ k 0 -δ V 1/6 (τ k 0 ) (29)
Therefore, the second inequality in (28) holds with

K = 6 √ pψǫ τ k 1 -τ k 0 -δ , K ′ = 6 √ p τ k 1 -τ k 0 -δ 1 √ γ 1 γ 2 + ϕ (ϕ + ψ) + ǫϕ
We now state the main result of this paper.

Theorem 1 Let Assumption 1 hold, e ∈ [0, 1) and (q * d ) p defined as in ( 14)-( 16). The closed-loop system ( 10)-( 12) initialized on Ω 0 such that V (τ 0 0 ) ≤ 1, satisfies the requirements of Proposition 1 and is therefore practically weakly stable with the closed-

loop state x(•) = [s(•), q(•)] and R = e -γ(t k f -t k ∞ ) (1 + K 1 + K 2 + ξ)/ρ where ρ = min{λ min (M (q))/2; γ 1 γ 2 }.
Proof: First we observe that conditions a) and d) of Proposition 1 hold when the hypothesis of the Theorem are verified. Thus in order to prove Theorem 1 it is sufficient to verify the conditions b), c) and e) of Proposition 1. It is noteworthy that during the transition phases I J k k some p ǫ -impact occurs (according to [START_REF] Dzonou | A sweeping process approach to inelastic contact problems with general inertia operators[END_REF] we have J ′ k = {1, . . . , p}). This means that we do not know which and how many of them are the constraints touched at each contact. However, in the neighborhood of the desired stabilization point situated on a surface of codimension p, only the corresponding p constraints enter the dynamics. In the sequel we shall also use the function V 1 (t, s) = 1 2 s(t) T M (q)s(t). b) Using that Ṁ (q) -2C(q, q)) is a skew-symmetric matrix, straightforward computations show that on R + \ k≥0 [t k 0 , t k f ] the time derivative of the Lyapunov function is given by

V (t) = -γ 1 s T s + 2γ 1 γ 2 qT q = -γ 1 || q|| 2 -γ 1 γ 2 2 ||q|| 2 On the other hand V (t) ≤ λ max (M (q)) 2 ||s|| 2 + γ 1 γ 2 ||q|| 2 ≤ γ -1 [γ 1 || q|| 2 + γ 1 γ 2 2 ||q|| 2 ]
where

γ -1 = max λ max (M (q)) 1 + 2γ 2 2γ 1 ; λ max (M (q))(γ 2 + 2) + 2γ 1 2γ 1 γ 2 > 0 Therefore V (t) ≤ -γ -1 V (t) on R + \ k≥0 [t k 0 , t k f ]. c) By definition V (t k- ℓ+1 )-V (t k+ ℓ ) = V 1 (t k- ℓ+1 )-V 1 (t k+ ℓ )+γ 1 γ 2 [(q T (t k- ℓ+1 ))q(t k- ℓ+1 )-(q T (t k+ ℓ ))q(t k+ ℓ )] (30)
On the other hand, straightforward computations show that

V 1 (t k- ℓ+1 ) -V 1 (t k+ ℓ ) = (t k ℓ ,t k ℓ+1 ) V1 (t)dt = γ 1 γ 2 (t k ℓ ,t k ℓ+1 ) s T p (t)(q * d ) p (t)dt -γ 1 (t k ℓ ,t k ℓ+1 ) s(t) T s(t)dt (31)
Furthermore,

(t k ℓ ,t k ℓ+1 ) s(t) T s(t)dt = (t k ℓ ,t k ℓ+1 ) || q(t)|| 2 + γ 2 2 ||q(t)|| 2 dt + γ 2 [(q T (t k- ℓ+1 ))q(t k- ℓ+1 ) -(q T (t k+ ℓ ))q(t k+ ℓ )] (32)
Therefore, inserting successively (32) in ( 31) and (31) in (30) one arrives at

V (t k- ℓ+1 ) -V (t k+ ℓ ) ≤ γ 1 γ 2 (t k ℓ ,t k ℓ+1 ) s T p (t)(q * d ) p (t)dt (33)

INRIA

In the sequel let us denote by S(q) the sum of all the components of the vector q.

Taking into account the definition ( 14) and the fact that (q d ) p and ( qd ) p are set to zero at t k+ 0 one obtains

(t k ℓ ,t k ℓ+1 ) s T p (t)(q * d ) p (t)dt = -ϕV 1/3 (τ k 0 ) • (t k ℓ ,t k ℓ+1 ) S( qp (t))dt+ γ 2 (t k ℓ ,t k ℓ+1 ) S(q p (t))dt Since ϕγ 2 V 1/3 (τ k 0 ) ≥ 0 and S(q p (t)) ≥ 0 it follows that t k ℓ+1 t k ℓ s T p (t)(q * d ) p (t)dt ≤ ϕV 1/3 (τ k 0 )[S(q p (t k ℓ )) -S(q p (t k ℓ+1 ))] Thus ℓ≥0 V (t k- ℓ+1 ) -V (t k+ ℓ ) ≤ γ 1 γ 2 ϕV 1/3 (τ k 0 )S(q p (t k 0 )) ≤ γ 1 γ 2 ϕV 1/3 (τ k 0 ) √ 3||q p (t k 0 )|| Recalling that t k 0 is an p ǫ -impact which means that ||q p (t k 0 )|| ≤ ǫ one obtains ℓ≥0 V (t k- ℓ+1 ) -V (t k+ ℓ ) ≤ K 1 V p1 (τ k 0 )
where

K 1 = √ 3γ 1 γ 2 ϕǫ > 0 and p 1 = 2
3 . e) First, let us compute the Lyapunov function's jumps at the instants t k ℓ , ℓ ≥ 1. Using the continuity of the real trajectory q(•) and the definition of the desired trajectory q d (•) on the I k phases (i.e. q d (t k+ ℓ ) = q d (t k- ℓ ), qd (t k+ ℓ ) = 0 = qd (t k- ℓ )) one gets

σ V (t k ℓ ) = V (t k+ ℓ ) -V (t k- ℓ ) = γ 1 γ 2 σ ||q|| 2 (t k ℓ ) + s T (t k+ ℓ )M ℓ s(t k+ ℓ ) -s T (t k- ℓ )M ℓ s(t k- ℓ ) 2 = T L (t k ℓ ) + γ 2 q(t k ℓ ) T M ℓ σ q (t k ℓ ) (34)
where M k ℓ denotes the inertia matrix M (q(t k ℓ )) and T L is the kinetic energy loss at the impact time t k ℓ . From equation ( 5) one has T L (t k ℓ ) ≤ 0 and equation (34

) becomes σ V (t k ℓ ) ≤ γ 2 q(t k ℓ ) T M ℓ σ q (t k ℓ ). Let us recall that M ℓ σ q (t k ℓ )
is the percussion vector (see [START_REF] Brogliato | Nonsmooth Mechanics[END_REF]). On the other hand in the X coordinates the percussion vector can be expressed as ∇F (X)λ. Writing the latter in the generalized coordinates introduced in Section 3 one obtains M ℓ σ q (t k ℓ ) = D T λ. In other words the generalized coordinates introduced in Section 3 coincide with the so called quasi-coordinates and the vector qtang is in fact

qn-m (i.e. σ q (t k ℓ ) = σ qm (t k ℓ ) 0 n-m
where 0 n-m denotes the nm vector with all its components equal zero). Therefore

σ V (t k ℓ ) ≤ γ 2 q(t k ℓ ) T M ℓ σ q (t k ℓ ) = γ 2 q p (t k ℓ ) T λ p = 0 (35)
where we have used (q d ) p (t k+ ℓ ) = 0 = (q d ) p (t k- ℓ ) and the last equality is stated using the complementarity relation entering the dynamics.

The Lyapunov function's jump corresponding to the first impact of each cycle can be computed as:

σ V (t k 0 ) = V (t k+ 0 ) -V (t k- 0 ) = γ 1 γ 2 σ ||q|| 2 (t k 0 ) + s T (t k+ 0 )M 0 s(t k+ 0 ) -s T (t k- 0 )M 0 s(t k- 0 ) 2 (36) • It is clear that t k 0 > τ k 1 implies q d (t k+ 0 ) = q d (t k- 0 )
and qd (t k+ 0 ) = 0 = qd (t k- 0 ). Thus, the computations for t k ℓ , ℓ ≥ 1 hold also for t k 0 .

• If t k 0 ≤ τ k 1 one has (q d ) p (t k- 0 ) = (q d ) p (t k+ 0 ) = 0 and ( qd ) p (t k- 0 ) = ( qd ) p (t k+ 0 ) = 0. Then the initial jump of each cycle is given by:

σ V (t k 0 ) = qd (t k- 0 ) T M 0 q(t k- 0 ) + γ 2 2 2 q(t k+ 0 ) T M 0 q(t k+ 0 ) -q(t k- 0 ) T M 0 q(t k- 0 ) + γ 2 q(t k+ 0 ) T M 0 q(t k+ 0 ) -q(t k- 0 ) T M 0 q(t k- 0 ) - 1 2 qd (t k- 0 ) T M 0 qd (t k- 0 ) + T L (t k 0 ) (37)
Since T L (t k 0 ) ≤ 0 the equation (37) rewrites as:

σ V (t k 0 ) ≤ λ max (M (q)) γ 2 ||( qd ) p (t k- 0 )|| • ||q(t k- 0 )|| + || q(t k- 0 )|| • ||(q d ) p (t k- 0 )|| + γ 2 2 2 ||q p (t k 0 )|| 2 + ||q p (t k- 0 )|| 2 + 2||(q d ) p (t k- 0 )|| • ||q n-p (t k- 0 )|| + 1 2 ||( qd ) p (t k- 0 )|| 2 + ||( qd ) p (t k- 0 )|| • || q(t k- 0 )|| (38)
Obviously || q(t k- 0 )|| = || q(t k- 0 )+ ( qd ) p (t k- 0 )|| and Lemma 1 combined with V (τ k 0 ) < 1 yields

|| q(t k- 0 )|| ≤ K + 2 λ min (M ) + γ 2 γ 1 + K ′ V 1/6 (τ k 0 ) Therefore σ V (t k 0 ) ≤ K 2 V p2 (τ k 0 ) + ξ where p 2 = 1 6 , ξ = 3 2 K 2 + γ 2 ǫK + γ 2 2 ǫ 2
2 and

K 2 = λ max (M (q)) 3KK ′ + 3 2 (K ′ ) 2 + γ 2 2 γ 2 γ 1 + 2γ 2 λ min (M (q))γ 1 + 2γ 2 γ 1 + (K ′ + K) 3 γ 2 γ 1 + 2 λ min (M (q)) + ǫγ 2 2 γ 2 γ 1 + 2 λ min (M (q)) + K ′
Defining α : R + → R + , α(ω) = ρω 2 we get α(0) = 0 and α(||[s(t), q(t)]||) ≤ V (t, s, q). Thus, Proposition 1 also yields

R = α -1 (e -γ(t k f -t k ∞ ) (1 + K 1 + K 2 + ξ)) = e -γ(t k f -t k ∞ ) (1 + K 1 + K 2 + ξ)/ρ
which ends the proof.

INRIA 8 Illustrative example

The main issues of the control scheme proposed in this paper are emphasized simulating the behavior of a planar two-link rigid-joint manipulator in presence of one unilateral constraint.

Dynamics Equation based on Lagrangian formulation

First, we will derive the dynamics of a two-link manipulator moving in the XOY plane.

In order to accomplish our goal let us introduce the following notations (see Figure 5):

• θ i -the joint angle of the joint i,

• m i -the mass of link i,

• I i -the moment of inertia of link i about the axis that passes through the center of mass and is parallel to the Z axis,

• l i -the length of link i,

• g -the gravitational acceleration. For the sake of simplicity we assume that the center of mass of each link is right at the middle of the link. We will also assume that the gravitational force acts in the negative direction of the Y axis.

First, let us choose θ = (θ 1 , θ 2 ) T as the generalized coordinates and denote T the kinetic energy, P the potential energy and L = T -P the Lagrangian function. Using the above notations, the potential and the kinetic energies of the system can be expressed as follows:

P = 1 2 m 1 gl 1 sin θ 1 + m 2 g l 1 sin θ 1 + l 2 2 sin(θ 1 + θ 2 ) , T = m 1 l 2 1 θ1 2 8 + m 2 2 l 2 1 θ2 1 + l 2 2 4 ( θ1 + θ2 ) 2 + l 1 l 2 ( θ2 1 + θ1 θ2 ) cos θ 2 + I 1 θ1 2 + I 2 ( θ1 + dotθ 2 ) 2
2 .

Substituting P and T in the formula of L the dynamical equation of the two-links manipulator:

d dt ∂L ∂ θ - ∂L ∂θ = 0
can be rewritten as

M (θ) θ + C(θ, θ) θ + G(θ) = 0 (39)
where M, C ∈ R 2×2 , G ∈ R 2 have the same meaning as in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF] and they are explicitly given by:

M = » M11 M12 M21 M22 - , 8 > > > > > > < > > > > > > : M11 = m1l 2 1 4 + m2 " l 2 1 + l 2 2 4 l1l2 cos θ2 « + I1 + I2 M12 = M21 = m2l 2 2 4 + m2l1l2 2 cos θ2 + I2 M22 = m2l 2 2 4 + I2 C = » C11 C12 C21 C22 - , 8 > > > > > > < > > > > > > : C11 = -m2l1l2 θ2 sin θ2 C12 = - m2l1l2 2 θ2 sin θ2 C21 = m2l1l2 2 θ1 sin θ2 C22 = 0 G = » G1 G1 - , 8 > < > : G1 = g 2 [l1(2m1 + m2) cos θ1 + m2l2 cos(θ1 + θ2)] G2 = m2gl2 2 cos(θ1 + θ2)
The admissible domain is the upper half plane y ≥ 0 and the unconstrained desired trajectory q nc (•) is given by a circle that violates the constraint. Precisely, the end effector must follow a half-circle, stabilize on the constraint (y = 0) and move on the constraint until the point where the circle q nc (•) re-enters the admissible domain. The lengths l 1 , l 2 of the manipulator's links are set to 0.5m, and their masses m 1 , m 2 are set to 1kg. The impacts are imposed using the parameter ϕ = 100 in ( 14)- [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF].

Some remarks concerning the simulation scheme

The numerical simulations are done with the Moreau's time-stepping algorithm of the SICONOS software platform [START_REF] Acary | An introduction to SICONOS[END_REF]. The choice of a time-stepping algorithm was mainly dictated by the presence of accumulations of impacts which render the use of eventdriven methods difficult [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF].

Let us set e = 0.7, γ 1 = 8, γ 2 = 7, 10 seconds the period of each cycle and 30 seconds the final simulation time. First, let us point out (Figure 6 (left)) the influence of ψ (i.e.

INRIA the choice of τ k 0 ) on the degree of deformation of the real trajectory w.r.t. the desired unconstrained one. As we have pointed out in Section 4 the deformation gets smaller when ψ > 0 decreases. It is noteworthy that the tangential approach corresponding to ψ = 0 lacks of robustness and is unreliable due to the nonzero initial tracking errors. For ψ = 0.01 in Figure 6 (right) we illustrate the trajectory of the system during the stabilization on ∂Φ = {(x, y) | y = 0}. The switches of the controller during the first cycle are depicted in Figure 7 (left). Clearly since the velocity jumps, the controller jumps as well. The Figure 7 (right) presents the variation of the contact force λ. Precisely, one sees that λ remains 0 during the free motion phases and it points out each impact during transients (better seen on I 0 since the impacts are more violent). The contact force λ is designed as a decreasing linear function during constrained motion phases Ω 2k+1 in order to allow a smooth detachment at the end of these phases. It is worth to mention that the magnitude of λ depends indirectly on V (τ k 0 ). Precisely, when V (τ k 0 ) approaches zero the system tends to a tangential stabilization on ∂Φ which implies larger values of t k 0 and consequently smaller length of [t k f , t k d ] and smaller magnitude of the contact force measured by λ (see Proposition 4).

Figure 8 shows that the tracking error described by the Lyapunov function rapidly decreases and remains close to 0. In other words the practical weak stability is guaranteed. On the zoom made in Figure 8 one can also observe the behavior of V (•) during the stabilization on ∂Φ, that is an almost decreasing function.

V (t 1+ 0 ) V (t 1+ 2 ) V (t 1- 1 ) V (t 1+ 1 ) V (t 1- 2 ) V (t 1- 0 )

The influence of the time-step on the closed-loop dynamics

It is noteworthy that the simulation results do not depend essentially on the time-step chosen for the scheme (since the Moreau's time-stepping algorithm converges [START_REF] Dzonou | A sweeping process approach to inelastic contact problems with general inertia operators[END_REF]) but, a smaller time-step allows to capture more precisely the behavior of the system. On the other hand, as can be seen in figure 9 the real trajectory and the lengths of each transition phase are almost unchanged starting with a sufficiently small time-step (h = 10 -3 ). We do not insist too much on the simulation results during free-motion phases since the smoothness of the system is guaranteed on these phases and the behavior of the system is clear. The most interesting phases from the numerical point of view are the transition (accumulation of impacts) phases. It is worth to clarify that the number of impacts during the transition phases is not so important and the major issue is the finiteness of these phases. To be more clear, in the next tables we summarize some numerical results when ψ = 0.01, e = 0.7, γ 1 = 35, γ 2 = 20. The period of each cycle is set to 5 seconds and the final simulation time rest 30 seconds. First, one can see that the length of the transition phase I 0 with respect to the time-step h does not vary significantly when the time-step decreases. Let us also denote by CP U the computing time necessary for the simulation (using an Intel(R) Core(TM)2 CPU 6300 1.86GHz) of one cycle. The evolutions of the number of impacts n i w.r.t. the restitution coefficient e and the time-step h are quite different. As expected, n i becomes larger when the restitution coefficient increases. Also, one can see that the accumulation of impacts can be captured with a higher precision when the time-step becomes smaller.

e \ h 10 -3 s 10 -4 s 10 -5 s 10 -6 s 0.2 In conclusion, reliable simulations with a reasonable CP U time can be performed with the Moreau's time-stepping scheme of the SICONOS platform, with a time-step h = 10 -4 s.

n i = 3 n i = 5 n i = 6 n i = 8 0.5 n i = 6 n i = 9 n i = 12 n i = 16 0.7 n i = 9 n i = 16 n i = 23 n i =

The influence of controller parameters on the closed-loop dynamics

In this paragraph we present several simulations with different values of the controller parameters in order to see their influence on the closed-loop dynamics. Some intuitive explanation will join the numerical results presented in the following. It is worth to precise that we keep for the next simulations a period of 5 seconds for each cycle and we fix the time-step at the value of 10 -4 seconds.

In order to explain more easily the results of this section we point out from the beginning that γ 1 is the coefficient of the velocity error and γ 1 γ 2 is the coefficient of the position error entering both in the controller and the Lyapunov function.

First one considers a fixed value of γ 1 = 25 and one studies the behavior of the closed-loop dynamics w.r.t. γ 2 variation. The main influence in this case can be seen during the transition phases I k . More precisely, the Figure 10 shows that diminishing the value of γ 2 the transition phases get larger. From a mechanical point of view decreasing γ 2 one decreases the stiffness of the spring entering the controller and therefore we get more jumps, higher first jump (see Figure 10 (right)) and longer transition phases. Obviously, this may lead to the system's instability (as can be seen in Figure 10 (left) for γ 2 = 2). The above discussion is also illustrated in the next table where we denote by H the height of the first jump. Next one considers a fixed value of γ 1 γ 2 = 800 and we repeat the simulation for different value of γ 1 . From mechanical point of view, this means that we fix the spring stiffness entering the controller and we point out the effect of the damper gain on the system dynamics. Precisely, decreasing the value of γ 1 we detect an increasing tracking error, see figure 11. In other words, the real trajectory approaches slower the desired trajectory. In the sequel we want to emphasize the behavior of the system during each phase of the motion. Let us consider the spring stiffness γ 1 γ 2 = 5 and the damper gain γ 1 = 10.

INRIA

Taking into account the response of the system to the variation of the spring stiffness and damper gain entering the controller this choice leads to one of the worst possible situation (small convergence speed, high jumps and long transients). In order to compensate the bad choice of controller gains we consider a low value of the restitution coefficient e = 0.3 and, since the initial tracking errors are relatively big, the impacts are imposed using a smaller parameter ϕ = 10 in ( 14)- [START_REF] Mcclamroch | Feedback stabilization and tracking of constrained robots[END_REF]. However, the high jumps and the small convergence speed allows to better identify each particular phase of the motion (see Figure 12 (left)). On the other hand, as depicted in Figure 12 (right), the accumulation of impacts phases are too long and the end-effector cannot accomplish its task in the fixed period set to 5 seconds. It is noteworthy that large values of spring stiffness and dumper gain improve the tracking control leading to an "almost tangential" contact. Nevertheless, as pointed out before, a tangential contact is not desired due to its unreliability.

Conclusions

In this paper we have proposed a methodology to study the tracking control of fully actuated Lagrangian systems subject to multiple frictionless unilateral constraints. The main contribution of the work is twofold: first, it formulates a general framework and second, it provides a complete stability analysis for the class of systems under consideration. It is noteworthy that even in the case of only one frictionless unilateral constraint the paper presents some improvements with respect to the existing works in the literature. Precisely, the stability analysis result is significantly more general than those presented in [START_REF] Bourgeot | Tracking control of complementarity lagrangian systems[END_REF] and [START_REF] Brogliato | On the control of finite-dimensional mechanical systems with unilateral constraints[END_REF] and, each element entering the dynamics (desired trajectory, contact force) is explicitly defined. Numerical simulations are done with the SICONOS software platform [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Acary | An introduction to SICONOS[END_REF] in order to illustrate the results.

  i∈I

F

  i (X(t)) = 0 where I ⊂ 1, . . . , m, card(I) = p.

Figure 1 :

 1 Figure 1: One degree of freedom: unitary mass dynamics

Figure 2 :

 2 Figure 2: Typical evolution of the Lyapunov function during one cycle of a weakly stable system.

Figure 3 :

 3 Figure 3: The closed-loop desired trajectory and control signals

  ) r which leads to the second part of definition (23). Furthermore, replacing (λ d ) p-r in (25) we get λ p-r = C 2 and b r + A 1 (λ-λ d ) r + A 2 (λ-λ d ) p-r ≥ 0 yields the first part of definition (23). To conclude, the solution of the LCP (24) is λ p = 0 C 2 ∈ R p and (λ d ) p is defined by (23).

Figure 5 :

 5 Figure 5: Two-link planar manipulator

Figure 6 :

 6 Figure 6: Left: The influence of ψ on the real trajectory's deformation for controller's gains set to γ1 = 8, γ2 = 7. Right: The trajectory of the end-effector on the transition phases when ψ = 0.1.

Figure 7 :

 7 Figure 7: Left: The switching controller on the first cycle; Right: Variation of the contact force λ

Figure 8 :

 8 Figure 8: Variation of the Lyapunov function for γ1 = 8, γ2 = 7; Zoom: Variation of the Lyapunov function during the phase I0

h 10 - 3 sINRIAFigure 9 :

 1039 Figure 9: Left: The variation of q1 during transition phase for h ∈ {10 -3 , 10 -4 , 10 -5 , 10 -6 }; Right: Zoom at the end of transition phase.

  29 0.9 n i = 23 n i = 40 n i = 64 n i = 81 0.95 n i = 32 n i = 67 n i = 108 n i = 161 However, a larger number of captured impacts does not change the global behavior of the simulated system and the transition phase ends almost in the same moment when h varies, see λ[I 0 ] in the first table.

Figure 10 :

 10 Figure 10: Left: The variation of q1(t) during 6 cycles for γ2 ∈ {2, 30, 55}; Right: The variation of q1(t) during I0 for γ2 ∈ {30, 55}.

Figure 11 :

 11 Figure 11: Smaller values of γ1 lead to larger tracking errors.

Figure 12 :

 12 Figure 12: Left: The system trajectory in the xOy-plane for γ1 = 10 and γ2 = 0.5; Right: The variation of q1(t) for γ1 = 10 and γ2 = 0.5.

  and (23) agree if C 1 < 0. In other words, replacing C 1 by C 1 (tt k d ) in (23) we assure a constrained motion on [t k f , t k d ) and the necessary conditions for detachment on [t k d , t k d + ǫ).

	(ii) Obviously (

RR n°6483

Morȃrescu & Brogliato

In particular it is implicitly assumed that the function F i (•) in (1) are linearly independent INRIA

Appendix 10.1 Computation details for inequality (29)

In order to prove inequality (29) we recall that

and

which leads to (29).


Justyna Nowakowska

Catherine Bastien

Brigitte Musch

Catherine Catherine Bastien Bastien

Genetic variability of Polish Scots pine provenances assessed with microsatellite markers Genetic variability of Polish Genetic variability of Polish Scots pine provenances Scots pine provenances assessed with assessed with microsatellite markers microsatellite markers

Results Results

MapInfo

  Every peak was examined for appropriate threshold scoring Every peak was examined for appropriate threshold scoring Every peak was examined for appropriate threshold scoring

  human activity during ages

Localization of the Scots pine provenances: 204. Augustów 208. Białowieża 504. Bolesławiec 306. Bolewice 105. Bytów 101. Goleniów 307. Gubin 606. Janów Lubelski 606. Józefów 654. Kluczbork 602. Kozienice 352. Krucz 303. Lipka 403. Łochów 404. Międzyrzec 107. Międzyzdroje 106. Miłomłyn 302. Niedźwiady 607. Niepołomice 205. Nowe Ramuki 401. Nowogród 104. Smolarz 601. Spała 206. Strzałowo Localization of the Scots pine provenances: 204. Augustów 208. Białowieża 504. Bolesławiec 306. Bolewice 105. Bytów 101. Goleniów 307. Gubin 606. Janów Lubelski 606. Józefów 654. Kluczbork 602. Kozienice 352. Krucz 303. Lipka 403. Łochów 404. Międzyrzec 107. Międzyzdroje 106. Miłomłyn 302. Niedźwiady 607. Niepołomice 205. Nowe Ramuki 401. Nowogród 104. Smolarz 601. Spała 206. Strzałowo 405. Strzelce 207. Supraśl 501. Syców 108. Wejherowo 305. Woziwoda 402. Wyszków 405. Strzelce 207. Supraśl 501. Syców 108. Wejherowo 305. Woziwoda 402. Wyszków -Provenance -Maternal Seed Zone Comment citer ce document :

Allele size Frequency Allele frequency > 0.03 (TG)17(AG)21 Results Results MapInfo Professional™ 4.1 MapInfo Professional™ 4.1 Geographical distribution of the most most frequent frequent alleles of alleles of 7.14 locus 7.14 locus Results Results Allele frequency distribution: SPAC 12.5 microsatellite loci (N=1385, 38 alleles) Allele frequency distribution: SPAC 12.5 microsatellite loci (N=1385, 38 alleles) Allele frequency > 0.04

		0.08 0.08 0.12 0.12
		0.07 0.07																			(GT)20(GA)10
		0.1 0.1
		06 0.06
	Frequency Frequency Frequency	0.03 0.04 0.08 0.08 0.05 0.06 0.06
		0.04 0.04
		0.02
		0.02 0.02
		0.01
		0 0 0
		173 173 119 119	175 175 121 177 177 123 179 179 125 181 181 127 183 183 121 123 125 127 129 185 185 129 131 187 187 131 133 189 189 133	191 191 135 135	193 193 137 137	195 195 139 139	197 197 141 141	199 199 143 143	201 201 145 145	203 203 147 147	205 205 149 149	207 207 151 151	209 209 153 153	211 211 155 155	213 213 157 157	215 215 159 159	217 217 161 161	219 219 163 163	221 221 165 165	223 223 167 225 225 169 227 227 171 229 229 173 230 230 167 169 171 173 175 232 232 175 177 234 234 177 179 236 236 179 181 238 238 181	240 240 183 183	242 242 185 185	244 244 191 191	246 246 193 193	260 260 195 195	262 262 197 197
													Allele size Allele size Allele size

Allele size Frequency (AT)10 Results Results MapInfo Professional™ 4.1 MapInfo Professional™ 4.1 Geographical distribution of the most most frequent frequent alleles of alleles of SsrPt- ctg4363 locus locus Provenance: Seed micro -zone:

								0.053			0.043		0.037	0.030	0.022			0.008		0
																									Supraśl	2 0 7
																									Syców	501
																									Woziwoda	305
																									Spała	6 0 1
																									Gubin	307
																									Janów Lub.	606
																									Kluczbork	654
																									Niedźwiady	302
																									Goleniów	101
		0.4 0.4																							Lipka Augustów	303 204
		0.35 0.35																							Strzałowo Józefów	206 606
		3 0.3																							Wejherowo Strzelce	108 405
	Frequency	0.2 0.25 0.15																							Wyszków Miłomłyn Smolarz Nowe Ramuki Nowogród Łochów	402 106 104 205 401 403
		0.1																							Bolesławiec Krucz	504 352
																									Międzyrzec	404
		0.05																							Bytów	105
																									Kozienice	602
		0	89 89	91 91	94 94	96 96	98 98	99 99	101 101	103 103	105 105	107 107	109 109	111 111	113 113	115 115	117 117	119 119	121 121	126 126	128 128	134 134	138 138	140 140	148 148 Bolewice 150 150 152 152 Międzyzdroje 166 166	306 107
														Allele size								Niepołomice	607
																									Białowieża	2 0 8

Neighbor-joining tree of Scots pine provenances in Poland based on microsatellite data Neighbor-joining tree of Scots pine provenances in Poland based on microsatellite data POPGENE v 1.32. POPGENE v 1.32. 108 405 352 101 107 105 305 303 302 104 307 504 501 607 Jan Józ 208 207 205 402 403 602 606 404 654 -Provenance -Maternal Seed Zone 107 Geographical distribution of the Scots pine linkage groups in Poland 204. Augustów 208. Białowieża 504. Bolesławiec 306. Bolewice 105. Bytów 101. Goleniów 307. Gubin 606. Janów Lubelski 606. Józefów 654. Kluczbork 602. Kozienice 352. Krucz 303. Lipka 403. Łochów 404. Międzyrzec 107. Międzyzdroje 106. Miłomłyn 302. Niedźwiady 607. Niepołomice 205. Nowe Ramuki 401. Nowogród 104. Smolarz 601. Spała 206. Strzałowo Geographical distribution of the Scots pine linkage groups in Poland 204. Augustów 208. Białowieża 504. Bolesławiec 306. Bolewice 105. Bytów 101. Goleniów 307. Gubin 606. Janów Lubelski 606. Józefów 654. Kluczbork 602. Kozienice 352. Krucz 303. Lipka 403. Łochów 404. Międzyrzec 107. Międzyzdroje 106. Miłomłyn 302. Niedźwiady 607. Niepołomice 205. Nowe Ramuki 401. Nowogród 104. Smolarz 601. Spała 206. Strzałowo 405. Strzelce 207. Supraśl 501. Syców 108. Wejherowo 305. Woziwoda 402. Wyszków 405. Strzelce 207. Supraśl 501. Syców 108. Wejherowo 305. Woziwoda 402. Wyszków

Acknowledgments Acknowledgments This work represents a part of genetic diversity research on conifer species financed by Polish General Directorate of the State Forests and French Ministry (DSUR-NGE-4C1-701 Project). Special thanks are addressed to Dr Sylvie Oddou (INRA Avignon, France) for valuable discussions and to Mrs. Bénédicte LeGuerroué and Mrs. Vanina Benoit -(INRA Orleans, France) for help in the laboratory task. This work represents a part of genetic diversity research on conifer species financed by Polish General Directorate of the State Forests and French Ministry (DSUR-NGE-4C1-701 Project). Special thanks are addressed to Dr Sylvie Oddou (INRA Avignon, France) for valuable discussions and to Mrs. Bénédicte LeGuerroué and Mrs. Vanina Benoit -(INRA Orleans, France) for help in the laboratory task.


domain: sdv sde

text:
T W Ellis

S Leguédois

) P B Hairsine

D J Tongway

Performance of a tree belt for capturing overland flow from agricultural land -implications for design

Introduction

There is growing interest in the strategic placement of tree belts in agricultural land to increase water use in areas with dryland salinity and waterlogging (Stirzaker et al., 2002;[START_REF] Ellis | An ecological optimality approach for predicting deep drainage from tee belts in alley farms in water limited environments[END_REF]. Depending on the soil type, tree belts placed on hillslopes can capture excess water moving from upslope by either intercepting shallow sub surface lateral flow (Stirzaker et al., 2002), or by capturing overland flow. These processes can also be utilised to improve the productivity of tree crops in arid to semiarid areas by increasing the water supply to them. On soils with a low infiltration capacity, the principle process for delivery of excess water from upslope is by overland flow. In these areas, the introduction of tree belts into agricultural landscapes is expected to result in a slow change in soil surface attributes that improve water capture and storage processes. We describe a field experiment, using a rainfall simulator, designed to investigate the performance of a tree belt for the interception of overland flow from grazed pasture. Our results are used to highlight the main parameters to be considered when designing tree belts specifically for this purpose.

Experimental design

The experimental site was located on red duplex soil (bulk density 1.35 -1.75 g cm -3 ) with 12 o slope on a grazing property near Boroowa (-34 o 22' S 148 o 42' E), New South Wales, Australia. The slope comprised 28 m of grazed pasture draining into 12 m of tree belt at the lower end. The 15 year-old Acacia spp. and Calistemon spp. tree belt, aligned perpendicular to the slope, was originally planted by the landholder for stock shelter and biodiversity reasons. Stock were excluded from the tree belt by a fence located near the drip line of the present canopy. A bare soil area associated with tree-pasture competition and stock traffic extended 5 to 7 m beyond the fence line into the pasture.

Details of the construction of the rainfall simulator used in this experiment are described in [START_REF] Motha | Tracer properties of eroded sediment and source material[END_REF] and [START_REF] Wilson | Effects of logging and fire on runoff and erosion on highly erodible granitic soils in Tasmania[END_REF].Three simulated rainfall events were applied: 45 mm hr -1 for 13 min duration; 45 and 75 mm hr -1 for 30 min duration, over the 40 x 15 m experimental area, split longitudinally so that duplicate pasture-tree belt sequences were measured. One of the sequences was retained as a single plot and measured for its combined behavior of pasture plus tree belt. The other sequence was divided into two separate plots: pasture only; and tree belt only, respectively.

For each plot, overland flow was directed into a flume and measured every 3 minutes during the rainfall simulations and during flow recession following rainfall cessation. Changes in flow patterns (spatial and temporal) and subsequent changes in surface micro topography were recorded for all events. Together with observations, flow hydrographs were used to estimate surface water storage (due to micro topography) during rainfall and recession periods.

Main results

• Overland flow from the first event (45 mm hr -1 ; 13 minutes) was not sufficient to 'break through' the 12 m section of tree belt, and all overland flow leaving the pasture infiltrated into the soil. This represented a 37% increase in rainfall depth within the tree belt.

• Approximately 50% of the overland flow generated by the second event (45 mm hr -1 ; 30 min) was intercepted and stored by the tree belt. This represented a 36% increase in rainfall depth within the tree belt.

• There was no measurable capture of overland flow resulting from the third event (75 mm hr -1 ; 30 min).

• The bare soil area was observed to dominate the runoff generation from the pasture component in the early phases of all three rainfall events, and infiltration in this zone was close to zero (as determined by 'before' and 'after' gravimetric soil water profiles).

• A significant backwater formed at the interface between pasture and tree belt where there was an abrupt increase in surface roughness. Sediment accumulation in this zone immediately upslope of the tree belt. The backwater was associated with small litter dams that were present throughout the tree belt. These structures were composed of pasture and tree litter and served to spread the surface flow and maximise the surface area in the tree belt that was available for infiltration.

• There appeared to be no appreciable differences in soil structure, between the pasture and tree belt plots, likely to be associated with interception and storage of water.

Conclusions

The tree belt significantly reduced the amount of overland flow leaving the experimental area, and the water captured represented a significant (greater than one third) increase in water supply to the trees. It is likely that tree belt designs can be optimised for this purpose and issues to be considered in this regard include:

1. The runoff coefficient of the upslope zone and how this would vary with management (e.g. crop or pasture).

2. The ratio of upslope length (crop/pasture) to downslope tree belt length.

3. Soil micro topography within the tree belt.

4. The nature, decomposition character and longevity of tree and pasture litter within the tree belt.

5. The water storage capacity within the tree belt.

6. Any prospect for soil structural changes within the tree belt (e.g. macropores) increasing the entry and storage of overland flow. These changes may be associated with older or native remnant tree belts.

The experiment highlighted the significant role of the bare soil zone, located upslope of the tree belt, as a source, rather than a sink, for overland flow. A modified tree belt design, where the fence is placed beyond this zone, may increase the potential for capture of overland flow by ground layer plants and accumulated litter.

In the absence of an extensive data set of the physical parameters required to describe the hydrological processes occurring within tree belts during overland flow events, it may be possible to develop simple empirical rules relating soil water capture potential to readily available attributes such as 2, 3 and 4 above.

Acknowledgements

We thank David Marsh for allowing us access to his property for the rainfall simulation experiment; Jim Brophy for his significant efforts in site preparation, erection and operation of the rainfall simulator; and work colleagues from CSIRO's Integrated Catchment Management Research Directorate for their assistance with monitoring and sample collection.