Harshal B Oza

Vincent Acary
email: vincent.acary@inrialpes.fr

Yury V Orlov
email: yorlov@cicese.mx

Sarah K Spurgeon
email: s.k.spurgeon@kent.ac.uk

Bernard Brogliato
email: bernard.brogliato@inrialpes.fr

Finite time tracking of

come

I. INTRODUCTION

Dynamical systems with non-smooth vector fields and with unilateral constraints have been studied in various disciplines of science and engineering such as electrical engineering and robotics amongst others [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF]. Establishing conditions for stability of such systems is a theoretical challenge due to the likely destabilizing effects of the unilateral constraints, resulting in state jumps. In this paper, planar controllable linear and feedback linearizable non-linear systems with unilateral constraints are studied. Unilaterally constrained systems can be studied as linear or non-linar complementarity systems [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF]. The closed-loop dynamics of these systems are strongly affected by hard non-linearities characterized by both the nature of the contact of two bodies as well as collision, with the constraint resulting in a restitution in one of the states of a planar system. The existing literature is abundant on the topic of stability of systems with jumps (see [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF], [START_REF] Morȃrescu | Trajectory tracking control of multiconstraint complementarity lagrangian systems[END_REF] and the references therein).

The main focus of this paper is on (i) finite time tracking of trajectories of the unilaterally constrained system to the desired periodic trajectories by using a second order sliding mode controller together with tuning of the controller parameters and (ii) digital implementation of the same. The theoretical motivation to study such a problem stems from several viewpoints. Firstly, asymptotic and finite time stability when impacts have a finite accumulation point have been studied for restitution coefficient in the interval [0, 1) [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF], [START_REF] Oza | Finite time stabilization of perturbed double integrator with jumps in velocity[END_REF]. However, a priori guarantee and tuning to achieve finite time convergence for the tracking of periodic trajectories is not studied for the general case where restitution in the plant state is beyond the interval [0, 1). Furthermore, robust finite time stabilization is theoretically an equally difficult challenge when compared to its asymptotic counterpart [START_REF] Morȃrescu | Trajectory tracking control of multiconstraint complementarity lagrangian systems[END_REF] provided that the effects of dissipative impacts tend to destabilize the system.

Secondly, research on robustness of the controller to disturbances at the time of state jumps is not rigourous. The present paper studies a similar class of systems to that presented in recent research [START_REF] Zheng | Finite time observer-based control of linear impulsive systems with persistently acting impatcs[END_REF] where a finite time observer is proposed to reject disturbances at the time of impacts. The presented result utilises tuning for a second order sliding mode (SOSM) controller [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF] such that it is robust to persisting disturbances in both the continuous and discrete-event phases of the trajectories, while requiring only upper bounds on the uncertainties to cause the tracking errors to converge to zero in finite time between successive impacts.

Finally, existing results on numerical methods for sliding mode synthesis [START_REF] Acary | Implicit euler numerical scheme and chattering-free implementation of sliding mode systems[END_REF], [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF] can be directly utilised for the tracking control of a unilaterally constrained system in the presence of impacts. A combination of the tuning of a robust SOSM synthesis that is based on the continuous time setting with the digital implementation, which matches the behaviour of the discrete time closed-loop system to that of the continuous time system at least in the case of no uncertainty, is theoretically novel for the underlying impact system.

From a practical viewpoint, the motivation stems from the need to propose tuning rules and digital implementation of SOSM for the presented class of systems. Following the recent advances in (i) Lyapunov methods for higher order sliding mode controllers and observers [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF], [START_REF] Santiesteban | Finite-time convergence analysis for "twisting" controller via a strict Lyapunov function[END_REF], [START_REF] Oza | Lyapunov based settling time estimate and tuning for twisting controller[END_REF], (ii) finite time convergence and tuning [START_REF] Oza | Lyapunov based settling time estimate and tuning for twisting controller[END_REF] and (iii) numerical schemes [START_REF] Acary | Implicit euler numerical scheme and chattering-free implementation of sliding mode systems[END_REF], [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF], a natural next step is to provide similar tuning and digital implementation for this problem. This problem is relevant in the area of biped robots and in general in fully actuated robots with unilateral constraints and thus has practical significance

The main contribution of the presented result is twofold. Firstly, the tuning rules for the controller are given for the presented tracking problem with an a priori guarantee of attaining a pre-specified settling time. The method differs from existing results [START_REF] Laghrouche | Higher order sliding mode control based on integral sliding mode[END_REF] that depend on open-loop optimal control. The drawback of the presented method is that it is more conservative. Nevertheless, it gives theoretical starting values for tuning that guarantee finite time tracking of the states of the hybrid system in the presence of disturbances in both continuous and discrete event (jump) phases. The second contribution is the 'numerical chatter-free' digital implementation of this SOSM controller for the tracking problem of a unilaterally constrained system resulting in reduced chattering in the presence of the aforementioned disturbances [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF].

II. PROBLEM STATEMENT Let the continuous system dynamics be given as follows:

ẋ1 = x 2 ẋ2 = u + ω(t) (x 1 , x 2 ) ∈ S (1)
where x 1 is position, x 2 is velocity, u is the control law to be synthesized, ω(t) is an external disturbance, and the constraint surface S is given by

S = {(x 1 , x 2 ) : x 1 = 0}. ( 2
)
The constraint surface gives rise to collision and induces resets in velocity x 2 (t). Let the velocity restitution (or impact event due to the unilateral constraint) be given according to the following Newtonian Reset Map:

x 2 (t + k ) = -ēx 2 (t - k ) + ω ∆ (x 1 , x 2 ) ∈ S (3)
where the notations x(t + k ) and x(t - k ) denote the right limit lim tt k k ) and the left limit lim

tt k x(t k ) respectively, ē ∈ [0, 1)
is the restitution parameter responsible for the energy absorbtion [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF], t k , k = 1, 2,. .., ∞ are the time instants at which these impulses occur and ω ∆ characterises the restitution uncertainty. Such an uncertainty ω ∆ represents a lack of complete knowledge of the linear reset map [START_REF] Oza | Finite time stabilization of perturbed double integrator with jumps in velocity[END_REF].

Let (x d 1 (t), x d 2 (t) = ẋd 1 ) represent a desired trajectory of the same structure as that of the unilaterally constrained plant, i.e., the desired position is a continuous function of time whereas the desired velocity is piecewise continuous and it undergoes restitution in the opposite direction whenever x d 1 (t) = 0. The following assumptions are made throughout.

1) The desired trajectory x d 1 (t), x d 2 (t) is a periodic function of time with period T d such that x 1 (t) = 0 ∀t = nT d and x 1 (nT d ) = 0, n = 0, 1, 2,. .., thereby causing the periodic velocity restitution x d 2 (nT d+ ) = -Ksign(x d 2 (nT d-)) with some positive constant K. 2) An impulse event occurs in infinitesimally small time.

The impact event is detected.

3) ess sup

(t 0 )|, |x 2 (t 0 )|, |x d 1 (t 0 )|, |x d 2 (t 0 )|},
where t 0 is the initial time, is known a priori and is finite so that R > K + M ∆ . The first assumption ensures that the chosen desired trajectory is periodic as the reset event in desired velocity x d 2 (t) causes the post reset velocity x d 2 (nT d+ ) to take the same value after every impact for all n = 0, 1, 2,. ... For example, setting

x d 1 (t) = 1 2 ᾱt 2 + βt + γ, where ᾱ = -2 K T , β = K, γ =
0 are scalars, gives a periodic solution in the state-space

(x d 1 , x d 2 )
. Then, upon using a finite time stable controller to achieve

x 1 = x d 1 , x 2 = x d 2 before the next impact at t k with settling time T tw < T d , the constraint surface becomes S = {(x 1 , x 2 ) : x 1 = x d 1 = 0}.
The second assumption is made simply to show the instantaneous nature of the impact event and to enable restitution in the desired velocity. The third assumption represents the upper bound on the uncertainty ω(t) and ω ∆ respectively. The last assumption dictates that the results presented in the paper are of a local nature 1 . The following tracking dynamics can be obtained by defining the error variables e

