Zhen Zhao

Caishan Liu

Wei Ma

Bin Chen

Experimental Investigation of the Painlevé Paradox in a Robotic System

Keywords: Painlevé paradox, robotic system, instability phenomenon, impulsive dynamics

 to determine the post-impact velocities of the system, and use an event-driven algorithm to perform numerical simulations. The qualitative comparisons between the numerical and experimental results are carried out and show good agreements. This study not only presents an experimental support for the shock assumption related to the problem of the Painlevé paradox, but can also find its applications in better understanding the instability phenomena appearing in robotic systems.

Introduction

It is well-known that the rigid body model for mechanical systems with unilateral constraints and friction may possess some singularities at which the dynamical equations will have multiple solutions or even no solution at all. The classical Painlevé example, where a planar slender rod slides on a rough surface, represents the simplest system with such singularities. Recently the interest in understanding the physical phenomena corresponding to the singular situations has witnessed a substantial increase . Rich information and a good overview on the subject can be found in the excellent book written by Brogliato [START_REF] Brogliato | Nonsmooth Mechanics[END_REF], which contains a wide variety of the problem-of-interest and a long list of references.

When Coulomb friction is coupled to unilateral constraints, rigid body models may have no solution for certain configurations. An important viewpoint adopted by many authors is that a shock should then exist at the contact point. Since the shock occurs in a special situation without normal velocity and friction is considered to be the main cause for its occurrence, different nomenclatures can be found in the literature, such as the impact without collision (IW/OC), frictional catastrophe, or tangential impact [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Wang | Two-dimensional rigid-body collisions with friction[END_REF][START_REF] Génot | New results on Painlevé paradoxes[END_REF][START_REF] Brach | Impact Coefficients and Tangential Impacts[END_REF]. According to the shock assumption, some crucial results can be deduced and the problem of the Painlevé paradox seems to be solvable. For example, recent development related to the time-stepping numerical method indicates that the singularity of rigid body model can be successfully avoided if the contact forces is allowed to be impulsive [START_REF] Stewart | Rigid-Body Dynamics with Friction and Impact[END_REF][START_REF] Stewart | Convergence of a Time-Stepping Scheme for Rigid-Body Dynamics and Resolution of Painleve's Problem[END_REF][START_REF] Anitescu | Time-stepping schemes for stiff multi-rigid-body dynamics with contact and friction[END_REF][START_REF] Anitescu | Time-stepping for three-dimensional rigidbody dynamics[END_REF], thus confirming a fact observed numerically in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]. Obviously the shock assumption is fundamental for the problem of Painlevé paradox, and the validation from experiments should be introduced. One aim of this paper is to develop an experimental method to support the shock assumption.

The Painlevé phenomenon was firstly discovered in the classical Painlevé's example (a slender rod that slides on a rough plane), and many excellent theoretical results are motivated from the simple system. However, using the example to serve as the experimental model creates a lot of difficulties. For instance, the coefficient of friction in the system has to be greater than 4/3 for the occurrence of the Painlevé paradox, a large value rarely found in practical materials [START_REF] Génot | New results on Painlevé paradoxes[END_REF]. Meanwhile, the sliding of a rod under gravity is also difficult to be implemented in practice. Therefore, it is crucial to search for a new example that can involve the paradox and can be easily implemented by experiments. Such examples have been presented in [START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF][START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF], where it has been shown that the Painlevé paradoxes may occur for arbitrarily small values of the coefficient of friction.

Recent studies have indicated that the Painlevé paradox may be a common phenomenon and could be found in a variety of different applications such as robotic manipulation, legged locomotion, and vehicle braking systems [START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF][START_REF] Wilms | The occurrence of Painlevé's paradox in the motion of a rotating shaft[END_REF][START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos. Part ii: Dynamics and modeling[END_REF][START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF]. Especially, according to our recent work for a robotic system that comes into contact with a moving belt [START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF], the Painlevé paradox will appear even though the coefficient of friction takes a very small value. In this paper, we will use the analytical model presented in [START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF] to develop an experimental setup for the demonstration of the dynamical behaviors associated with the Painlevé paradox.

According to the theoretical analysis shown in [START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF], the Painlevé paradox can be found by setting the robotic system within a paradoxical configuration, which depends on the practical value of the coefficient of friction. By conducting the shock assumption into the paradoxical situation of the robotic system, and using the Darboux-Keller's method [START_REF] Zhao | The Numerical Method for Three-Dimensional Impact with Friction of Multi-Rigid-Body System[END_REF][START_REF] Kane T R | Dynamics: Theory and Applications[END_REF][START_REF] Stronge | Swerve during three-dimensional impact of rough rigid bodies[END_REF][START_REF] Keller | Impact with friction[END_REF][START_REF] Bhatt | Partitioning the parameter space according to different behaviors during three-dimensional impacts[END_REF][START_REF] Batlle | Rough balanced collisions[END_REF][START_REF] Darboux | Etude Géométrique sur les percussions et le choc des corps[END_REF], it is also shown that the tangential impact will take a peculiar property; the relative motion of the contact point must be immediately brought into a stick in tangential direction. Therefore, the shock assumption can be verified by observing whether there is a tangential stick in paradoxical situations. In our experiments, two laser-vibrameters with precise accuracy will be used to carefully measure the changes of the contact relative velocities of the robotic system.

Using numerical methods to reproduce the dynamical behaviors is significant for the analysis of mechanical systems. It is obvious that a jump rule related to the tangential impact should be provided for systems involving paradoxical situations. Since the tangential impact will excite a bouncing motion that can make the system contain differ-ent modes of motion such as slipping phases, collisions with friction, and flying without contact, an event-driven algorithm will be used to perform the numerical simulations. Qualitative comparisons between the numerical and experimental results are carried out.

Based on the numerical analysis and the experimental observation, two instability phenomena exhibited in robotic systems are discovered: one that is the bouncing motion induced by a tangential impact, the other that is related to a slip-stick motion that appears in the robotic system without the occurrence of the Painlevé paradox. The first case indicates that the problem of the Painlevé paradox will extremely influence the controllability of the robotic system, as pointed out by Brogliato [START_REF] Brogliato | Some Perspectives on the Analysis and Control of Complementarity Systems[END_REF]. The second case shows that the stick-slip motion is periodic when the rail moves with a constant velocity.

The organization of the study is as follows. Section 2 presents the description of the model for a two-link manipulator. The key theoretical results developed in [START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF] for the paradoxical situation will be reviewed in section 3. The experimental setup and the phenomena associated with the paradoxical situations will be exhibited in section 4. Section 5 provides the comparison between the numerical and experimental results. We conclude in section 6 with a summary and the potential application of our study.

A two-link manipulator with unilateral constraint

This section will first conduct the model of a two-link robotic system, and then find the condition for the occurrence of the Painlevé paradox by using an LCP (Linear Complementarity Problem) method, a theory for non-smooth dynamics established by Moreau [START_REF] Moreau | Les liaisons unilatérales rt le principe de Gauss[END_REF][START_REF] Moreau | Mécanique Classique, tome II[END_REF] and then extended into the multibody systems by Pfeiffer and Glocker [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF].

The manipulator is shown in figure .1, which consists of two identical rods with length We set an inertial coordinate frame Oxy attached at the joint O, and suppose that a local inertial coordinate frame (B, t, n) with origin located at contact point B, is defined such that n is normal to the contact surface and (t, n) forms a right-handed coordinate system. (F t , F n ) represent the contact forces in the tangential and normal direction. The components in Oxy for the contact point B can be expressed by the generalized coordinates.

x = x t x n = l(sin θ 1 + sin θ 2 ) -l(cos θ 1 + cosθ 2 ) (1)
These kinematics can yield the contact Jacobian matrix K that relates velocities and accelerations of the contact point B to the generalized coordinates through the relations

ẋ = K T q
(2)

ẍ = K T q + S (3)
Where

K = K 1 K 2 = l cos θ 1 l sin θ 1 l cos θ 2 l sin θ 2 , S = S 1 S 2 = -l( θ2 1 sin θ 1 + θ2 2 sin θ 2 ) l( θ2 1 cos θ 1 + θ2 2 cos θ 2 )
The governing equations for the system with persistent contact can be written as

q = M -1 KF + M -1 (-R + W) (4)
Where

M = 4ml 2 /3 ml 2 cos(θ 1 -θ 2 )/2 ml 2 cos(θ 1 -θ 2 )/2 ml 2 /3 F = F t F n q = θ1 θ2 W = τ 1 -τ 2 -3mgl sin θ 1 /2 τ 2 -mgl sin θ 2 /2 R = ml 2 θ2 2 sin(θ 1 -θ 2 )/2 ml 2 θ2 1 sin(θ 1 -θ 2 )/2
The substitution of Equation (3) into Equation (4) leads to

ẍ = QF + K T M -1 (-R + W) + S (5)
where

Q = K T M -1 K = Q 11 Q 12 Q 21 Q 22 ( 6
)
is a matrix that depends only on the configuration of the system.

Fact 1. Q is a symmetric and positive definite matrix, since M is a symmetric positive matrix and K is full-rank in most cases (iE; K may have singularities only in some extreme configurations).

The relative velocity between the contact point and the belt is

ẋr = x t -v t = l( θ1 cos θ 1 + θ2 cos θ 2 ) -v t (7)
If ẋr = 0, the manipulator will slip on the moving belt, otherwise the tip sticks on the belt. Defining a velocity-dependent coefficient of friction µ, in which µ = µ 0 for ẋr < 0 and µ = -µ 0 for ẋr > 0, the relationship between the tangential and the normal contact force by the Coulomb's frictional law can be expressed as

F t = µ • F n (8)
During the slip mode, the dynamical equations in tangential and normal directions take the following form.

ẍn = A(q, µ) • F n + B(q, q) (9) ẍt = C(q, µ) • F n + D(q, q) (10)
where

A(q, µ) = µQ 21 + Q 22 , B(q, q) = K T 2 M -1 (-R + W) + S 2 C(q, µ) = µQ 11 + Q 12 , D(q, q) = K T 1 M -1 (-R + W) + S 1
Combining the Signorini complementarity condition ( ẍn ≥ 0, F n ≥ 0 and ẍn • F n = 0) with equation ( 9) gives the standard formulation of an LCP whose unknown is F n and whose matrices (here a scalar) are A(q, µ) and B(q, q). Obviously, negative values of A will make this LCP possess multiple solutions or no solution at all. As illustrated by many authors a so-called impact without collision occurs because of the configurations in which A < 0 and B < 0.

More interestingly, by observing the ingredients of the coefficient A, we can find that the paradoxical situation just depends on the configuration of the system and the velocitydependent coefficient of friction. Therefore, for a given coefficient of friction, the Painlevé paradox is allowed to occur only when the system takes the paradoxical configurations, for which A < 0. In other words, we can rely on the practical value of the coefficient of friction to determine the initial configuration that can make the Painlevé paradox appear.

3 The properties of the tangential impact and the impact rule

By conducting the shock assumption into the paradoxical situations, the theoretical description related to the properties of the tangential impact and the impact rule for the robotic system have been presented in [START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF]. In this section some key results will be introduced.

The properties of the tangential impact

In order to consider the coupling between normal and tangential motion, the experience of an impact with friction must be carefully investigated. By assigning to the impact duration a very short but not infinitesimal time, Darboux [START_REF] Darboux | Etude Géométrique sur les percussions et le choc des corps[END_REF] and Keller [START_REF] Keller | Impact with friction[END_REF] developed a method that yields a set of differential equations with respect to the normal impulse, a 'time-like' independent variable. These nonlinear differential equations describe the impulsive behaviors (the impact dynamics), such that the singularities of impacts due to friction can be successfully avoided [START_REF] Bhatt | Partitioning the parameter space according to different behaviors during three-dimensional impacts[END_REF][START_REF] Batlle | Rough balanced collisions[END_REF]. This method has been extended for the investigation of the properties of the tangential impact [START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF].

Let us set the impact duration as [t 0 , t f ] and divide this short time into much smaller intervals [t i , t i+1 ]. Integrating equation ( 5) and ignoring the contribution of the finite forces on [t i , t i+1 ], yields the following differential equations of motion:

     d ẋt = Q 11 dP t + Q 12 dP n d ẋn = Q 21 dP t + Q 22 dP n (11)
where

dP t = t i+1 t i F t dt, dP n = t i+1 t i F n dt
are the changes of tangential and normal impulses on [t i , t i+1 ], respectively. Now applying the relationship dP t /dP n = µ defined by the Coulomb's frictional law into equation ( 11) directly leads to

     d ẋt = (µQ 11 + Q 12 )dP n d ẋn = (µQ 21 + Q 22 )dP n (12)
where

d ẋt = t i+1 t i ẍt dt, d ẋn = t i+1 t i ẍn dt
are the changes of the tangential and normal velocities, respectively.

During the interval [t i , t i+1 ], the normal impulse is a strictly monotone function of time. This permits to think of it as a 'time-like' independent variable, and to perform a time-scale of the shock dynamics. Thus, equation ( 12) is a set of first order ordinary differential equations with respect to dP n , which varies like a time variable. Combining equation [START_REF] Wang | Two-dimensional rigid-body collisions with friction[END_REF] with the condition for the occurrence of the Painlevé paradox, we can deduce the property of the tangential impacts, which is elucidated by the following theorem.

Theorem.1 The impulsive process induced by the Painlevé paradox will first result in a normal compressional process, and then immediately bring the relative tangential velocity of the contact point to zero. After that the tangential motion of the contact point will stick on the contact surface, while the normal motion of the contact point will continue to be compressional until the normal velocity equals zero. Then an expansion phase in the normal direction is carried out to make the impact finish.

Proof:

We will use fact 1 for the properties of matrix Q, the condition for the occurrence of Painlevé paradox, and the Coulomb's frictional law to prove the property of the tangential impact.

Matrix Q can be thought of as a constant matrix due to the little change of the configuration in an impulsive process. So the elements in Q satisfy the following relationships based on fact 1:

Q 22 > 0, Q 11 > 0, Q 11 Q 22 > Q 2
During the slip mode, the condition for the occurrence of Painlevé paradox permits to write that

A = µQ 21 + Q 22 < 0, -µQ 21 > Q 22
In the case ẋt < 0, we have µ = µ 0 > 0, so Q 21 < 0. Therefore,

-µQ 21 Q 11 > Q 22 Q 11 > Q 2 21 and (µQ 11 + Q 12 ) > 0
Similarly, if ẋt > 0, we can obtain the following inequality

(µQ 11 + Q 12 ) < 0
Thus, these two cases ( ẋt > 0 and ẋt < 0) can be expressed by using a uniform inequality:

ẋt (µQ 11 + Q 12 ) < 0
The above inequality indicates that the magnitude of the tangential velocity always decreases when the sliding friction is sustained. Since the condition A < 0 cannot be removed if the tangential velocity is not set to zero, the tangential motion of the contact point must reach the point of slip stopping.

Before slip stops, we can find from the second equation of ( 12) that the contact point will take a normal velocity penetrating into the contact surface, with an increment for its magnitude (A < 0 and the initial value equal to zero). Due to the coupling between the normal and tangential impulses, a compressional process in the normal direction is generated once the paradox appears. Before slip stopping, A < 0 is always satisfied, so that the magnitude of the normal velocity will continue to increase. So the tangential slip at the contact point stops in the compressional process where the normal velocity is not equal to zero.

Once the tangential speed vanishes, the tangential motion might stick on the contact surface or continue to slip, depending on the property of dry friction. Setting d ẋt = 0 in equation [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], one can obtain the following inequality

|dP t /dP n | = |Q 12 /Q 11 | < µ 0
Usually the static coefficient of friction, µ s is larger than the sliding coefficient µ 0 . The above inequality indicates that contact forces must enter into the interior of the friction cone once the tangential velocity disappears. So stick should occur at the instant when the tangential velocity vanishes. During the tangential impact, no other additional impulses are applied on the system. The stick mode can be preserved until the impulsive process finishes. Therefore, we can conclude that the relative tangential speed after tangential impact must equal zero. This well coincides with the maximum dissipation principle stated in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] by Moreau. Once stick occurs, dP n and dP t cannot be connected linearly by the coefficient of friction µ, and must satisfy the relationship defined by the first equation in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] by setting d ẋt = 0. This relationship can make the normal velocity of the contact point decrease and reach zero. After that an expansion phase will occur in order to release the energy accumulated in the compressional phase. The expansion phase can be governed by an impact law such as the Poisson's or Stronge's laws [START_REF] Stronge | Swerve during three-dimensional impact of rough rigid bodies[END_REF]. ✷ We can summarize the process of the tangential impact as follows: the friction will first result in a compressional motion in the normal direction and then bring the relative tangential motion from slip to stick. Then, sticking motion will persist until the contact constraint is released. The normal motion at the contact point will continue to be compressional from the instant of stick appearance to the time when the normal velocity vanishes. After that an expansion phase is carried out to make the impact terminate.