1 (t) = x 1 (t) -x d 1 (t), e 2 (t) = x 2 (t) -x d 2 (t): ė1 = e 2 ė2 = u(e) + ω(t) -ẋd 2 (t) (x 1 , x 2 ) ∈ S e 2 (t + k ) = -ēx 2 (t - k ) + ω ∆ -x d 2 (t + k ) ( x 1 , x 2 ) ∈ S (4) When ẋd 2 (t)
is bounded, the first line of (4) coincides with the standard perturbed double integrator [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF] and the twisting controller drives the errors e 1 , e 2 to the origin in arbitrarily small finite time. Hence, it is guaranteed that the discrete event t k and nT d coincide. Furthermore, due to chatter-free numerical integration, the coincidence of t k and nT d holds true when the uncertainty ω = 0 as in this case numerically integrated states match the continuous time counterpart identically even in sliding. However, t k does not coincide with nT d when ω is present since the discrete time states differ from the continuous time ones by the order dictated by the sampling time. When t k coincides with the impacting time of the desired trajectories i.e. t k = nT d , then the impact law satisfies e 2 (t

+ k ) = -ēx 2 (t - k )+ω ∆ -Ksignx d 2 (t - k ).
The aim of the paper is to propose (i) a switched state feedback synthesis, (ii) the corresponding tuning rules to establish the gains of the controller to give an a priori guarantee of finite time convergence of the states of system (4) to the desired trajectory and (iii) to propose a digital implementation of the proposed continuous time synthesis.

III. FINITE TIME TRACKING OF PERIODIC TRAJECTORY

The motivation to introduce the following definition is to accommodate the uncertainty at the time of impacts.

Definition 1: A unilaterally constrained system is said to possess finite time impact attenuation in the closed-loop if for every ε > 0, there exists an open neighborhood Ω R around the origin such that every solution φ (t) of ( 4) that satisfies φ (0) = p, also satisfies (i) φ (t) = 0 for all t ∈ [nT d +T tw , (n+ 1)T d ) and for some T tw < T d and (ii)

|φ (t)| < ε for all t ∈ [nT d , nT d + T tw ) for all p ∈ Ω R, where n = 0, 1, 2,. .., ∞.
Remark 1: The uncertainty ω ∆ at the time of jumps in the actual plant velocity makes the error states deviate from the desired periodic trajectory at every jump and hence finite time stability can not in general be proven for all t ≥ 0. This is the main motivation for the definition of finite time impact attenuation. It should be noted that the differentiability of the solution φ (t) is not required. As will be shown in this section, this is because it is possible to construct a switched synthesis for which an absolutely continuous solution φ (t) in the sense of Fillipov [START_REF] Filippov | Differential Equations with Discontinuous Righ thand Sides[END_REF] is shown to exist for all t ∈ [nT d , (n + 1)T d ) where n = 0, 1, 2,. .., ∞. Let the switched control synthesis be given as follows:

u(e 1 , e 2 ) = Le, (e 1 , e 2 ) / ∈ Ω R ; Φ(e), (e 1 , e 2 ) ∈ Ω R , (5)
where,

L = -l 1 -l 2 T , Φ(e) = -µ 1 sign(e 2 ) -µ 2 sign(e 1 ) Ω R = {(e 1 , e 2 ) : µ 2 |e 1 | + 1 2 e 2 2 ≤ R}. (6)
Let the following a priori tuning rules be utilised:

µ 1 > max N, 2δ T s 2 √ 1-η 2 + N µ 2 > max R 2 , µ 1 -N η , µ 1 + N, ρ R 2(1-ρ) , ρ, β µ 1 l 1 = α λ 2 , l 2 = (α + 1) λ (7) where λ ≥ 2µ 2 N+(α+1) 2 N 2 r , δ > r 1 (β +1) β -1 , if r 0 > r 1 (β +1) β -1 ; r 0 (β +1) β -1 , if r 0 < r 1 (β +1) β -1 , λ = max 2µ 2 N+(α+1) 2 N 2 r , κ , r 2 = ρR µ 2 µ 2 (α + 1)|e 1 0 | + 2µ 2 κ |e 2 0 | e -κ T s 1 + 2α 2 κ 2 e 2 1 0 + (α 2 + 1)e 2 2 0 e -κ T s 1 ≤ R 1 , (8)
T s 1 = ε 1 T d , T s 2 = (1 -ε 1 )T d if r 0 > δ ; T s 1 = 0, T s 2 = ε 1 T d otherwise; (9)
r 1 = √ 2R, r ∈ (0, r 2 ), R 1 = R -r, α > 1 is an arbitrary scalar, β > 1 is a tuning variable, η ∈ (0, 1 β ),
R is an arbitrary positive constant representing the switching boundary, ε 1 ∈ (0, 1), ρ ∈ (0, 1) are arbitrary scalars, r 0 = (e 1 0 ) 2 + (e 2 0 ) 2 is the upper bound on the Euclidian norm of the system initial conditions where e 1 0 = e 1 (t 0 ), e 2 (t 0 )+M ∆ = e 2 0 . It should be noted that there always exists a solution κ such that (9) is satisfied (and can be obtained using numerical optimization routines). The basis for this is that lim κ→∞ κ 2 e -κ t = 0 for all t ≥ 0. The quantity on the left hand side of the inequality (9) can be viewed as a function of κ. The parameter µ 2 is a fixed entity on the left hand side of (9) due to the order in which the parameters µ 1 , µ 2 , l 1 , l 2 are tuned. Furthermore, parameter µ 1 is a function of T s 2 and parameters l 1 , l 2 are functions of T s 1 .

In the following, a switched discontinuous controller is utilised and hence, the vector field of the closed-loop system is discontinuous. Existing results [START_REF] Ye | Stability analysis of systems with impulse effects[END_REF] are not applicable and a new proof is required. An a priori guarantee stemming from ( 7) is utilised for obtaining finite time convergence to the origin between the impacts. In this case, the stability of the system can at best be defined in the sense of Definition 1 due to the uncertainty at the time of the jumps. The following result is in order.

Theorem 1: Assume that the trajectory x d 1 (t), x d 2 (t) is designed such that assumptions (1), ( 2), ( 3) and (4) hold true. Then, the closed-loop system (4), (5) possesses finite time impact attenuation if the tuning rules [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF] are employed for determining the gains l 1 , l 2 , µ 1 , µ 2 .

Proof: Let the neighbourhood Ω R around the origin be chosen a follows:

Ω R = {max{|x 1 (t 0 )|, |x 2 (t 0 )|, |x d 1 (t 0 )|, |x d 2 (t 0 )|} ≤ R} (10)
The tuning rules [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF] ensure that the error variables (e 1 (t), e 2 (t)) ∈ R 2 settle to the origin (0, 0) of the continuous dynamics in less time than T d for all (x 1 (t 0 ), x 2 (t 0 ), x d 1 (t 0 ), x d 2 (t 0 )) ∈ Ω R when the parameters e 10 , e 20 are chosen greater than R in [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF]. See [START_REF] Oza | Lyapunov based settling time estimate and tuning for twisting controller[END_REF] for detailed mathematical proof. What remains to be proven is the existence of finite time impact attenuation.

The solutions of the closed-loop system are given in [START_REF] Filippov | Differential Equations with Discontinuous Righ thand Sides[END_REF]. The upper bounds N, M ∆ on disturbances ω(t), ω ∆ are sufficient to compute the largest time required by the trajectory to enter the level set Ω R of the Lyapunov function V (e 1 , e 2 ) = µ 2 |e 1 |+ 1 2 e 2 2 in finite time when a linear controller is being used. Although the control u, when switched to the second order sliding mode controller, is discontinuous on the axes e 1 (t) = 0, e 2 (t) = 0, the continuous flow can be integrated piece-wise in each of the four quadrants

G 1 = {e : e 1 > 0, e 2 > 0} ,G 2 = {e : e 1 > 0, e 2 < 0} , G 3 = {e : e 1 < 0, e 2 < 0} ,G 4 = {e : e 1 < 0, e 2 > 0} .
(11) The integrals [START_REF] Filippov | Differential Equations with Discontinuous Righ thand Sides[END_REF] prove uniform boundedness of the states as a function of time in the neighbourhood B r 0 where B r 0 is a ball centered at the origin with a radius r 0 > 0. See [START_REF] Oza | Lyapunov based settling time estimate and tuning for twisting controller[END_REF] for detailed computation of the integrals in [START_REF] Filippov | Differential Equations with Discontinuous Righ thand Sides[END_REF]. Let the trajectory be analysed in the period t ∈ [nT d , (n + 1)T d ) for any n. The scalar e 0 = r 0 = e 2 10 + e 2 20 is the upper bound on the norm of the periodic initial conditions (or equivalently jumps) e(nT d ) because R > K + M ∆ .

By definition of the desired trajectory x d 2 , the time T d between the two jumps is known. The switched controller u(e 1 , e 2 ) is tuned as per [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF] to obtain an arbitrarily small settling time T tw < T d for an arbitrary initial condition. Since K < R, all post-impact errors e(nT d+ ) settle at the origin in less time than T tw due to the fact that the upper bound M ∆ on ω ∆ is embedded in the tuning rules in terms of worst case value of e 20 . Hence, it holds that for any arbitrary initial condition e(nT d+ ) ∈ Ω R where n = 0, 1, 2,. .., ∞, the continuous flow (12) of the system (4), (5) between the impacts settles at the point e 1 = 0, e 2 = 0 before the next jump event t = (n + 1)T d-. Then it follows that the jumps in both the actual and the desired velocities occur at the same time t = t k = nT d . Hence, the trajectory of (4), (5) always passes through the point (e 1 , e 2 ,t) = (0, 0, nT d ), where n = 0, 1, 2,. .., ∞.

In the presence of uncertainty ω ∆ , the system (4),(5) possesses finite time impact attenuation in the sense of Definition 1 as there exist finite scalars ε > 0, δ > 0, T tw < T d such that for all (e 1 (nT d+ ), e 2 (nT d+ )) <

e 1 (t) e 2 (t) =                    αe 1 0 + 1 λ e 2 0 e -λ t - e 1 0 + 1 λ e 2 0 e -α λ t + t 0 1 λ e -λ (t-τ) -e -αλ (t-τ) ω(τ)dτ - αλ e 1 0 + e 2 0 e -λ t + α λ e 1 0 + e 2 0 e -α λ t + t 0 -e -λ (t-τ) + αe -αλ (t-τ) ω(τ)dτ     , if e ∈ Ω R ;
e At e(0) + t 0 e A (t-τ) Bu(τ)dτ, otherwise. [START_REF] Filippov | Differential Equations with Discontinuous Righ thand Sides[END_REF] where,

e = e 1 e 2 T , A = 0 1 0 0 , B = 0 1 , u(t) = -(µ 1 -N)sign(e 2 ) -µ 2 sign(e 1 ) (13)
ε and (e 1 (nT d+ ), e 2 (nT d+ )) ∈ Ω R, the inequalities (i)

(e 1 (T tw ), e 2 (T tw ) = (0, 0) for all t ∈ [nT d +T tw , (n+1)T d ) and (ii) e(t) < δ , δ = ε + |M ∆ | hold true ∀t ∈ [nT d , nT d + T tw ].
Remark 2: The above theorem proves Lyapunov stability (and not asymptotic stability) of the proposed tracking problem. However, it proves finite time stability in the presence of impacts ∀t ≥ 0 in the case of gains µ 1 , µ 2 , l 1 , l 2 approaching infinity as this gives T tw → 0 (see [START_REF] Oza | Lyapunov based settling time estimate and tuning for twisting controller[END_REF] for theoretical proof).

IV. CHATTERING-FREE DIGITAL IMPLEMENTATION

In this section, a numerical chattering-free implementation of the second-order sliding mode control is presented. This implementation of sliding mode control has been developed previously on non-impacting systems [START_REF] Acary | Implicit euler numerical scheme and chattering-free implementation of sliding mode systems[END_REF], [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF] and is mainly based on an implicit (backward) Euler discretization of the unperturbed state dynamics with impact

ẋ1 (t) = x 2 (t), ẋ2 (t) = u(t), (x 1 , x 2 ) ∈ S x 2 (t + k ) = -ēx 2 (t - k ), (x 1 , x 2 ) ∈ S (14)
together with the twisting controller [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF]

u(t) ∈ -µ 2 sign(x 1 (t) -x d 1 (t)) -µ 1 sign(x 2 (t) -x d 2 (t)). (15)
The notation in terms of the inclusion in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF] outlines the fact that the function sign(.) is considered as a multi-valued mapping

sign(x) =    1, x > 0, -1, x < 0, [ -1, 1], x = 0. (16)
The time discretization of the impacting state dynamics follows the lines of Moreau's time-stepping scheme [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]. As usual in complementarity theory [START_REF] Cottle | The Linear Complementarity Problem[END_REF], the symbol x ⊥ y means that x T y = 0. Moreau's time-stepping scheme for the nonsmooth dynamics with the impact law (3) reads as

       x 1,k+1 = x 1,k + hx 2,k+1 , x 2,k+1 = x 2,k + hu k+1 + λ k+1 , if x 1,k + hx 2,k ≤ 0, then 0 ≤ x 2,k+1 + ēx 2,k ⊥ λ k+1 ≥ 0, otherwise λ k+1 = 0. ( 17
)
where h = t k+1t k is the sampling time. As in the previous work on implicit implementation of sliding mode controllers [START_REF] Acary | Implicit euler numerical scheme and chattering-free implementation of sliding mode systems[END_REF], [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF], the computation of the control input u k+1 is given by an implicit evaluation of the control law [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF],

u k+1 ∈ -µ 2 sign(x 1,k+1 -x d 1 (t k+1 )) -µ 1 sign(x 2,k+1 -x d 2 (t k+1 )). (18)
Contrary to event-driven schemes, the time discretization ( 17), ( 18) is solved in a single step without resorting to an event-location procedure. When an impact is not forecast (x 1,k + hx 2,k > 0), the value of the impulse λ k+1 vanishes and it remains to solve the generalized equation

   W x k+1 = x k + hBs k+1 , e k+1 = x k+1 -x d (t k+1 ) s k+1 ∈ sign(e k+1 ) (19)
with

W = 1 -h 0 1 , B = 0 0 -µ 2 -µ 1 . (20)
Rearranging (19) produces, e k+1 = q + Ms k+1 s k+1 ∈ sign(e k+1 )

where q = W -1 x kx d (t k+1 ) and M = hW -1 B. The problem (21) appears as a Mixed Complementarity Problem (MCP), a well-known problem in the Mathematical Programming community [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF] for which numerous solvers are available. When an impact is forecast (x 1,k + hx 2,k ≤ 0), the problem is slightly more complicated as it is necessary to solve

           e k+1 y k+1 = p + N s k+1 λ k+1 s k+1 ∈ sign(e k+1 ) 0 ≤ y k+1 ⊥ λ k+1 ≥ 0 (22) with p = q H T (q + x d (t k+1 )) + ēx 2,k , H = 0 1 , (23)
and

N = M W -1 H H T M H T W -1 H . ( 24
)
Problem ( 22) is an MCP that can be easily solved. Once one of the problems (21) or ( 22) is solved, the control is computed according to

u k+1 = Bs k+1 . (25)
V. NUMERICAL SIMULATIONS In this section, numerical simulations are presented using the digital implementation developed above. Given a constant control input u k+1 for the time-step [t k ,t k+1 ], the simulations are performed by solving the dynamics with disturbances on [t k ,t k+1 ] as follows:

ẋ1 (t) = x 2 (t) + ω(t), ẋ2 (t) = u k+1 , (x 1 , x 2 ) ∈ S x 2 (t + k ) = -ēx 2 (t - k ) + ω ∆ , (x 1 , x 2 ) ∈ S (26)
For the sake of simplicity, Moreau's time-stepping ( 17) is used with the finer time-step h ≪ h to mimic the behaviour of the continuous time plant. This integration could also be performed with an accurate event-driven procedure using, for example, a zero-order hold integration before and after the impact. The numerical simulations are obtained with the SICONOS open-source package 2 dedicated to nonsmooth dynamical systems [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF]. For all the simulations presented below, the initial conditions have been chosen as (x d 1 (0), x d 2 (0), x 1 (0), x 2 (0) = (0, 4, 0, 3). The desired trajectories x d (t) have been designed such that

x d 1 (t) = x d 2 (0)(1 - t T d )t x d 2 (t) = ẋd 1 (t), t ∈ [0, T d ), (27)
The tuning rules [START_REF] Brogliato | Chattering-free digital slidingmode control with state observer and disturbance rejection[END_REF] can then be employed to compute the controller gains as follows: µ 1 = 3.787, µ 2 = 9.165. The value of K is chosen such that K = x d 2 (0) and the desired period is chosen such that T d is greater than 10 (T d = 12 in the presented simulations). The tuning rules (7) are then applied with T d = 10.

Figure 1 shows the finite settling time behavior of the closed-loop system with the proposed implicit digital implementation without disturbances for the sampling rate h = 0.1. It should be noted that a large error occurs after each impact due to the inelastic impact ( ē = 0.5). Nevertheless, the control algorithm is able to track the desired trajectory in less than ten seconds and without numerical chattering. The next simulation is performed with the introduction of disturbances in both the continuous time-dynamics and discrete event dynamics. The disturbances are modeled:

ω(t) = M sin(ϕt), ω ∆ = M ∆ random(t), (29)
with the frequency ϕ = 0.3 and the upper-bound M = 1, M ∆ = 0.4. The symbol random(t) denotes a function that generates random values in [-1, 1]. Figure 2 shows the behavior for the same sampling rate h = 0.1. Between impacts, the control counteracts the continuous time disturbances ω(t). To complete this section on the numerical simulations, Fig. 3 shows the behavior of the system for h = 0.1 without disturbances and with the control defined using a standard (explicit) digital implementation of the control as follows:

u k = -µ 2 sign(x 1,k -x d 1 (t k )) -µ 1 sign(x 2,k -x d 2 (t k )). ( 30
)
It should be noted that the numerical chattering due to the explicit implementation of the controller is not observed in the implicit implementation. It is noted, however, that increasing the sampling frequency will reduce the chattering magnitude and increase the frequency of chattering when an explicit integration scheme is used for the presented problem, numerical simulation of which is not presented due to want of space.

VI. CONCLUSIONS AND FUTURE WORK

Tuning rules are established when a second order sliding mode controller is utilized to achieve finite time tracking for a class of unilaterally constrained planar systems in the presence of external disturbances in continuous and discreteevent phases. Numerical chatter-free digital implementation is also presented. These results, while being based on strong theoretical developments, appear to make the latest advances in the area closer to industrial application. The results presented also have the potential for extension. Stability of periodic orbits can be studied via a method that computes the Poincaré maps and their Jacobian [1, Ch. 7,8]. This method requires explicit solutions of the closed-loop systems. The presented result is not a valid candidate for the aforementioned Poincaré based method as it relies on a mathematical proof involving geometric series [START_REF] Oza | Lyapunov based settling time estimate and tuning for twisting controller[END_REF]. It is then an interesting open problem to render the presented method a valid candidate for the aforementioned Poincaré impact map results and to develop a new methodology to compute the Jacobian of Poincaré maps when a finite time convergent SOSM synthesis is used.

  t≥0 |ω| ≤ N and ess sup k≥0 |ω ∆ | ≤ M ∆ where N, M ∆ are a priori known positive scalars.

  and x d (t + nT d ) = x d (t) for t ∈ [0, T d ), n ∈ N. The coefficient of restitution ē is set to 0.5. The aforementioned tuning steps have been carried out rendering the tuning variables as follows: ρ = 0.99, r = 0.16, R = 1.5, ε 1 = 0.3, η = 0.45, β = 2.2, α = 1.01.

2Fig. 1 .

 1 Fig. 1. Finite settling time behavior of the system (14),(15) between the impacts with the implicit digital implementation. a) position x 1 . b) desired position x d 1 . c) velocity x 2 . d) desired velocity x d 2 . u) control u. Sampling time h = 0.1.

Fig. 2 .

 2 Fig. 2. Finite settling time behavior of the system (14),(15) between the impacts with the implicit digital implementation and disturbances. a) position x 1 . b) desired position x d 1 . c) velocity x 2 . d) desired velocity x d 2 . u) control u. Sampling time h = 0.1.

Such an upper bound is generally known a priori for a large class of mechanical systems.

This work was partially supported by EPSRC via research grant EP/G053979/1 and by CONACYT via grant 165958. Yury V.


Benoît Escrig
email: escrig@enseeiht.fr

Random Matrix Theory applied to the Estimation of Collision Multiplicities

Keywords: multi-packet reception, collision multiplicity, model order selection, IEEE 802.11-based networks

. This estimation step is a key task in multi-packet reception approaches and in collision resolution techniques. The two techniques are proposed for IEEE 802.11 networks but they can be used in any OFDM-based system. The techniques are based on recent advances in random matrix theory and rely on eigenvalue statistics. Provided that the eigenvalues of the covariance matrix of the observations are above a given threshold, signal eigenvalues can be separated from noise eigenvalues since their respective probability density functions are converging toward two different laws: a Gaussian law for the signal eigenvalues and a Tracy-Widom law for the noise eigenvalues. The first technique has been designed for the white noise case, and the second technique has been designed for the colored noise case. The proposed techniques outperform current estimation techniques in terms of mean square error. Moreover, this paper reveals that, contrary to what is generally assumed in current multi-packet reception techniques, a single observation of the colliding signals is far from being sufficient to perform a reliable estimation of the collision multiplicities.

I. INTRODUCTION

Collisions are known to degrade the throughput of random access wireless networks, such as ad hoc networks 1 . A collision occurs when two or more user nodes access the channel simultaneously. Over the last decades, Medium Access Control (MAC) protocols have been designed with the rationale that all data packets from the colliding user nodes are lost when a collision occurs because the signals from all users mix.

It is is possible to retrieve part of the data packets that are involved in a collision, with approaches such as the capture effect [START_REF] Lang Tong | Signal processing in random access[END_REF]- [START_REF] Durvy | Self-Organization Properties of CSMA/CA Systems and Their Consequences on Fairness[END_REF]. Sophisticated techniques based on Multi-User Detection (MUD) allow the decoding of more than one data packets. MUD receivers have been successfully implemented in a wide range of application areas [START_REF] Verdú | Multiuser Detection[END_REF]. In this context, the number of colliding user is often needed to efficiently parameterize the MUD receivers.

Another approach is often invoked when a collision occurs. It consists in triggering a collision resolution (CR) mechanism where transmissions from the colliding users are re-scheduled in order to avoid another collision [START_REF] Lang Tong | Signal processing in random access[END_REF], [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF]- [START_REF] Dimic | Medium access control -physical cross-layer design[END_REF]. In this context, CR mechanisms operate more efficiently when the number of colliding user is known. For instance, one could increase the contention window of a CR protocol with respect to the estimate of the number of colliding nodes. Note that our intention is not to propose a new CR mechanism but rather to provide a new parameter to improve the tuning of the mechanisms.

So, the purpose of this study is to estimate the number of colliding users, i.e., the collision multiplicity (CM) [START_REF] Escrig | Estimation of Collision Multiplicities in IEEE 802.11based WLANs[END_REF]. We focus here on the following scenario. The receiver at the destination node is processing a collision signal that consists of a mixture of signals from all users plus some Additive White Gaussian Noise (AWGN). The destination node is the node toward which all users involved in the collision are intending to send data 2 . So, from the observations of mixtures, the destination node has to estimate the number of original signals. This problem is a typical Model Order Selection (MOS) problem. MOS problems arise in the signal processing area and related areas such as signal array processing, radar, and sonar processing.

The MOS techniques are all based on the following rationale: the mixture of signals and noise can be decomposed into a noise subspace and a signal subspace, and the dimension of the signal subspace equals the number of signals.

In order to perform this separation, the following steps are implemented. First, T observations of the mixture are gathered by the processing node. Then, the sample covariance matrix (SCM) of these observations is computed and an eigenvalue decomposition of this matrix is performed. These observations are obtained over T time-slots. The MOS techniques use the property that signal eigenvalues are much higher than noise eigenvalues, provided that the Signal to Noise Ratio (SNR) is high enough. When T eigenvalues are available and when K eigenvalues are significantly higher than the T -K remaining eigenvalues, the conclusion is that there are K signals in the mixture.

The proposed approach has also been motivated by the following two observations: (i) many CM estimation techniques are based on the assumption that signal samples are uniformly distributed over a finite alphabet, i.e., signal samples are modulation symbols such as PSK or QAM symbols [START_REF] Zheng | Multipacket Reception in Wireless Local Area Networks[END_REF]- [START_REF] Talwar | Blind Separation of Synchronous Co-Channel Digital Signals Using an Antenna Array. Part I. Algorithms[END_REF], and (ii) in many existing techniques, the number of observations is not much greater that the number of signals [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF], [START_REF] Özgül | Wireless Access with Blind Collision-Multiplicity Detection and Retransmission Diversity for Quasi-Static Channels[END_REF].

Our point concerning (i) is to design estimation techniques that could be used in the context of Gaussian distributed samples. In this paper, we focus on the wireless standards that implement Orthogonal Frequency Division Multiplex (ODFM) transmissions; so signal samples can be considered as being Gaussian distributed and not uniformly distributed. 3 The issue that derives from (ii) is to clarify whether the number of observations can be made as small as K + 1 or K + 2 as in [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF], [START_REF] Özgül | Wireless Access with Blind Collision-Multiplicity Detection and Retransmission Diversity for Quasi-Static Channels[END_REF] or not. Indeed, assuming that T and K are on the same order is in stark contrast with the typical assumptions that are used in signal processing [START_REF] Stoica | Model-Order Selection : A review of information criterion rules[END_REF]- [START_REF] Wax | Detection of Signals by Information Theoretic Criteria[END_REF].