The impact rule for the tangential impact

Assigning a duration [t 0 , t f ] to the tangential impact, we can split such an impulsive process into three periods. The first one is a sliding compressional period denoted as [t 0 , t 1 ], where t 1 is related to the instant of stick appearing. The second one is a sticking compressional period denoted as [t 1 , t 2 ], where t 2 corresponds to the instant when the normal velocity vanishes. The third one is a sticking restitution period denoted as [t 2 , t f ], which describes the expansion process of the normal motion, defined by using the Poisson's law for normal impact.

Sliding compressional period, [t 0 , t 1 ]
Let us set ẋ0 t as the initial speed of slip at t 0 . At t = t 1 , the tangential velocity ẋ1 t will be equal to zero. Thus, the change of the normal impulse can be obtained by integrating the first equation of ( 12)

P 1 n = - ẋ0 t µQ 11 + Q 12 ( 13
)
By considering the initial value of the normal velocity ẋ0 n = 0, we can obtain the normal velocity ẋ1 n at t 1 by integrating the second equation in [START_REF] Wang | Two-dimensional rigid-body collisions with friction[END_REF] and by using expression (13):

ẋ1 n = (µQ 21 + Q 22 )P 1 n = - µQ 11 + Q 12 µQ 21 + Q 22 ẋ0 t ( 14
)
Sticking compressional period, [t 1 , t 2 ]
The stick at the contact point implies the following relationship:

d ẋt = 0, dP t /dP n = -(Q 12 /Q 11 ) (15)
Combining [START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF] with the second equation of [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], one can deduce the differential equation for the normal motion in the stick mode:

d ẋn = (- Q 2 12 Q 11 + Q 22 )dP n (16)
Due to ẋ2 n = 0 at the end of this period, we can easily obtain the change of the normal impulse P 2 n at t 2 :

P 2 n = Q 11 ẋ1 n Q 2 12 -Q 11 Q 22 = - Q 11 (µQ 11 + Q 12 ) (Q 2 12 -Q 11 Q 22 )(µQ 21 + Q 22 ) ẋ0 t ( 17
)
Sticking restitutional period, [t 2 , t f ] This period represents an expansion process of the normal motion. The expansion impulse P r n can be obtained by using the Poisson's coefficient e p e p = P r n P c n where P c n = P 1 n + P 2 n is the compressional impulse in the normal direction. So we have P r n = e p P c n = e p (P 1 n + P 2 n ) At the beginning of this period, the normal speed ẋ2 n = 0. Meanwhile, the normal motion during this period will be governed by equation ( 16) since stick in tangential direction is preserved. Thus, at the end of this period the normal speed ẋf n can be expressed as

ẋf n = (- Q 2 12 Q 11 + Q 22 )P r n ( 18
)
Clearly the post-impact velocity in the normal direction, ẋf n is not equal to zero except for e p = 0 for the tangential impact without any initial normal velocity. In other words, after the tangential impact the contact point will leave the contact surface with a certain velocity in the normal direction. Nevertheless, the tangential velocity at the contact point ẋf t must be equal to zero when the tangential impact finishes. By integrating equation ( 4) and neglecting the contribution of the finite forces, the changes of the generalized velocities of the system due to the tangential impact can also be calculated:

∆ θ1 ∆ θ2 = M -1 K P t P n ( 19
)
where ∆ θ1 and ∆ θ2 represent the changes of the generated velocities, respectively. P t and P n are the total impulses in the tangential and normal directions.

4 The dynamical behavior related to the paradoxical situation

In this section, we will present the experimental setup for the robotic system that corresponds to the analytical model described in above section. According to the coefficient of friction estimated from experiments, we will firstly determine the paradoxical configuration that can make the Painlevé paradox appear. By setting the system with the paradoxical configuration and initially establishing a contact constraint, we can observe the dynamical behavior associated with the paradoxical situation. The paradoxical phenomena will be demonstrated by showing the velocities of the contact point. The influence of the rail's speed on the bouncing motion generated due to the tangential impact will be exhibited. In addition, how the dynamical behaviors of the robotic system evolve into a shock from a normal configuration to the singular one will be demonstrated experimentally, and the stick-slip phenomena appearing in the robotic system will be investigated.

The description of experiment

The experimental setup of the robotic system formed by using two identical aluminium cylindrical bars (m = 0.12kg, l = 0.21m) connected with revolute joints, see figure 2. The upper revolute joint is used to connect the system with a fixed bracket that contains a slot to make the height of the system adjustable. A semi-spherical head made of plastic material is installed on the contact end in order to make the contact achieve relatively uniform conditions during the motion. The belt used in the analytical model is replaced by a steel rail, which moves along the horizonal direction, dragging a rope by hand.

A passive torque is provided by a torsional spring mounted at the middle revolute To consider the effects of joint friction, we use a uniform coefficient c to represent the damping torques acting on the two identical revolute joints. So the torques τ 1 and τ 2 at the revolute joints can be approximated as