We shall show that the proposed techniques outperform current techniques in terms of Mean Square Error (MSE) and that the number of observations T should be much higher than the number of signals K in order to obtain satisfactory MSE performance.

The rest of the paper is organized as follows. The introduction end ups with a Related Work subsection. The system model is introduced in Section II and some results on eigenvalue statistics are stated in Section III. The CMETs are described in Section IV. Simulation results are presented in Section V and a conclusion is drawn in the last section.

Related work

The first contributions, in the field of eigenvalue-based MOS techniques, were developed by Bartlett [18] and Lawley [START_REF] Lawley | Tests of significance of the latent roots of the covariance and correlation matrices[END_REF]. They propose a subjective setting of the threshold that is used to separate signal eigenvalues from noise eigenvalues. This approach is still used in some CR mechanisms [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF], [START_REF] Özgül | Wireless Access with Blind Collision-Multiplicity Detection and Retransmission Diversity for Quasi-Static Channels[END_REF] in order to minimize the number of observations T . The algorithm starts with T = 1 and the number of observations is incremented by one each time the smallest eigenvalue is higher than a noise threshold. As soon as the smallest eigenvalue crosses the noise threshold at step T , that means that the T th highest eigenvalue is a noise eigenvalue, and hence that there are T -1 signals. Information theoretic criteria such as the Akaike Information Criterion (AIC) and the Minimum Description Length (MDL), developed by Wax and Kailath [START_REF] Wax | Detection of Signals by Information Theoretic Criteria[END_REF], have then been proposed in order to alleviate the limiting constraint imposed by the subjective setting of the separation threshold. The criteria are usually composed of a function that depends on the maximum likelihood estimator of the parameters of the model and a term that adjusts the first component to the context. This second term depends on the parameters of the system such as the number of samples per observation. The MDL have been used over several decades in several areas. An interesting review of this criterion can be found in [START_REF] Nadakuditi | Sample Eigenvalue Based Detection of High-Dimensional Signals in White Noise Using Relatively Few Samples[END_REF].

The MDL is a consistent estimator of the number of signals when the number of observations, T , is fixed and when the number of samples per observations, denoted m, is such that m → +∞, provided that T is much larger than the number of signals K. Two limitations have recently been pointed out for this criterion.

First, the MDL have been shown to be inconsistent as the variance of the noise tends toward zero [START_REF] Ding | Inconsistency of the MDL: On the Performance of Model Order Selection Criteria with Increasing Signal-to-Noise Ratio[END_REF]. Second, the MDL is based on the distribution of the signal eigenvalues. It uses the property that sample eigenvalues, i.e., the eigenvalues that are obtained from the eigenvalue decomposition, are symmetrically distributed around the population eigenvalues, i.e., the theoretical eigenvalue. Basically, all the above-mentioned techniques are based on this property about signal eigenvalues. This centrality assumption makes sense in the T fixed, m → +∞ case. However, when T and m are on the same order, even when T, m → +∞, the previous property is no longer valid.

Recent advances in the Random Matrix Theory (RMT) field have brought into light several properties on the distribution of both signal and noise eigenvalues [START_REF] Johnstone | High Dimensional Statistical Inference and Random Matrices[END_REF], [START_REF]Multivariate analysis and Jacobi ensembles: Largest eigenvalue[END_REF]. These properties have been used in order to design new MOS techniques in the context where T, m → +∞ [START_REF] Nadakuditi | Sample Eigenvalue Based Detection of High-Dimensional Signals in White Noise Using Relatively Few Samples[END_REF], [START_REF] Nadakuditi | Fundamental Limit of Sample Generalized Eigenvalue Based Detection of Signals in Noise Using Relatively Few Signal-Bearing and Noise-Only Samples[END_REF]. These new techniques have been shown to outperform the classical MDL estimator.

Our purpose in this paper is to use these new properties in order to design new CM estimation techniques. Note that we are not in the context where T, m → +∞ but rather in the T fixed, m → +∞ case. The point is that current estimation techniques, whether they are based on the RMT or not, are based on the assumption that the number of observations T is much larger than the number of signals K. In our context, we want to minimize T with respect to K, so we are dealing with a context where T is in the order of K. Since the new RMT-based MOS techniques are performing in a wider range of parameter values, we believe that they can be considered as relevant candidates for our objective.

II. SYSTEM MODEL

In this paper, we take in interest collision resolution algorithms and, more precisely, the estimation of the number of stations that are involved in a collision. We consider the scenario where K stations are simultaneously transmitting data packets to the same destination node (see Fig. 1).

We assume that the destination code and the colliding stations are all equipped with a single antenna. This is a worst case scenario. When the destination node is equipped with several antennas, observations are gathered more rapidly. The K colliding stations are transmitting OFDM frames of m samples each. The OFDM signal samples are Gaussian distributed with zero mean and variance unity. Moreover, the OFDM signals from the colliding users are assumed to be uncorrelated. This assumption makes sense since each transmission is affected by both a different Doppler shift and a different Doppler spread. The destination node receives T observations. So it is assumed that the destination can trigger transmissions from the colliding nodes. Note that T must be larger than K to allow an efficient estimation of K. This will be discussed in the simulation section. The destination triggers transmissions from the colliding users by sending a feedback frame. This frame serves also as a synchronization frame for user transmissions. Hence, the K users can be assumed to be coarsely synchronized in time. So this system model is similar to a source separation problem when K signals are impinging on a T sensor array. The T × 1 observation vectors y i can be written as

y i = Ax i + b i , i = 0, 1, . . . , m (1)
where x i ∼ CN K (0, R x ) are K ×1 complex Gaussian vectors with zero mean and covariance matrix R x , b i ∼ CN T (0, R b ) are T × 1 complex noise vectors that are Gaussian distributed with zero mean and noise covariance matrix R b . In the case of a white noise, we have that R b = σ 2 I T where σ 2 denotes the noise variance and I T is the T ×T identity matrix. The channel matrix A is considered as an unknown T × K non-random matrix. For our study, we assume that the coefficients of A are modeled as circularly symmetric Gaussian coefficients with power unity (Rayleigh fading). The channels gains are assumed to have constant values over the duration of the frame and change randomly from one frame to another. This corresponds to the typical block-fading channel assumption.

The observations y i can be whitened by the following transformation

y † i = R -1/2 b y i provided that the noise covariance matrix R b is known a priori and is nonsingular, where R +1/2 b
is the Hermitian nonnegative definite square root of R b . The transformation simply reduces to a normalization step in the case of a white Gaussian noise. The covariance matrix R of the snapshots y i is given by

R = AR x A H + R b = Ψ + R b
with H denoting the complex conjugate, the signal and noise vectors being independent. We assume that the channel matrix A is full rank and that the signal covariance matrix R x is nonsingular so that the rank of Ψ is min(K, T ). Hence, when T ≤ K, there are T non-zero eigenvalues in the matrix Ψ and when T > K, there are K non-zero eigenvalues. This property is used in [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF], [START_REF] Özgül | Wireless Access with Blind Collision-Multiplicity Detection and Retransmission Diversity for Quasi-Static Channels[END_REF] where the number of observations, T , is incremented until R has zero eigenvalues. When the whitening transformation is applied, the covariance matrix R † of the whitened observations, denoted R † , is defined as

R † = R -1/2 b RR -1/2 b = R -1/2 b ΨR -1/2 b + I T The population eigenvalues R † , denoted λ 1 ≥ λ 2 ≥ ... ≥ λ T , are such that λ i > 1 for 1 ≤ i ≤ min(K, T ) (2) λ i = 1 for min(K, T ) < i ≤ T (3)
In most approaches, T is much larger that K so the conditions in ( 2) and ( 3) should be written as 1 ≤ i ≤ K and K < i ≤ T respectively. However, in some CR algorithms such as the ones in [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF] and in [START_REF] Özgül | Wireless Access with Blind Collision-Multiplicity Detection and Retransmission Diversity for Quasi-Static Channels[END_REF], ( 2) and ( 3) are used in this way, since they are initialized with T = 1 and T is incremented by one each time the T th smallest eigenvalue is not detected has being a noise eigenvalue. In any case, the estimation of K can be performed from the multiplicity of the λ i equalling one. When R and R b are known, and when the rank of R -1 b Ψ is K, the CM estimation can be easily performed from the multiplicity of the λ i equalling one. Otherwise, when R and R b are unknown and have to be estimated, the sample covariance matrix (SCM) of the observations y i , denoted R, must be computed

R = 1 m m i=1 y i y H i
The SCM of the noise, denoted R b , must be also computed using

R b = 1 N N j=1 b j b H j
Note that the b j , 1 ≤ j ≤ N are independent noise-only samples. We assume that we can get noise-only samples by using idle time slots. In that case, the estimation of K is based on the eigenvalue decomposition of the SCM R

† R † = 1 m m i=1 y † i (y † i ) H = R -1/2 b R R -1/2 b (4)
The sample eigenvalues of the SCM matrix R † are denoted λ1 ≥ λ2 ≥ ... ≥ λT .

III. SOME RESULTS ON EIGENVALUE STATISTICS

In the following, we investigate new results in RMT. In particular, we provide new characterizations of the eigenvalues of SCMs. We first begin with the signal-free case and then address the signal bearing case.

A. Signal-free Case

We start with the signal-free case where no user is transmitting, so K = 0. As a consequence, we have that

R = R b R † = I T (5)
and all population eigenvalues λ i are all equal to one.

From the system model and (5), the vectors y i are Tdimensional complex Gaussian independent vectors such that

y i ∼ CN T (0, R b ).
From Property 6.1, we get that m R is a T -variate complex Wishart matrix with m degrees of freedom and covariance matrix R b , i.e.,

m R ∼ CW T (m, R b )
Similarly, since we have that y † i ∼ CN T (0, I T ), we have that m R † ∼ CW T (m, I T ) From Property 6.2, we now characterize the eigenvalues of R in the signal-free case and in the presence of a white Gaussian noise [START_REF] Nadakuditi | Sample Eigenvalue Based Detection of High-Dimensional Signals in White Noise Using Relatively Few Samples[END_REF], [START_REF] Johnstone | High Dimensional Statistical Inference and Random Matrices[END_REF], [START_REF] Perry | Minimax Rank Estimation for Subspace Tracking[END_REF], [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF].

Corollary 3.1: In the signal-free case, when y i ∼ CN T (0, λI T ), the Empirical Spectral Distribution (ESD) of the T × T random matrix R † converges almost surely to the Marcenko-Pastur law in ( 14) (see Fig. 2). When the snapshots y i are whitened and become y † i , a new distribution for the ESD must be considered [START_REF] Nadakuditi | Fundamental Limit of Sample Generalized Eigenvalue Based Detection of Signals in Noise Using Relatively Few Signal-Bearing and Noise-Only Samples[END_REF].

From Property 6.3, we now characterize the eigenvalues of R † in the signal-free case in the presence of an arbitrarily colored Gaussian noise [START_REF] Johnstone | High Dimensional Statistical Inference and Random Matrices[END_REF], [START_REF] Nadakuditi | Fundamental Limit of Sample Generalized Eigenvalue Based Detection of Signals in Noise Using Relatively Few Signal-Bearing and Noise-Only Samples[END_REF]- [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF].

Corollary 3.2: In the signal-free case, when y † i ∼ CN T (0, I T ), the ESD of the T × T random matrix R † converges almost surely to the modified Marcenko-Pastur law in [START_REF] Nadakuditi | Sample Eigenvalue Based Detection of High-Dimensional Signals in White Noise Using Relatively Few Samples[END_REF] (see Fig. 3). In the signal-free case, we are interested in the eigenvalues of R † or, equivalently, in the eigenvalues θ that satisfy

Rv = θ R b v (6)
We can rewrite (6) as follows

m Rv = ( m N θ)N R b v (7) Since m R ∼ CW T (m, R b ) and N R b ∼ CW T (N, R b ), (7)
is similar to [START_REF] Johnstone | High Dimensional Statistical Inference and Random Matrices[END_REF]. So, using Property 6.5, the largest eigenvalue λ1 that satisfies [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF] is Tracay-Widom (TW) distributed and we can derive the following property. Property 3.1: The largest eigenvalue λ1 that satisfies (6) is TW distributed, i.e.,

P log( m N λ1 ) -µ(T, m, N ) σ(T, m, N ) ≤ x → T W C (x)
The pdf of the largest eigenvalue λ1 is represented in Fig. 4. Note that this characterization uses explicitly the double Wishart setting that has been motivated by the need to whiten the observations when the additive Gaussian noise of the channel is not white. Another and simpler characterization of λ1 has also be proposed in [START_REF] Johnstone | High Dimensional Statistical Inference and Random Matrices[END_REF], [START_REF]Multivariate analysis and Jacobi ensembles: Largest eigenvalue[END_REF].

Property 3.2: In the signal-free case, the whitened snapshots y † i are N T (0, I T ) and the largest eigenvalue

λ1 of the SCM R † is Tracy-Widom distributed. When T, m → ∞ such that T /m → c ∈ (0, ∞), P m λ1 -µ T,m σ T,m ≤ x → T W C (x)
where When the snapshots y † i are N T (0, σ 2 I T ), the convergence limit of m λ1 is σ 2 ( √ T + √ m) 2 . This corresponds to the non-normalized case. Note that the convergence rate to the

µ T,m = ( √ T + √ m) 2 σ T,m = ( √ T + √ m)( 1 √ T + 1 √ m ) 1/3
T W C (x) distribution function is O(T -1/3
). When the parameters m and T are not so large, which is practically the case when we want to reduce the number of observations, the convergence rate to the TW distribution is rather O(T -2/3 ) provided that the mean and standard deviation have been modified appropriately .

B. Signal Bearing Case

When there are K signals and when T → ∞, the limiting ESD of R † still converges to a Marcenko-Pastur distribution. This is because all eigenvalues are equally weighted by the ESD so the impact of K signals vanishes when T → +∞. Now there are K signal eigenvalues and T -K noise eigenvalues, and their respective characterization are different. 4 Noise eigenvalues are still distributed according to a TW distribution. The situation is a little bit more complicated for the signal eigenvalues since the characterization depends on a threshold τ . In the case of a white Gaussian noise, when no whitening transfomation in applied, the threshold is defined as τ = 1+

√ c [START_REF] Nadakuditi | Sample Eigenvalue Based Detection of High-Dimensional Signals in White Noise Using Relatively Few Samples[END_REF]. The threshold takes into account the c 1 ratio in Property 6.3, when dealing with a colored noise [START_REF] Nadakuditi | Fundamental Limit of Sample Generalized Eigenvalue Based Detection of Signals in Noise Using Relatively Few Signal-Bearing and Noise-Only Samples[END_REF]. When the i th largest signal eigenvalue λi is strictly higher than τ , the signal eigenvalue converges to a limit different from that in the signal-free case [START_REF] Nadakuditi | Sample Eigenvalue Based Detection of High-Dimensional Signals in White Noise Using Relatively Few Samples[END_REF]. In that case, the signal eigenvalue is Gaussian distributed, i.e., P λiµ i 4 This case is often referred to as a "spiked" model in RMT.

σ i ≤ x → G(x)
where

µ i = λ i (1 + c λ i -1 ) σ i = λ i √ T 1 - c (λ i -1) 2 and G(x) = x -∞ 1 √ 2π exp (- y 2
2 )dy

In Fig. 5, the pdf of the smallest signal eigenvalue is represented. Otherwise, when the i th largest signal eigenvalue λi is strictly lower that τ , the signal eigenvalue is distributed according to a TW distribution, i.e.,

P log( m N λi ) -µ(T -i, m, N ) σ(T -i, m, N ) ≤ x → T W C (x)
Note that there is a specific phenomenon when the i th largest signal eigenvalue exactly equals to the threshold [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF]. Hence, when K ≪ T , the signal eigenvalues exhibit, on rescaling, fluctuations described by the TW distributions when there are strictly below the threshold τ . So the distributions obtained for the signal-free case (K = 0) closely approximate the distribution of the signal eigenvalues. These results have an impact on the design of CM estimation techniques since signal eigenvalues strictly below τ are considered as being noise eigenvalues. Note that adding observations, i.e., increasing T , cannot be an option to tackle this problem since τ grows with T . However, when the signal eigenvalues are above the threshold, a reliable estimation of the CM is possible.

IV. COLLISION MULTIPLICITY ESTIMATION TECHNIQUES

We first present the CM estimation method based on the distribution of noise eigenvalues. Then we review two approaches based on the distribution of signal eigenvalues.

A. A method based on the distribution of noise eigenvalues

We assume that T observations are available at the destination node. The SCM R † is computed from (4). Then, an eigenvalue decomposition of R † is performed and the resulting eigenvalues are sorted in descending order λ1 ≥ λ2 ≥ ... ≥ λT . The method relies on the property that the i th largest noise eigenvalue is TW distributed, provided that all signal eigenvalues are located above the detectability threshold. So for every l = 0, 1, . . . , T -1, we test the null hypothesis that "there are exactly l colliding signals" against the alternative hypothesis that there are "at least l + 1 colliding signals". The test is performed by computing the test statistic S λ (l)