τ 1 = -c θ1 (20)
τ 2 = k(θ 1 -θ 2 -α 0 ) + c( θ1 -θ2 ) ( 21
)
where k is the stiffness of the torsional spring, and α 0 is the initial angle of the spring. Two laser-Doppler vibrometers (OFV-303-353) with a controller(OFV-3001) are used to measure the rail speed and the velocity of the contact point, in which laser signals are sent to track the movements of the sensitive papers attached in the rail and the contact head. The experimental signals are transferred into a laptop through an A\D card with 10kHz sample rate. The sketch of the experimental system is depicted in figure 3.

By using a simple slide experiment for the homogenous contact surface, we estimate the coefficient of friction as µ = 0.6. Meanwhile, the contact constraint will allow the following geometric relationship to exist

l(cos θ 1 + cos θ 2 ) = H
where H is the height of the system (see figure 1), and θ 1 and θ 2 are the angles related to the upper joint and the middle joint, respectively. Based on the above relationship and one the condition A < 0, a function of θ 1 with respect to H and µ can be obtained for the occurrence of the Painlevé paradox. For the given coefficient of friction µ = 0.6 with different value of H, Table 1 presents the allowable scope of θ 1 for the occurrence of the Painlevé paradox.

Table 1: The paradoxical configurations in different height

The coefficient of friction µ = 0.6 θ 1 (deg) none 41.5∼48.5 33∼42. [START_REF] Brogliato | Nonsmooth Mechanics[END_REF] Since the paradoxical situations appear under the condition that the robotic system should be in a paradoxical configuration with an initial condition of slip, we can release the tip of the system with an approximately zero height on the moving rail in order to generate the slip mode. In this case, the absolute value of the velocity of the tip is equal to zero, such that a relative slip between the tip and the rail can be established. According to the property of the tangential impact, the absolute tangential velocity of the tip should immediately approach the one of the moving rail if a shock exists. In the following, we will present the experimental results for the robotic system with different configurations by changing the joint angles and the height of the system. In particular, the observation of the stick phenomena associated with the paradoxical situations will be emphasized.

The experimental results

Table 2: The configurations investigated in experiments

The height H = 0.3775m θ 1 (deg) -15 [START_REF] Klein | Zu Painlevés Kritik der Coulombschen Reibungsgesetze[END_REF] Let us set the robotic system with a fixed height H = 0.3775m, then repeat experiments for the system with joint angle θ 1 that takes different values among the scope of θ 1 ∈ (-37 • , 37 • ) (the possible values that can make the tip of the robotic system touch on the moving rail). The negative value of θ 1 corresponds to the situation where rod 1 slopes to the left side against the vertical line passing through the fixed revolute joint. Table 2 presents the configurations of the system that are investigated in experiments, in which "Y" represents the paradoxical configurations, while "N" represents the ones of non paradox appearing. Figure 4 shows the experimental curves for the rail's speed and the tangential velocity of the contact point for the system with a paradoxical configuration of θ 1 = 32 • . When the tip touches the moving rail with zero velocity at t = 4.84s, the first vertical line shown in figure 4 indicates that the tip immediately approaches the value of the rail speed. This is related to a sticking phenomenon corresponding to a tangential impact appearing in the paradoxical configuration. After this event the tip bounces on the moving rail and sequential collisions appear at t = 4.92s, 5.02s, 5.10s, etc.

We also found from experiments that the magnitude of the bouncing motion will be significantly influenced by the rail's speed. Figure 5 presents the experimental results by setting the robotic system in the same configuration as in the previous experiments, while the rail's speed is changed to v t = -0.5m/s. A tangential impact appears at the measure time t = 2.12s when the tip touches the rail, and then the subsequent impacts occur at the instant t = 2.21, 2.29, .... The comparison between figures 4 and 5 clearly shows that the magnitude of the tangential velocity is enlarged due to the increase of the rail's speed.

If the rail moves slowly, the bouncing motion induced by the tangential impact will be of low magnitude and even disappear when the velocity of the rail is lower than a

Tangential impact

Normal velocity of contact point

Tangential velocity of contact point

Velocity (m/s) Time (s)

Figure 5: The normal and tangential velocity of the contact point (H = 0.3775m, θ 1 = 32 • , v t = -0.5m/s) certain threshold. Figure 6 shows the experimental results for the robotic system with the same configurations as the previous two experiments, while setting the rail moving with v t = -0.075m/s. Clearly after the tangential impact the tip of the robotic system will stick on the moving rail, even though there is a peak in the curve of the tangential velocity because of the tangential impact. This phenomenon can be elucidated from the viewpoint of the system's energy. According to equations ( 18) and ( 19), the normal and tangential impulses are proportional to the rail's speed. So the robotic system can gain more energy from the tangential impact in the situation of the rail moving fast. Thus the sequential collisions can be enlarged to make the magnitude of the bouncing motion increase. If the rail moves much more slowly, the robotic system cannot gain enough energy to overcome the contact force generated by the gravity and the torsional spring. In this case, the bouncing motion cannot be observed experimentally.

For the system with the configuration of θ 1 = 30.5 • , the bouncing motion due to the tangential impact can also be observed by experiments when the rail takes a relatively high speed (as shown in figure 17). The experimental results related to the two cases of the system with θ 1 = 30.5, 32 • will be used in the following section to verify the numerical simulations.

When θ 1 < 29 • , the joint angle θ 2 will be greater than θ 1 if the tip can touch on the rail for the system with height H = 0.3775m, so that the initial torque of the torsional spring applied at the middle joint will change its direction and then influences the dynamical behavior of the system. Figure 7 presents the experimental results for the system with a paradoxical configuration by setting joint angle θ 1 = 25 • . When contact is established, a tangential impact appears at the contact point (the tangential velocity of the tip immediately approaches the one of the rail), and then an oscillation for the tangential motion

θ 1 = 32 • , v t = -0.075m/s)
is induced. However, the contact point doesn't leave away from the rail surface after the tangential impact and no bouncing motion can be observed. The reason for that is because the system with the configurations of θ 2 > θ 1 can only obtain very little energy from the tangential impact that is not enough to overcome the effects of the gravity and the torsional spring. Therefore, the contact point will stick on the moving rail.

Although the above experimental investigations indicate that the dynamical behavior of the system with paradoxical situations seems to be well modeled as a tangential impact as the relative tangential velocity could immediately approach zero, it is clear in physical meaning that the tangential motion with a normal configuration evolves into a shock for the paradoxical configurations should correspond to a process with gradual change in tangential velocity. Figure 8 to 11 present the experimental results obtained by setting the system with θ 1 = -15, 7, 15, 21 • , respectively. This corresponds to the situation where the initial configuration of the system gradually approaches the boundary of the singular region. From the experimental curves (except for the case of θ 1 = -15 • ), we can find that friction will decrease the slip velocity and finally bring the contact point to stick on the moving rail. In particular, the duration from the beginning of slip to the occurrence of sticking becomes shorter when the configuration of the system is near to the paradoxical situation. If the configuration is very close to the boundary of the paradoxical region (the case of θ 1 = 21 • ), the duration for the stop of the relative tangential motion is about t = 0.13s. Once the configuration of the system enters into the paradoxical region, the duration for the tangential motion of the contact point changing from slip to stick will be less than 0.01s (as shown for the case of θ 1 = 25 • ), a short timescale that can be connected with an impact process. The above experimental phenomena well agree with the mechanism of the Painlevé paradox generated due to the coupling between friction and the configuration of the system, and further confirm that the dynamical behavior of the system in the Painlevé paradox could be described by a shock.