S λ (l) = log( m N λl+1 ) -µ(T -l, m, N ) σ(T -l, m, N ) (8)
and comparing it to a threshold τ α . The threshold τ α is defined as

τ α = T W -1 C (1 -α) ( 9
)
where α is some significance level. More precisely, a significance level α is set so that the probability that the null hypothesis is detected by chance is α. So the null hypothesis is accepted when we have

S λ (l) < τ α (10)
The sequence of tests begins with l = 0. If the test S λ (0) is passed, i.e., when S λ (0) < τ α , there are no colliding users. Otherwise, there is at least one signal. The procedure proceeds to the next step with l = 1. The tests are performed subsequently until the test is passed. When the test is passed at step l = q, the eigenvalue λq+1 is detected as being TW distributed. So λq+1 is a noise eigenvalue, and hence there are exactly q colliding signals. Once the number of colliding users has been determined, the destination node stops the CM estimation process and proceeds to the next processing block. The procedure is summarized in Algorithm 1. The

Algorithm 1 TWIT algorithm Compute R † Perform the eigenvalue decomposition of R † Sort the eigenvalues λi , i = 1, . . . , T of R † KTWIT ← 0 and T est ← F alse while T est = F alse and KTWIT < T do µ ← µ(T -KTWIT , m, N ) σ ← σ(T -KTWIT , m, N ) T est ← {σ -1 [log(m λ KTWIT+1 /N ) -µ] < τ α } if T est = F alse then KTWIT ← KTWIT + 1
else break end if end while mean µ(x, y, z) and the standard deviation σ(x, y, z) in the algorithm are defined in Property 6.4. Note that this criterion has been originally designed for arbitrary (or colored) noise [START_REF] Nadakuditi | Fundamental Limit of Sample Generalized Eigenvalue Based Detection of Signals in Noise Using Relatively Few Signal-Bearing and Noise-Only Samples[END_REF]. That is the reason why the algorithm uses the eigenvalues of R † .

B. Methods based on the distribution of signal eigenvalues

We present two CMETs that are all based on the Gaussian distribution of signal eigenvalues. We first present the wellknown MDL criterion, and then address another criterion based on recent advances in RMT.

1) The MDL criterion: The MDL criterion has been defined in [START_REF] Wax | Detection of Signals by Information Theoretic Criteria[END_REF]. This criterion has been used for decades in the area of signal array processing, and other related domains. The MDL estimator KMDL for the CM K is defined as

KMDL = argmin k=1,...,T {MDL(k)} where MDL(k) = -m(T -k) log g(k) a(k) + 1 2 k(2T -k) log(m)
where

g(k) = T i=k+1 λ 1 T -k i a(k) = 1 T -k T i=k+1 λi
where the λi denote the eigenvalues of R † with 1 ≤ i ≤ T , ordered in descending order. This estimator is consistent in the m → ∞ sense. One of the reason why the MDL criterion has been widely used over the past two decades comes from its robustness to model mismatch, in particular when the underlying assumptions of snapshots and noise Gaussianity can be relaxed [START_REF] Xu | Detection of number of sources via exploitation of centro-symmetry property[END_REF]- [START_REF] Fishler | Estimation of the number of sources in unbalanced arrays via information theoretic criteria[END_REF].

2) The criterion based on recent results on signal eigenvalues: This new criterion has been designed using the last results in RMT [START_REF] Nadakuditi | Sample Eigenvalue Based Detection of High-Dimensional Signals in White Noise Using Relatively Few Samples[END_REF]. We denote this estimator SEMOS (Signal Eigenvalue-based Model Order Selection). The estimator KSEMOS for the CM K is defined as, for complex data

KSEMOS = argmin k=0,...,T -1 1 2 ( 1 c ) 2 q 2 k + 2(k + 1)
The test statistic q k is defined as

q k = (T -k) g s (k) a s (k) -(1 + c) × T (11)
where

g s (k) = T i=k+1 λ2 i a s (k) = ( T i=k+1 λi ) 2 V. SIMULATION RESULTS
We study the performance of three estimation techniques: the MDL criterion, the TWIT, and the SEMOS estimator. The methods are evaluated in the context of Rayleigh blockfading channels. The channel coefficients A ij in [START_REF] Escrig | Estimation of Collision Multiplicities in IEEE 802.11based WLANs[END_REF] IEEE 802.11 standard [START_REF]IEEE 802.11n standard, Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) specifications: Amendment 5: Enhancements for Higher Throughput[END_REF] and use the signal parameters listed in Table I.

From the signal parameters, we get that the number of samples per OFDM frame is m = (1 + GI) × N sub × N OFDM = 61440. For the sake of simplicity, we set that N = m, i.e., N = 61440. The performance of estimation techniques have been evaluated over 10,000 Monte Carlo trials (N Simu = 10000) using the MATLAB software.

The three estimation techniques are compared in Fig. 6 to 10 in the case of a white Gaussian noise. The results have been obtained for both a variable T and a variable SNR. The significance level α for the TWIT is set to 0.01 [START_REF] Nadakuditi | Fundamental Limit of Sample Generalized Eigenvalue Based Detection of Signals in Noise Using Relatively Few Signal-Bearing and Noise-Only Samples[END_REF]. Figures 6 and7 represent the simulation results for K = 3 and K = 4 respectively. The simulation results are similar with K = 3 and K = 4 so we focus on the latter case in the following. When the number of observations and the SNR are increasing, the curves are converging toward the same position so it is difficult to compare them. So we choose instead to compare the estimation techniques with respect to the MSE ǫ 2 between the estimates K and K

ǫ 2 = 1 N Simu NSimu i=1 | Ki -K| 2 ( 12
)
where Ki is the estimate of K for the i th trial. This is shown in Fig. 8 to 10. Since the curves are plotted on a log scale, only non-zero values are plotted. So when the experimental MSE is zero, the corresponding point of the curve is not represented. The simulation results show that the estimation techniques perform better at high SNR and when T increases. The number of observations should be higher than almost three times the CM K in order to get a relative MSE lower than 10%. Similarly, the SNR should be higher than 0 dB in order to get a similar performance level, for any values of T . The two RMT-based techniques, i.e., the SEMOS estimator and the TWIT, outperforms the MDL criterion. However, the differences between the simulation curves decrease when the SNR or T is increasing. Moreover, the SEMOS estimator outperforms the TWIT. This can be due to the fact that the former technique used all the eigenvalues in the computation of the estimator while the latter only uses one eigenvalue to detect the number of signals. Note, however, that the TWIT has not been designed with the purpose of outperforming the SEMOS technique. Rather, the TWIT has been designed to cope with colored Gaussian noises. This is illustrated in the next figure. estimator and the TWIT. The simulation parameters are the same as the ones that have been used in previous simulations. The noise is now a colored noise. For this purpose, a white Gaussian noise with variance unity is passed through a filter with N f coefficients f i , 1 ≤ i ≤ N f given in Table II.

We have that

N f i |f i | 2 = 1
so the colored Gaussian noise at the output of the filter is also unit variance. The simulation results in Fig. 11 show that the SEMOS estimator is unable to estimate correctly the CM while the TWIT is still performing well.

A. Discussion

The simulation results show that the estimation techniques perform better when the number of observations T is much larger than the number of signals K. This is in stark contrast with previous studies, such as the ones implemented in [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF] and [START_REF] Özgül | Wireless Access with Blind Collision-Multiplicity Detection and Retransmission Diversity for Quasi-Static Channels[END_REF] where T is set to a value not greater than K + 2. In other approaches [START_REF] Zheng | Multipacket Reception in Wireless Local Area Networks[END_REF], [START_REF] Huang | Cross-Layer Multi-Packet Reception Based Medium Access Control and Resource Allocation for Space-Time Coded MIMO/OFDM[END_REF], a single observation is required to perform the MUD of the colliding users. These references do not explicitly claim that the number of signals is estimated using only one observation. In this context, the blind separation technique designed in [START_REF] Talwar | Blind Separation of Synchronous Co-Channel Digital Signals Using an Antenna Array. Part I. Algorithms[END_REF], is implemented and the number of sources (users) is assumed to have been estimated using an MDL-like criterion.

VI. CONCLUSION

In this paper, two new CM estimation techniques, the TWIT and the SEMOS estimator, have been proposed. The methods are based on eigenvalue statistics. In the TWIT, eigenvalues are tested in descending order, from the largest to the lowest. The first eigenvalue λq that is considered as being Tracy-Widom distributed allows the CM estimation by K = q -1. In the SEMOS estimator, a criterion is computed just as in the MDL criterion approach. These CM estimation techniques have been shown to outperform the typical MDL criterion. The SEMOS estimator outperforms the TWIT in the presence of a white Gaussian noise but it is inefficient in the presence of a colored noise. Moreover, simulation results have shown that a large number of snapshots T is needed in order to allow a good estimation of K in terms of MSE. Furthermore, the number of snapshots must be significantly higher than the number of colliding users K (T ≫ K). These settings are similar to the settings that are used in MOS techniques for signal array processing.

The impact of these results is twofold. First, some CR techniques such as the network-assisted diversity multiple access (NDMA) [START_REF] Zhang | Collision Resolution in Packet Radio Networks Using Rotational Invariance Techniques[END_REF], [START_REF] Özgül | Wireless Access with Blind Collision-Multiplicity Detection and Retransmission Diversity for Quasi-Static Channels[END_REF] cannot be implemented in IEEE 802.11 networks notably because these CR techniques are based on the assumption that T can be made as small as K +1 or K +2. Second, some multiple packet reception (MPR) protocols for IEEE 802.11 networks that use the blind user separation in [START_REF] Talwar | Blind Separation of Synchronous Co-Channel Digital Signals Using an Antenna Array. Part I. Algorithms[END_REF] appear to be rather questionable since they assume that the CM estimation can rely on a single observation of collided request-to-send (RTS) frames. Even if the destination node is equipped with four antennas (T = 4), our simulation results have shown that the receiver at the destination node needs many more snapshots in order to provide a good estimation of K. This paper has pointed out a strong constraint in the design of MPR techniques. It revealed that a single observation of the colliding signals is far from providing enough information to estimated the number of colliding nodes.

Further investigations are now needed in order to fully characterized the performance of the proposed CM estimation techniques in typical operating conditions. The obtained results will allow the implementation of these techniques in current or future standards.

Fig. 1 .

 1 Fig. 1. Collision Scenario with K = 3 colliding nodes.

Fig. 2 .

 2 Fig. 2. Distribution of the eigenvalues of a matrix R in the signal-free case with T = 500 and m = 2000 (λ = 1).

Fig. 3 .

 3 Fig. 3. Distribution of the eigenvalues of a matrix R † in the signal-free case with T = 500 and m = 2000.

Fig. 4 .

 4 Fig. 4. Probability density function of the largest eigenvalue of a matrix R † in the signal-free case with m = 61440 and different values for T and the noise variance σ 2 .

Fig. 5 .

 5 Fig. 5. Probability density function of the smallest signal eigenvalue of a matrix R † in the signal-bearing case with K = 4, i.e. λ5 , with m = 61440, a noise variance of σ 2 = 0.1, and different values for T . The solid curves with markers are the simulation results and the dotted curves are the Gaussian distribution.

Fig. 6 .

 6 Fig. 6. Estimates of K = 3 with the TWIT and the MDL criterion for a variable number of snapshots T , 5 ≤ T ≤ 15, and different SNR values.

Fig. 7 .

 7 Fig. 7. Estimates of K = 4 with the TWIT and the MDL criterion for a variable number of snapshots T , 5 ≤ T ≤ 15, and different SNR values.

Fig. 8 .

 8 Fig. 8. MSE of two estimation techniques (TWIT and the MDL criterion) for K = 4, a variable number of snapshots T , 5 ≤ T ≤ 15, and different SNR values.

Fig. 11

 11 Fig. 11 represents the performance of both the SEMOS

Fig. 9 .

 9 Fig. 9. MSE of two estimation techniques (the SEMOS estimator and the MDL criterion) for K = 4, a variable number of snapshots T , 5 ≤ T ≤ 15, and different SNR values.

Fig. 10 .

 10 Fig. 10. MSE of two estimation techniques (the SEMOS estimator and the TWIT) for K = 4, a variable number of snapshots T , 5 ≤ T ≤ 15, and different SNR values.

Fig. 11 .

 11 Fig. 11. MSE of two estimation techniques (the SEMOS estimator and the TWIT criterion) for K = 4, a variable number of snapshots T , 5 ≤ T ≤ 15, and different SNR values, in the presence of a colored noise.

  are circularly symmetric Gaussian coefficients with zero mean and power unity and the coefficients are randomly changing from one observations to another. User stations are transmitting OFDM signals. The signals are constructed according to the

  Number of OFDM symbols per OFDM block N OFDM 48

	TABLE I
	OFDM SIGNAL PARAMETERS
	Modulation	BPSK
	Number of sub-carriers N sub	1024
	Guard Interval GI	1/4

In infrastructure-based networks, multiplexing techniques allow a fair bandwidth allocation among the users without any risk of collision.

There are scenarios in which a collision occures between several source nodes transmitting toward several destination nodes. These scenarios are not addressed in this paper

Note that signal samples in Code Division Multiple Access (CDMA) systems are also uniformly distributed since modulation symbols are multiplied by +1/ -1 spreading codes.

APPENDIX

A. Complex Wishart Matrices Property 6.1: Consider a T × m matrix X with m samples drawn from a T -dimensional complex Gaussian law with zero mean and covariance matrix R X , denoted CN T (0, R X ). The random matrix XX H is a T -variate complex Wishart matrix with m degrees of freedom [START_REF] Johnstone | High Dimensional Statistical Inference and Random Matrices[END_REF]. This is denoted

where CW T (m, R X ) denotes complex Wishart law, parameterized accordingly.

B. Empirical Spectral Distribution (ESD) of a Matrix

The ESD of an T ×T Hermitian matrix X, denoted F X (x), is the distribution function of the eigenvalues of X, i.e., for

where λ 1 , . . . , λ T are the population eigenvalues of X and 1 A (x) is the indicator function of the set A, i.e., 1 A (x) = 1 if x ∈ A, and 1 A (x) = 0 otherwise. The Hermitian property is required to ensure that all eigenvalues of X belong to the real line. Property 6.2: Let X ∈ C T ×m be a random matrix with independent and identically distributed (iid) entries X ij such that X ij has zero mean and variance λ.

converges almost surely to a non-random distribution function with density f c (x) given by

The probability density function (pdf) f c (x) is the Marcenko-Pastur density. Property 6.3: Let B ∈ C T ×m be a random matrix defined as

where X ∈ C T ×m is a random matrix with iid entries that have zero mean and power unity, and where T is an T × T Hermitian non negative definite matrix. As T, m → ∞ with T /m → c ∈ (0, ∞), the ESD of B converges almost surely to a non-random distribution function with density f c (x) given by

where T /N → c 1 ∈ (0, 1) as T → +∞.

The probability density function (pdf) f c (x) is a modified Marcenko-Pastur density.

C. Distribution of the Largest Eigenvalue of Wishart Matrices

We now investigate the distribution of the largest eigenvalue of Wishart matrices. Consider A ∼ CW T (m, R b ) and B ∼ CW T (N, R b ), two independent random matrices where m, N > T . We consider the generalized eigenproblem constructed from A and B

where v denotes the eigenvector corresponding to the eigenvalue θ. This is referred to as a double Wishart setting [START_REF] Johnstone | High Dimensional Statistical Inference and Random Matrices[END_REF]. Property 6.4: The largest eigenvalue satisfying [START_REF] Wax | Detection of Signals by Information Theoretic Criteria[END_REF], denoted λ (D) 1 , is distributed according to a Tracy-Widom (TW) distribution when m, N → ∞ as T → ∞ with m, N > T , i.e., P W (λ

where W (θ) denotes the logit transformation of θ, i.e., W (θ) = log[θ/(1θ)] and T W C (x) is the TW distribution function for complex data.

The centering µ(T, N, m) and rescaling σ(T, N, m) in ( 18) are parameterized as follows

Note that the covariance matrix R b has no effect on the distribution of the eigenvalue λ (D)

1 . A proof of Property 6.4 can be found in [START_REF]Multivariate analysis and Jacobi ensembles: Largest eigenvalue[END_REF] for R b = I T . We now rewrite [START_REF] Wax | Detection of Signals by Information Theoretic Criteria[END_REF] as

and derive a new property for the following generalized eigenproblem Av = θBv [START_REF] Johnstone | High Dimensional Statistical Inference and Random Matrices[END_REF] This is referred to as a single Wishart setting. From ( 19) and [START_REF] Bartlett | A note on the multiplying factors for various χ 2 approximations[END_REF], we get the following property.