If the initial configuration is far away the singular region, as shown in figure 8 for the system with a configuration of taking θ 1 with a negative value and in figure 12 for the robotic system with the initial configuration of H = 0.25m, neither sticking phenomena nor the bouncing motion can be found in the robotic system.

In some cases, the stick-slip phenomenon is also found in our experiments for the robotic system. By setting the initial configuration of the system to H = 0.314m and θ 1 = 50 • , figure 13 shows that the friction will bring the tip into a stick, and then make the slip resume. Such a process can be repeated to render the robotic system unstable on the contact surface.

The comparison between experimental and numerical results

The experimental results presented in the above section have confirmed that the Painlevé paradox does induce a tangential impact at the contact point, and then makes the robotic system behave in a more complex way, with slip phases, stick phases, flight without contact phases, as well as tangential and normal impacts with friction. In this section, we will use an event-driven algorithm to perform the numerical simulations.

Firstly, we will carefully estimate the parameters used for the simulations. According to property of the collision between plastic and steel materials, the coefficient of the restitution can be set as e p = 0.1. A single pendulum system (shown in figure 14 is established to obtain its frequency f , and the value of the stiffness of the torsional spring k can be calculated by using the following expression:

k = 1 3 ml 2 (2πf ) 2 - l 2 mg (22)
Based on f = 4.6Hz measured from experiments, we have k = 1.3N • m/rad. Since the damping coefficient c is difficult to obtain from the experiments, a fitting method is used for its estimation. By setting c with different values to perform the simulations, we choose the best one among the different values of c as the damping coefficient that can make the corresponding numerical simulation better coincide with the experimental results. According to the numerical experiments, we find that the change of c limited in the range of [0.003, 0.007], has little influence on the numerical results. So the damping coefficient is chosen as c = 0.005Nms/rad for the following simulations.

Roughly speaking, the dynamical behaviors of the system are governed by equations ( 4) and ( 5) with the variations of the contact forces, which depend on the mode of motion at the contact point. For instance, the contact forces can be set to zero for the flying mode, while in the case of preserved contact, the contact forces should be determined by using LCP's equations and the Coulomb's friction law. If the paradoxical situations appear in the rigid body model, the simulation can be continued by setting equation ( 4) with new initial conditions that can be obtained from the tangential impact rule expressed in [START_REF] Ivanov | Singularities in the dynamics of systems with non-ideal constraints[END_REF] and [START_REF] Brach | Impact Coefficients and Tangential Impacts[END_REF]. Similar process is also carried out for the collisions with friction [START_REF] Zhao | The Numerical Method for Three-Dimensional Impact with Friction of Multi-Rigid-Body System[END_REF], in which the changes of the velocities of the system are obtained by integrating the impulsive differential equations expressed in [START_REF] Wang | Two-dimensional rigid-body collisions with friction[END_REF]. Since there are no accumulations of events, event-driven schemes are well-suited to the numerical integration of this nonsmooth system, see e.g. [START_REF] Brogliato | Numerical simulation of finite dimensional multibody nonsmooth dynamical systems[END_REF].

Figure 15 presents the comparisons between the experimental and numerical results by setting the system to the initial configuration of H = 0.3775m and θ 1 = 32 • . It is clear that the model used in simulation can well reproduce the qualitative behaviors of the system. The discrepancies appearing in figure 10 are partly due to the unmodeled effects existing in the experimental setup, such as the vibration of the sensitive papers, the clearances in the revolute joints, and the errors of the physical parameters. Other sources of discrepancies are the non-uniformity of the rail's speed in experiments, since it is assumed to be a constant in the numerical simulations. , v t = -0.2m/s) Figure 16 shows the comparisons of the normal velocities between the experimental and numerical results. The vibration of the sensitive paper induced by the bouncing motion will much influence the accuracy of the measurements. Nevertheless, the qualitative behavior of the system can still be captured through the numerical simulations. Keeping the system with the same height H = 0.3775m as in the previous case, we can change the initial configuration by adjusting the joint angle from θ 1 = 32 • to 30.5 • . The results obtained from the simulation also agree well with the experimental results (shown in figure 17).

As mentioned in the above subsection, the phenomenon of stick-slip motion can be found in experiments even though there is no paradox appearing. By setting the initial configuration of the system with H = 0.314m and θ 1 = 55 • , figure 18 presents the curves corresponding to the tangential velocities of the contact point obtained from simulations and experiments, respectively. The numerical results indicate that the stick-slip is a pure periodic motion when the rail speed takes a constant value, but this characteristic will be slightly destroyed due to the non-uniformity of the rail's motion. Summarizing the comparisons between the numerical and experimental results presented in above, we can conclude that the dynamical behaviors of the system can be well captured qualitatively by using the rigid body model. Even in the paradoxical situation, the simulations can be led by using the tangential impact rule to reinitialize the dynamical equations. The numerical results show that such a tangential impact induced by the Painlevé paradox can be well governed by the impact rule presented in [START_REF] Liu | The Bouncing Motion Appearing in a Robotic System with Unilateral Constraint[END_REF].

Summary and conclusion

This study mainly developed an experimental setup to demonstrate the phenomenon of the Painlevé paradox that appears in a robotic system coming into contact with a moving rail. According to the experimental results, we have verified that a shock is truly related to the problem of the Painlevé paradox, and takes a particular property that a tangential stick appears at the contact point.

Two kinds of instability phenomena for the robotic system are observed from the experiments. The first one is induced by the tangential impact, in which the robotic system will bounce or stick on the moving rail, depending on the value of the rail's speed. The other form of the instability is the stick-slip motion, which appears in the system does not involve the Painlevé paradox.

Based on the careful estimation of the physical parameters, an event-driven algorithm is used to perform numerical simulations. By setting the robotic system with paradoxical situations, comparisons between the numerical and experimental results are carried out and show good agreement. This illustrates that the rigid body model can well reproduce the complexly dynamical behaviors, and the Painlevé paradox can be overcome by using the tangential impact rule obtained from the Darboux-Keller's shock dynamics. The present work not only provides a basis for the theoretical results associated with the problem of the Painlevé paradox, but also may be useful for the design of feedback controllers in robotic systems.

Figure 1 :

 1 Figure 1: Two-link manipulator contacting with a constantly moving belt l and mass m and comes into contact with a moving belt with velocity v t . The external torques τ 1 and τ 2 are applied on the joints O and A. H is the height from the fixed point O to the rough surface. The joint angles θ 1 and θ 2 are selected as the generalized coordinates of the system when unconstrained by contact, and their positive values are assigned along the counter-clockwise direction.We set an inertial coordinate frame Oxy attached at the joint O, and suppose that a

Figure 2 :

 2 Figure 2: The physical model of the experimental setup joint. To consider the effects of joint friction, we use a uniform coefficient c to represent the damping torques acting on the two identical revolute joints. So the torques τ 1 and τ 2 at the revolute joints can be approximated as

Figure 3 :

 3 Figure 3: The sketch of the experimental system

  15 21 25 30.5 32 Paradox N N N N Y Y Y

Figure 4 :

 4 Figure 4: The relative velocity of the contact point in tangential direction ( H = 0.3775m, θ 1 = 32 • ,v t = -0.16m/s)

Figure 6 :

 6 Figure 6: The tangential velocity of the contact point (H = 0.3775m, θ 1 = 32 • , v t = -0.075m/s)

Figure 7 :

 7 Figure 7: The tangential velocity of the tip with H = 0.3775m and θ 1 = 25 •

Figure 8 :

 8 Figure 8: The tangential velocity of the tip with H = 0.3775m and θ 1 = -15 •

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: The tangential velocity of the tip with H = 0.3775m and θ 1 = 7 •

Figure 12 :Figure 13 :Figure 14 :

 121314 Figure 12: The tangential velocity of the tip without the Painlevé paradox (H = 0.25m, θ 1 = 69.4 • )

Figure 16 :

 16 Figure 16: Experimental and numerical results for the normal velocity (H = 0.3775m, θ 1 = 32 • , v t = -0.24m/s)

Figure 17 :

 17 Figure 17: The comparison of the tangential speed between experimental and numerical results(H = 0.3775m, θ 1 = 30.5 • , v t = -0.2m/s)

Figure 18 :

 18 Figure 18: Comparison between the experimental and numerical results for the Stick-Slip Phenomena under the configuration taken as H = 0.314m and θ 1 = 55 •

Simulation: Tangential velocity of contact point Experiment: Tangential velocity of contact point Experiment: Rail velocity Tangential impact Collision Figure

  15: Experimental and numerical results for the tangential speed (H = 0.3775m, θ 1 = 32

		0.8
		0.6
		0.4
	Velocities (m/s)	-0.2 0 0.2
		-0.4
		-0.6
		-0.8	1.8	2	2.2	2.4	2.6	2.8	3	3.2	3.4	3.6	3.8
								Time (s)

Acknowledgements

We would like to express our thanks to Dr Bernard Brogliato for his valuable comments and the suggestions to this paper. We also appreciate the assistance from Dr. Zhili Sun, Jiyan Huo and Ming Qiang to the work of experiment. The research was supported by Chinese National Science Foundation (Grant No. 60334030, 10642001, 10772002).


Linear projection equations arise in numerous fields and play a significant role in mathematical programming. The linear complementarity problem [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF] (LCP) u ≥ 0, (M u + q) ≥ 0, u T (M u + q) = 0 is equivalent to a special (LPE) with Ω = {u ∈ R n | u ≥ 0} [15]. Various constrained least squares problems [START_REF] Stoer | On the numerical solution of constrained least squares problems[END_REF] and convex quadratic programming problems can be tranlated into a linear projection equation [START_REF] Bruck | An iterative solution of a variational inequality for certain monotone operators in Hilbert space[END_REF] in which Ω is a general orthant [START_REF] He | A projection and contraction method for a class of linear complementarity problem and its application in convex quadratic programming[END_REF].

The complementarity problem has been studied starting with the works of Cottle, Dantzig [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF], and Lemke [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF][START_REF] Lemke | Equilibrium points of bimatrix games[END_REF] and has been developed by many others. There is already a substantial number of algorithms for solving linear projection equations [START_REF] Bruck | An iterative solution of a variational inequality for certain monotone operators in Hilbert space[END_REF][START_REF] Dafermos | An iterative scheme for variational inequalities[END_REF][START_REF] Fang | An iterative method for generalized complementarity problems[END_REF][START_REF] Harker | A damped-Newton method for the linear complementarity problem[END_REF][START_REF] He | A projection and contraction method for a class of linear complementarity problem and its application in convex quadratic programming[END_REF][START_REF] He | A globally linearly convergent projection and contraction method for a class of linear complementarity problems[END_REF][START_REF] He | A new method for a class of linear variational inequalities[END_REF][10][START_REF] Kojima | A polynomial-time algorithm for a class of linear complementarity problems[END_REF][START_REF] Mangasarian | Solution of symmetric linear complementarity problems by iterative methods[END_REF][START_REF] Mizuno | A new polynomial time algorithm for a linear complementarity problems[END_REF][START_REF] Pang | Iterative methods for variational and complementarity problems[END_REF][START_REF] Pang | Variational inequality problems over product sets: applications and iterative methods[END_REF][START_REF] Stoer | On the numerical solution of constrained least squares problems[END_REF], especially for linear complementarity problems and linear constrained least squares problems. Let Ω * = {u * | u * is a solution of (LPE)} be the solution set of ( 1) and [START_REF] Dafermos | An iterative scheme for variational inequalities[END_REF] e(u) := u -P Ω [u -(M u + q)].

be the "error" by which a given point u fails to satisfy [START_REF] Bruck | An iterative solution of a variational inequality for certain monotone operators in Hilbert space[END_REF]. In the projection and contraction method of [START_REF] He | A projection and contraction method for a class of linear complementarity problem and its application in convex quadratic programming[END_REF][START_REF] He | A globally linearly convergent projection and contraction method for a class of linear complementarity problems[END_REF], the vector ( 5)

g(u) = M T e(u) + (M u + q)
is used as the search direction. The recursion (6)

u k+1 = P Ω [u k -ρ(u k )g(u k )]
with [START_REF] He | A projection and contraction method for a class of linear complementarity problem and its application in convex quadratic programming[END_REF] ρ(u) = e(u) 2 (M T + I)e(u) 2 produces a sequence {u k } ⊂ Ω, which satisfies

(8) u k+1 -u * 2 ≤ u k -u * 2 -ρ(u k ) e(u k ) 2 .
The main advantages of this method are its simplicity and ability to handle the linear projection equation (1) while some other algorithms (e.g. [START_REF] Luenberger | Introduction to Linear and Nonlinear Programming[END_REF]) can only solve special cases of [START_REF] Bruck | An iterative solution of a variational inequality for certain monotone operators in Hilbert space[END_REF]. Each iteration of this method consists essentially of only two matrix vector products and two projections of a vector onto Ω. Therefore the method allows the optimal exploitation of the sparsity of the matrix M and may thus be efficient for large sparse problems [10]. Since the method is easy to parallelise, it may be even more favorable for parallel computation. However, for ill-conditioned problems, the search direction (5) may lead to a very slow convergence. Our objective in this paper is to find better search directions and thereby to construct more efficient methods for solving problem [START_REF] Bruck | An iterative solution of a variational inequality for certain monotone operators in Hilbert space[END_REF]. Throughout this paper we assume that Ω * / = ∅ and that the projection onto Ω is simple to carry out (e.g. when Ω is a general orthant, a box, a sphere, a cylinder or a subspace).

The paper is organized as follows. In Sect. 2 we illustrate our motivation. The main theorem is proved in Sect. 3. In Sect. 4 some new methods are presented and their contractive properties are shown. In Sect. 5 we prove convergence. Finally, in Sect. 6, we give some extensions and conclusions.

We use the following notation. For u ∈ R n , the component i is denoted by u i .

A superscript such as in u k refers to a specific vector and usually denotes iteration index. The Eucliden norm and the max-norm will be denoted by • and • ∞ , respectively. Throughout this paper, G denotes a positive definite matrix and u G denotes (u T Gu) 1 2 .

Motivation

In order to illustrate our motiviation, we let

f (u) = 1 2 u T M u + q T u
and M be symmetric positive definite, and consider the following unconstrained optimization problem

(9) min u∈R n f (u).
Solving problem ( 9) is equivalent to finding a zero point of ∇f (u). The search direction in classical methods for unconstrained optimization is

(10) d(u) = Q∇f (u)
with different matrices Q. If Q = I, we obtain the direction of the steepest descent method. Setting Q = [∇ 2 f (u)] -1 yields the direction of Newton's method. When

Q = σI + (∇ 2 f (u)) -1 or Q = [σI + ∇ 2 f (u)] -1
for some nonegative value of σ > 0, the search direction can be regarded as some combination of steepest descent (σ very large) and Newton's method (σ = 0). For the convex constrained optimization problem

(11) min{f (u) | u ∈ Ω}, the Kuhn-Tucker theorem tells us that u * is a minimum if u * ∈ Ω and it satisfies (u -u * ) T ∇f (u * ) ≥ 0 for all u ∈ Ω.
This means that u * is a zero point of the function e(u). Note that in the case of

Ω = R n , e(u) = ∇f (u).
Since most search directions in unconstrained optimization are constructed from ∇f (u), a natural question is whether we can build useful search directions for the constrained optimization problem [START_REF] Kojima | A polynomial-time algorithm for a class of linear complementarity problems[END_REF] based on e(u). Further, for problem (1), if we take

(12) d(u) = Qe(u)
as the search direction, which step length should be taken?

The main theorem

The following theorem plays an important role in our new methods.

Theorem 1. Let u * ∈ Ω * . Then (13) (u -u * ) T (I + M T )e(u) ≥ e(u) 2 + (u -u * ) T M (u -u * ) ∀u ∈ R n .
Proof. Since Ω ⊂ R n is a closed convex set and u * ∈ Ω, we know by the properties of a projection on a closed convex set [14, Appendix B] that

{v -P Ω (v)} T {P Ω (v) -u * } ≥ 0 ∀v ∈ R n .
By setting v := u -(M u + q) we obtain ( 14)

{e(u) -(M u + q)} T {P Ω [u -(M u + q)] -u * } ≥ 0. Since P Ω (•) ∈ Ω, it follows from (2) that (15) (M u * + q) T {P Ω [u -(M u + q)] -u * } ≥ 0.
Adding ( 14) and ( 15) we get ( 16)

{e(u) -M (u -u * )} T {u -u * -e(u)} ≥ 0 and it follows that (u-u * ) T (I + M T )e(u) ≥ e(u) 2 + (u -u * ) T M (u -u * ).
A similar but weaker result of Theorem 1 was given in [START_REF] He | A new method for a class of linear variational inequalities[END_REF]. The above proof is an improved version of the one in [START_REF] He | A new method for a class of linear variational inequalities[END_REF]. We point out that inequality [START_REF] Mizuno | A new polynomial time algorithm for a linear complementarity problems[END_REF] is sharp. This can easily be seen by setting M = I. We then obtain {e(u)-(u-u * )} T {u-u * -e(u)} = 0. This implies that also the result of Theorem 1 is tight.

Remark. The methods in [START_REF] He | A projection and contraction method for a class of linear complementarity problem and its application in convex quadratic programming[END_REF] and [START_REF] He | A globally linearly convergent projection and contraction method for a class of linear complementarity problems[END_REF] take g(u) as the search direction. It was shown that

(u -u * ) T g(u) ≥ e(u) T (M u + q).
But only under the assumption that u ∈ Ω can we prove

e(u) T (M u + q) ≥ e(u) 2 . Therefore, -g(u) is a descent direction of F (u) = 1 2 u -u * 2 at u ∈ Ω. However,
here the assertion (13) in Theorem 1 is true for all u ∈ R n . Although, as in [START_REF] He | A new method for a class of linear variational inequalities[END_REF], -(I + M T )e(u) can be taken as a descent direction of F (u) for all u ∈ R n , Theorem 1 offers us the possibility to construct better search directions.

The methods and their contractive properties

In this section, based on Theorem 1, we give some new methods for solving linear projection equations and show their contractive properties. The iterative scheme of these methods is ( 17)

u k+1 = u k -ρ(u k )Qe(u k )
with different matrices Q and steplengths ρ(u).

Method 1.

(for symmetric M ≥ 0)

Q = I, ρ(u) = e(u) 2 e(u) T (I + M )e(u)
.

Method 2.

(for symmetric M > 0)

Q = M -1 , ρ (u) = e(u) 2 e(u) T (I + M -1 )e(u)
.

Method 3.

(for symmetric M > 0)

Q = I + M -1 , ρ (u k ) = e(u) 2 (I + M -1 )e(u) 2 M . Method 4.
(for M ≥ 0 but not necessarily symmetric)

Q = (I + M ) -1 , ρ (u) = 1.
The first method can be viewed as an extention of the steepest descent method for unconstrained optimization, because we take e(u) as the search direction and e(u) is the residue of the projection equation. Obviously, each iteration of this method consists essentially of only a projection to Ω and the computation of M u and M e(u).

The second method can be viewed as an extention of Newton's method for unconstrained optimization. As in Method 1, each iteration of this method consists essentially of only a projection to Ω and the computation of M u and M -1 e(u).

Method 3 can be regarded as a combination of steepest descent and Newton's method for unconstrained optimization. Method 4 can be viewed as an extention of the Levenberg-Marquardt method for unconstrained optimization.

Theorem 2. The sequence {u k } generated by each method of methods 1-4 for (LPE) satisfies

(18) u k+1 -u * 2 G ≤ u k -u * 2 G -ρ(u k ) e(u k ) 2 -2ρ(u k ) • (u k -u * )M (u k -u * ) ∀u * ∈ Ω * where G = ⎧ ⎪ ⎨ ⎪ ⎩ I + M in Method 1 (I + M )M in Method 2 M in Method 3 (I + M T )(I + M ) in Method 4 Proof. First, u k+1 -u * 2 G = (u k -u * ) -ρ(u k )Qe(u k ) 2 G = u k -u * 2 G -2ρ(u k )(u k -u * )GQe(u k ) + ρ 2 (u k )e(u k ) T Q T GQe(u k ).
Note that in all cases

GQ = I + M T and ρ(u) • e(u) T Q T GQe(u) = e(u) 2 .
Using [START_REF] Lemke | Equilibrium points of bimatrix games[END_REF] the theorem is proved.

The sequence {u k } generated by these methods does not necessarily lie in Ω. In general, (because G / = I and the projection is an orthogonal projection with respect to the Euclidean norm), we can not prove that

{u k } satisfies P Ω [u k+1 ] -u * G ≤ u k -u * G even if u k ∈ Ω.
Note that in all of these methods, the steplength ρ is bounded below. Therefore, there is a c > 0, so that the sequence {u k } generated by each of these methods satisfies

(19) u k+1 -u * 2 G ≤ u k -u * 2 G -c e(u k ) 2 ∀u * ∈ Ω * .
Due to [START_REF] Stoer | On the numerical solution of constrained least squares problems[END_REF] and the fact that each iteration performs a projection onto Ω (in the computation of e(u)), we call these methods the projection and contraction methods (PC methods).

Convergence

The PC methods in this paper generate an infinite sequence {u k }, which is not necessarily contained in the feasible set Ω, but, will be asymptotically feasible as e(u k ) → 0, and, in fact converges to a solution of (LPE).

Theorem 3. If the sequence {u k } satisfies [START_REF] Stoer | On the numerical solution of constrained least squares problems[END_REF], then it converges to a solution point u * .

Proof. See [START_REF] He | A new method for a class of linear variational inequalities[END_REF], Theorem 3.

In [START_REF] He | A new method for a class of linear variational inequalities[END_REF] we have also proved that an inequality of the form [START_REF] Stoer | On the numerical solution of constrained least squares problems[END_REF] implies the global linear convergence in the case that Ω is an orthant. Theorem 4. If the sequence {u k } satisfies [START_REF] Stoer | On the numerical solution of constrained least squares problems[END_REF] and

Ω = {u | u ≥ 0}, then {u k } converges to a solution point u * ∈ Ω * globally linearly.
Proof. See [START_REF] He | A new method for a class of linear variational inequalities[END_REF], Theorem 4.

The iterative scheme of the fundamental projection method (see [START_REF] Ciarlet | Introduction to matrix numerical analysis and optimization[END_REF]) (20)

u k+1 = P Ω [u k -(M u k + q)],
belongs to the class of steepest descent methods. In the case that M is positive semidefinite, we denote the largest and the smallest eigenvalue of the matrix M by λ max (M ) and λ min (M ), respectively. If 0 < δ ≤ λ min (M ) ≤ λ max (M ) ≤ 2δ, then the sequence {u k } generated by the fundemental projection method (20) satisfies ( 21)

e(u k+1 ) = u k+1 -P Ω [u k+1 -(M u k+1 + q)] = P Ω [u k -(M u k + q)] -P Ω [u k+1 -(M u k+1 + q)] ≤ (I -M )(u k -u k+1 ) = (I -M )e(u k ) ≤ (1 -δ) e(u k ) and (22)
u k+1 -u * = P Ω [u k -(M u k + q)] -u * = P Ω [u k -(M u k + q)] -P Ω [u * -(M u * + q)] ≤ (I -M )(u k -u * ) ≤ (1 -δ) u k -u * .
In the following we prove that the PC methods in this paper have similar properties. Under the assumptions that M is positive definite and M ≤ 2, the seqeunce { e(u k ) } generated by Method 1 or Method 4 is monotonically decreasing. More precisely, we have the following theorem:

Theorem 5. Let M be positive semidefinite and symmetric. If λ max (M ) ≤ 2, then the sequence {e(u k )} generated by Method 1 or Method 4 for (LPE) satisfies

(23) e(u k+1 ) ≤ e(u k ) . Moreover, if δ ≤ λ min (M ) ≤ λ max (M ) ≤ 2 -δ for some δ > 0, then the sequence {e(u k )} satisfies (24) e(u k+1 ) ≤ (1 - δ 3 ) e(u k ) .
Proof. First, in Method 1, under the assumptions we have 1 3 ≤ ρ(u k ) ≤ 1. By using

u k+1 = u k -ρ(u k )e(u k ),
we get

e(u k+1 ) = u k -ρ(u k )e(u k ) -P Ω [u k+1 -(M u k+1 + q)] = (1 -ρ(u k ))e(u k ) + u k -e(u k ) -P Ω [u k+1 -(M u k+1 + q)] ≤ (1 -ρ(u k ))e(u k ) + P Ω [u k -(M u k + q)] -P Ω [u k+1 -(M u k+1 + q)] ≤ (1 -ρ(u k )) e(u k ) + (I -M )(u k -u k+1 ) ≤ (1 -ρ(u k )) e(u k ) + ρ(u k ) (I -M ) • e(u k ) . Since 0 ≤ λ min (M ) ≤ λ max (M ) ≤ 2, it follows that I -M ≤ 1 and e(u k+1 ) ≤ e(u k ) . Moreover, if δ ≤ λ min (M ) ≤ λ max (M ) ≤ 2 -δ, then I -M ≤ 1 -δ, and it follows that e(u k+1 ) ≤ (1 -ρ(u k )) e(u k ) + (1 -δ)ρ(u k ) e(u k ) ≤ (1 - δ 3 ) e(u k ) .
In Method 4, using

u k+1 = u k -(I + M ) -1 e(u k ) we get e(u k+1 ) = u k -(I + M ) -1 e(u k ) -P Ω [u k+1 -(M u k+1 + q)] = e(u k ) -(I + M ) -1 e(u k ) + u k -e(u k ) -P Ω [u k+1 -(M u k+1 + q)] ≤ (I + M ) -1 M e(u k ) + P Ω [u k -(M u k + q)] -P Ω [u k+1 -(M u k+1 + q)] ≤ (I + M ) -1 M e(u k ) + (I -M )(u k -u k+1 ) = (I + M ) -1 M e(u k ) + (I -M )(I + M ) -1 e(u k ) .
Let M = U T ΣU be the Schur normal form of M with Σ = diag(σ 1 , . . . , σ n ). Then

e(u k+1 ) ≤ (I + M ) -1 M e(u k ) + (I -M )(I + M ) -1 e(u k ) = (I + Σ) -1 ΣU e(u k ) + (I -Σ)(I + Σ) -1 U e(u k ) .
Since 0 ≤ σ i ≤ 2, it follows that

σ i 1 + σ i + |1 -σ i | 1 + σ i ≤ 1 and e(u k+1 ) ≤ (I + Σ) -1 ΣU e(u k ) + (I -Σ)(I + Σ) -1 U e(u k ) ≤ U e(u k ) = e(u k ) . Moreover, if δ ≤ σ i ≤ 2 -δ, it is straightforward to prove that σ i 1 + σ i + |1 -σ i | 1 + σ i ≤ 1 - δ 3 and it follows that e(u k+1 ) ≤ (1 - δ 3 ) e(u k ) .
The following theorem contrasts the convergence proofs for the fundamental projection methods, in that no advance knowledge of the largest and smallest eigenvalue of M is required. Theorem 6. Let M be positive definite, then all four PC methods are globally linearly convergent.

Proof. From (18) (Theorem 2) we have

u k+1 -u * 2 G ≤ u k -u * 2 G -2c(u k -u * ) T M (u k -u * ). Since M is positive definite, there exists a τ > 0, such that (u -u * ) T M (u -u * ) ≥ τ u -u * 2 G . It follows that u k+1 -u * 2 G ≤ (1 -2cτ ) u k -u * 2
G , which implies that {u k } converges to u * globally and linearly.

Extensions and conclusions

Let α > 0 be a constant. It is easy to see that problem (1) is equivalent to the following problem (25)

(LPE α ) u = P Ω [u -α(M u + q)].
Therefore, instead of taking M and q, we can use αM and αq in our PC methods. In addition, for some parameter γ, 0 < γ < 2, with the same direction d(u) and its relevant steplength ρ(u), the iterative scheme (26) u k+1 = u kγρ(u k )d(u k ) produces a sequence {u k }, which satisfies

(27) u k+1 -u * 2 G ≤ u k -u * 2 G -γ(2 -γ)ρ(u k ) e(u k ) 2
and thus also converges to a solution point u * . A close look at the inequalities used in the proof of Theorem 1 shows that the best choice of γ should be ≥ 1.

The search directions of the presented methods are the extensions of those in unconstrained optimization. If Ω = R n , then the search direction in [7] is

g(u) = (M T + I)(M u + q).
In this special case, the search directions in our new methods are d 1 (u) = (M u + q), (in Method 1) d 2 (u) = M -1 (M u + q), (in Method 2) d 3 (u) = (I + M -1 )(M u + q), (in Method 3) d 4 (u) = (I + M ) -1 (M u + q), (in Method 4)

respectively. Although we have only proved linear convergence, by comparing analogous direction for unconstrained optimization we are convinced that the directions d(u) developed in this paper are better than the direction g(u) in the original PC methods [START_REF] He | A projection and contraction method for a class of linear complementarity problem and its application in convex quadratic programming[END_REF][START_REF] He | A globally linearly convergent projection and contraction method for a class of linear complementarity problems[END_REF]. We believe that the use of Newton-like directions will lead to a substantial improvement in computational efficiency. Developing these methods for nonlinear problems is a topic of further research .

Acknowledgement. The author would like to thank the anonymous referees for several suggestions for improvement.