Puneet Kumar

Mathias Perrollaz

Stéphanie Lefèvre

Christian Laugier

Learning-Based Approach for Online Lane Change Intention Prediction

Predicting driver behavior is a key component for Advanced Driver Assistance Systems (ADAS). In this paper, a novel approach based on Support Vector Machine and Bayesian filtering is proposed for online lane change intention prediction. The approach uses the multiclass probabilistic outputs of the Support Vector Machine as an input to the Bayesian filter, and the output of the Bayesian filter is used for the final prediction of lane changes. A lane tracker integrated in a passenger vehicle is used for real-world data collection for the purpose of training and testing. Data from different drivers on different highways were used to evaluate the robustness of the approach. The results demonstrate that the proposed approach is able to predict driver intention to change lanes on average 1.3 seconds in advance, with a maximum prediction horizon of 3.29 seconds.

I. INTRODUCTION

Studies of the causes of road accidents show that 57% of them are solely due to driver factors [START_REF] Lum | Interactive highway safety design model: accident predictive module[END_REF]. Advanced Driver Assistance Systems (ADAS) can help drivers to understand traffic situations better, take actions to make driving more comfortable, improve the traffic flow, limit the energy consumption, avoid accidents, or mitigate their consequences. For example, predicting in advance other drivers's intentions at an intersection or predicting whether a vehicle on the highway will change lanes can help a driver to understand the situation better and to avoid accidents.

This work focuses on lane change intention prediction on highways. The goal is to predict whether the ego vehicle will make a lane change and in which direction, based on sensor data. Mathematically, a vehicle trajectory can be defined by the deterministic function:

Φ : Φ(t) = x(t), y(t), dx dt (t), dy dt (t), d 2 x d 2 t (t), d 2 y d 2 t (t) (1)
where the outputs of the function Φ represent the longitudinal position, lateral position, longitudinal velocity, lateral velocity, longitudinal acceleration and lateral acceleration of the vehicle respectively [START_REF] Boubezoul | Vehicle trajectories classification using support vectors machines for failure trajectory prediction[END_REF]. The trajectory segment in the interval [T 1 T 2 ] is a lane change trajectory segment iff it intersects with the curve formed by the road markings. The time where the lane change occurs is defined as the time at which the trajectory segment and the lane markings intersect, as depicted in Fig. 1. In this context the problem we tackle P. Kumar is currently with École Centrale de Paris, Center for Visual Computing, and GALEN team of Inria Saclay. The work presented in this paper was done when he was with Inria Grenoble Rhone-Alpes.

M. Perrollaz, S. Lefèvre, and C. Laugier are with Inria Grenoble Rhone-Alpes, 655 av. de l'Europe -Montbonnot, 38334 Saint Ismier Cedex, France. {puneet.kumar, mathias.perrollaz, stephanie.lefevre, christian.laugier} at inria.fr In this paper a novel approach based on the combination of multiclass Support Vector Machine (SVM) and Bayesian filtering is proposed. The algorithm was tested on a real dataset collected using a passenger vehicle. The performance of the approach was evaluated by looking at the prediction horizon, the rate of false alarms, and the rate of missed detections.

The remaining of this paper is organized as follows. Section 2 reviews related work, Section 3 discusses the different aspects and constraints of the problem and then describes our approach to solve the lane change prediction problem within the discussed constraints. Section 4 presents the results. Section 5 concludes and outlines future work.

II. RELATED WORK

The term intention prediction can have several interpretations. Different terminologies such as behavior prediction/recognition/identification, situation assessment/prediction, intention prediction/estimation have been used in the literature to qualify the same type of problems. Examples of such problems are: predicting whether a driver will change lanes, predicting whether a vehicle will stop at the red traffic light etc.

In [START_REF] Mandalia | Using support vector machines for lane change detection[END_REF] the lane change prediction problem is addressed using SVMs with feature vectors consisting of the variances of the features (speed, steering angle etc.). In [START_REF] Mccall | Lane change intent analysis using robust operators and sparse bayesian learning[END_REF] a sparse Bayesian classifier is used for lane change intent analysis. The idea is to use time series data describing the vehicle's surrounding, the driver's head motion, and the vehicle's internal state to create a feature vector for classification. The results show that the inclusion of the driver's state (head motion) in the feature vector increases the prediction horizon. The extension of this work proposed in [START_REF] Morris | Lane change intent prediction for driver assistance: on-road design and evaluation[END_REF] uses Relevance Vector Machine (RVM), a Bayesian extension of SVM. It relies on information from ACC (Adaptive Cruise Control) radar, LDW (Lane Departure Warning) camera, SWA (Side Warning Assist) radars, and head tracking camera. This approach is able to predict lane changes up to 3 seconds in advance. Another work [START_REF] Aoude | Behavior classification algorithms at intersections and validation using naturalistic data[END_REF] deals with intentions at road intersections and considers two types of on-road agents: compliant and violating. Two approaches were compared, the first one using an SVM-based binary classification along with Bayesian filtering and the second one using Hidden Markov Models (HMMs). The results showed that the SVM-BF based approach performed better. In [START_REF] Oliver | Graphical models for driver behavior recognition in a smartcar[END_REF] Coupled HMMs were trained to create models of seven different driving maneuvers. The results showed that maneuvers could be predicted on average 1 second before they actually started. In [START_REF] Meyer-Delius | Probabilistic situation recognition for vehicular traffic scenarios[END_REF] a probabilistic approach using HMMs is used for situation modeling and recognition. A situation is defined as a distribution over sequences of states having a meaningful interpretation and a HMM is used to characterize each situation. The approach was tested on real data for highway driving with three situations: passing, aborted passing, and following. Similarly to this approach, in [START_REF] Tay | Analysis of dynamic scenes: application to driving assistance[END_REF] a Hierarchical HMM was used for behavior recognition. The idea is to model behaviors in two layers. High level behaviors, such as: go straight, turn left, turn right, and overtake are treated as the hidden states of an HMM in the upper layer. For each high-level behavior there exists a HMM at the lower layer representing the sequence of transitions of the corresponding behavior. In [START_REF] Hülnhagen | Maneuver recognition using probabilistic finite-state machines and fuzzy logic[END_REF] a Probabilistic Finite State Machine (PFSM) and fuzzy logic are used for maneuver recognition. The input variables (velocity, steering angle etc.) are fuzzified in order to estimate the basic elements (braking, halt, start etc.) constituting maneuvers. Finally a Bayesian filter is used to find the probability distribution of the basic elements of the PFSM for maneuver recognition.

III. PROPOSED APPROACH

Our approach uses lane information, speed, and steering angle for lane change intention prediction in highway scenarios. We formulate the lane change intention prediction problem as a multiclass classification problem with the three following families/classes of trajectories: left lane change, right lane change, and no lane change. Other situations like multiple lane changes can be seen as combinations of these three classes. A kernel-based large margin classifier known as Support Vector Machine (SVM) is used to address this problem. Major motivations behind the use of SVMs are: (a) a driver's state may lie in a high dimensional feature space [START_REF] Wipf | Driver intent inference annual report[END_REF] and a kernel maps the input data from a low dimensional space into a high dimensional space converting a nonlinear classification problem at low dimension into a linear classification problem at high dimension, (b) SVM is a maximum margin classifier, therefore it is expected to classify similar trajectories belonging to different classes reliably, (c) the objective function of an SVM is convex, therefore the solution is a global optimum, and (d) several studies (e.g. [START_REF] Aoude | Threat assessment for safe navigation in environments with uncertainty in predictability[END_REF]) have shown that SVMs give promising results when applied to intention prediction problems for ADAS.

A Bayesian Filter (BF) is used on top of the multiclass classifier in order to improve the reliability of the predictions. It is expected that the smoothing introduced by the filter will reduce the rate of false alarms and missed detections. A generalized Bradley-Terry model [START_REF] Huang | Generalized bradley-terry models and multi-class probability estimates[END_REF] is used to obtain probabilistic outputs from the SVM. This probabilistic output of the SVM is fed to the BF providing a filtered and smooth output. The state transition matrix of the BF is learned from real data. The output of the BF is used for final lane change intention prediction.

Details about the proposed approach are provided in the rest of this section: (a) lane tracker, (b) multiclass probabilistic estimates using SVM, and (c) Bayesian Filter.

A. Lane tracker

To make the overall approach independent of localization devices (e.g. GPS) and digital maps (GIS), a vision and IMU based lane tracker is used to capture information about the position of the ego vehicle with respect to the road (see device in Fig. 3). In general, the lane tracking problem can be seen as a state estimation problem with the state vector consisting of the lane parameters defining the road geometry [START_REF] Goldbeck | Lane detection and tracking by video sensors[END_REF]. Here we use a particle filter and follow the following steps:

1) Define parametric equations describing the lane position and geometry with respect to the vehicle. The lane is assumed to have a shape similar to a clothoid curve. 2) Prediction step of the particle filter: the state vector consisting of the parameters of the parametric equations is recursively estimated using the past state vector and odometry information. 3) Image capturing and ridge detection: a ridge filter provides the low level features used in the update stage of the particle filter. 4) Update step of the particle filter: the weight of each particle is adapted to reflect how well it matches the ridge features.

B. Multiclass probabilistic estimates using SVM

Support Vector Machine (SVM) is a supervised learning algorithm that uses the concept of margin maximization for the purpose of classication. The basic idea is to find a hyperplane that maximizes the separation between the datapoints of different classes. It uses the concept of kernels to project data from low dimensions to higher dimensions thus converting a nonlinear problem in low dimensions into a linear problem in higher dimensions. Given a dataset S = {(x 1 , y 1 ), ..., (x m , y m )} ∈ (χ × γ) m , with χ = ℜ n and γ = {-1, 1}, the convex objective function for the binary SVM is defined as [START_REF] Cristianini | An introduction to support vector machines and other kernel-based learning methods[END_REF], [START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF]: One simple approach for the multiclass extension of the above defined binary SVM is to learn k different ws for k classes using a one-vs-all approach.

w * = min w,ζ 1 2 w T w + C m m ∑ i=1 ζ i s.t. y i (w T x i + b) ≥ 1 -ζ i , ∀i (2)

C. Feature Selection

The discriminative capabilities of any classifier greatly depends on the selection of the feature vectors x i . We use feature vectors consisting of the following meta-features:

• the lateral position of the vehicle w.r.t. a lane (l)

• the steering angle of the vehicle w.r.t. the road (φ )

• the first derivative of l • the first derivative of φ These features were selected because their patterns are remarkably different for each of the three classes we are interested in. Since a lane change is a continuous process and can last for several seconds, the feature vector x i must contain information concatenated within a time span in order to try to capture the continuity. To do so, a sliding window approach is used (see Fig. 2). A window is selected around the time when the vehicle reaches the lane markings for a lane change, all the meta-features within that window are labeled as the ground-truth data-points representing the corresponding lane change. This window size should be large enough to capture a lane change process completely. Within this window a sub-window is created and different meta-features at different times within this sub-window are concatenated to form the final feature vector of fixed length. This subwindow is moved within the window in order to generate different feature vectors representing lane changes. Mathematically, if (t1 + t2) represents the window size, f represents the sub-window size, f ps represents the data frames collected per second, then the feature vector at time i can be represented as:

x i = [l i , li , Φ i , Φi ] ∈ ℜ n ( 3
)
where

t = 1 f ps n = 4 × f l i = [l i-f ×t , l i-( f -1)×t , . . . , l i-t , l i ] li = [ li-f ×t , li-( f -1)×t , . . . , li ] Φ i = [φ i-f ×t , φ i-( f -1)×t , . . . , φ i ] Φi = [ φi-f ×t , φi-( f -1)×t , . . . , φi ]
We use a radial basis function kernel, K (x i , x j ) = exp(-g x i -x j 2 ), with g as the kernel parameter and d > 0. The approach proposed by [START_REF] Huang | Generalized bradley-terry models and multi-class probability estimates[END_REF] is used to get the probabilistic outputs from the multiclass SVM. This approach extends the Bradley-Terry model for paired individual comparisons into paired team comparisons.

D. Bayesian Filter

The Bayesian filtering algorithm is the recursive form of Bayes rule which gives an inverse relationship between posterior, prior and likelihood. For lane change intention prediction, the equation to compute the probability of a maneuver M t can be written as:

P(M t |Z 0:t ) ∝ P(Z t |M t )× [ ∑ M t-1 =L,R,N P(M t |M t-1 ) × P(M t-1 |Z 0:t-1 )] (4)
where Z t corresponds to the observations at time t and L, R, N represent left, right and no lane changes respectively. In our approach the likelihood term in the above equation, P(Z t |M t ), is taken as the probabilistic output of the SVM. The state transition probabilities P(M t |M t-1 ) are learned offline using training data. For all the three maneuvers we get a 3 × 3 state transition matrix with nine unknown state transition probabilities. The corresponding state transition diagram for all the three states or classes L, R, N is shown in Fig. 2. Since in our case all the maneuvers M t in our training data can be labeled, the state transition matrix is simply the normalization of the co-occurrences (counts):

S ML (i, j) = S(i, j) ∑ k S(i, k)
where, S(i, j) is the number of transitions from maneuver i to maneuver j in the training data and S ML (i, j) is the maximum likelihood probability for i → j transition .

The posterior P(M t |Z 0:t ) is used for the final lane change intention prediction.

IV. RESULTS

A dataset with 139 lane changes was collected with two different drivers on a highway near Grenoble, France, using our Lexus LS600h experimental platform (see Fig. 3). The vehicle is equipped with a TYZX stereo camera with 22 cm baseline, 62 • field of view, 512 × 320 pixels resolution and focal length of 410 pixels. The training dataset consists of 22 left lane changes, 24 right lane changes, and 24 trajectories with no lane change. The testing dataset consists of a total of 69 lane changes (left + right). The parameters of the SVM were set to g = 0.0625, C = 8, t1 = 2, t2 = 2, f = 32, fps = 32, and final classification rate = 5 classifications per second. The final classification rate is lower than the fps because 32 classifications per second is more than what is necessary for a real-time application. On the other hand, for training we used all the data acquired at 32 fps because data at higher fps better approximates the lane change event. The reasons for choosing the particular values for the other parameters are discussed in subsequent sections.

A. Evaluation metrics

The following terminologies and evaluation metrics are used:

• Ground truth: the moment at which the vehicle reaches the lane markings.

B. Classification Results

Fig. 4 shows the ground truth for 5 lane changes (3 left + 2 right) belonging to the testing data. The final predicted classes for the same data using SVM alone and SVM + BF are shown in Fig. 5. We observe many false alarms if the SVM alone is used for the lane change intention prediction. Tab. I gives the performance results in terms of the average precision, recall and prediction time for the SVM and SVM+BF approaches. We can see from the table that the precision is improved from 0.2857 to 0.7154 by adding the Bayesian Filter. The average prediction time is decreased by 0.022 seconds, which is acceptable with such a remarkable improvement on the precision. Fig. 6 displays the histogram of prediction time for 54 lane changes. Most of the lane changes are predicted almost 1.3 seconds before the ego vehicle crosses the painted line and that the maximum prediction horizon reaches 3.29 seconds.

C. Finding the best parameters

As mentioned earlier there are mainly five parameters involved during the implementation of the proposed approach (SVM+BF) and each parameter has a strong effect on the performance. These parameters are: g, C, t1, t2, and f. Simultaneously finding the best parameters is a five dimensional grid search problem and since the values of each parameter may vary in a large range it is computationally very expensive to do. Instead a simplified strategy is used:

1) Fix t1, t2, f and find the best C and g using grid search with exponentially growing sequences of C and g [START_REF] Hsu | A practical guide to support vector classification[END_REF], for example C = (2 -3 , 2 -1 , . . . , 2 3 ) and g = (2 -4 , 2 -2 , . . . , 2 2 ). While fixing t1, t2 and f it was made sure that these parameters gave good results when tested on different scenarios. 2) Fix C and g to their best values, fix t1 and t2, and find the best value for f . 3) Fix C, g, and f to their best values and then iterate over a range of t1 and t2 to find their best values. Fig. 6 shows performances for different values of the parameters. The plots show that choosing bad parameters can lead to situations with very poor precision or recall.

D. Robustness evaluation

To test the effect of different drivers on the performance of the proposed approach, the following cases were considered with our two drivers x and y:

Fig. 1 .

 1 Fig. 1. Lane change scenario

Fig. 2 .

 2 Fig. 2. Left: Sliding window approach for ground-truth feature vector generation. Right: State Space Diagram for different maneuvers

Fig. 3 .

 3 Fig. 3. The Lexus experimental platform used for real data collection

  • Classes: Class 1 = right lane change Class 2 = no lane change Class 3 = left lane change • Prediction point: the time at which a lane change is predicted. • Prediction time: the time difference between the ground truth and the corresponding prediction point. In the evaluation, the prediction point nearest to the ground truth is always chosen (worst case). If the prediction point occurs later than the ground truth then this is considered as false prediction. • Precision, recall and F1-Score: the precision is the fraction of trajectories where the classifier correctly declared class m out of all instances where the classifier declared class m. The recall is the fraction of events where the classifier correctly declared class m out of all of the cases where the true class was m. The F1-Score is the harmonic mean of the precision and the recall.

Fig. 4 .

 4 Fig. 4. Ground truth testing data with 5 lane changes (3 left + 2 right)

• Case 1 :

 1 Training with x and testing with y. • Case 2: Training with y and testing with x. • Cases 3(a) and 3(b): Training with x and testing with different driving scenarios with x. Tab. II shows the performance for each case. The recall is always 1. The F-1 score, precision and prediction time are almost the same, which shows the robustness of the system.

Fig. 5 .

 5 Fig. 5. Comparison of lane change prediction using SVM and SVM+BF. (Top Left + Top Right): Probability distribution for all the three classes using SVM, and lane change prediction based on this distribution. (Bottom Left + Bottom Right): Probability distribution for all the three classes using SVM+BF, and lane change prediction based on this distribution.

Fig. 6 .

 6 Fig. 6. (Top Left): Histogram of prediction time. (Top Right): Plot of precision, recall and average prediction time for left and right lane changes for the selection of the best (C,g). (Bottom Left): average precision, recall and prediction time for finding the best f. (Bottom Right): average precision, recall and prediction time for finding the best (t1, t2). Time is expressed in seconds.

TABLE I COMPARISON

 I OF SVM AND SVM+BF (TESTED ON 69 LANE CHANGES)

	Cases	Precision	Recall Avg prediction time (sec)
	SVM	0.2857	1	1.2947
	SVM + BF	0.7154	1	1.2718

V. CONCLUSIONS

This paper proposed a solution to lane change intention prediction based on a combination of a multiclass SVM classifier and Bayesian filtering. When tested on real data, the algorithm is able to predict the lane changes accurately (recall = 1) on average 1.3 seconds before they occur and to differentiate between left lane changes and right lane changes. Future work will focus on the reduction of the number of false alarms. In our current system most of the false alarms are caused by inaccuracies of the lane tracker, which sometimes fluctuates and jumps. Therefore in the future the lane tracker should be improved. Further improvements could be obtained by incorporating additional information to our prediction framework, such as the distance to the vehicle in front or the speed difference with the vehicle in front. Indeed, these are useful cues to estimate whether the driver intends to overtake the vehicle in front. Also an extended evaluation with a larger dataset and several drivers is to be done.


Mounia Haddouni
email: mounia.haddouni@inria.fr

Vincent Acary
email: vincent.acary@inria.fr

Jean-Daniel Beley
email: jean-daniel.beley@ansys.com

Comparison of index-3, index-2 and index-1 DAE solvers for nonsmooth multibody systems with unilateral and bilateral constraints

Keywords: nonsmooth mechanics, unilateral contact, impacts, event-driven schemes, constraints, DAE, numerical solvers, LCP, NSCD

This work is devoted to the study of some numerical solvers for the numerical time integration of nonsmooth multibody systems with unilateral or bilateral constraints. In the framework of event-detecting time-stepping schemes (a.k.a. eventdriven schemes), an index-3 Differential Algebraic Equation (DAE) has to be solved between two events with a constant number of bilateral constraints. In this paper, we compare several solvers for index-1, index-2 and index-3 DAEs in the specific context of unilateral contact and impact. These solvers will be compared in terms of drift of the constraints which is an important feature when we have to update the index sets of the active unilateral constraints. Their implementation and efficiency with fixed and adaptive time-step strategies will also be reviewed. Finally, we discuss the way they handle cases when the constraints are not sufficiently differentiable, which is of utmost importance for practical applications.

Introduction

The simulation of multi-body systems composed of rigid bodies or flexible bodies, such as the one presented in Fig. 1, with joints, unilateral contact and Coulomb's friction requires the development of sophisticated numerical schemes (see [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Studer | Numerics of Unilateral Contacts and Friction. -Modeling and Numerical Time Integration in Non-Smooth Dynamics[END_REF] for a review). In the framework of event-detecting time-stepping schemes (a.k.a. event-driven schemes), an index-3 Differential Algebraic Equation (DAE) has to be solved between two events with a constant number of bilateral constraints [START_REF] Eich-Soellner | Numerical methods in multibody dynamics[END_REF][START_REF] Schwerin | Multibody System Simulation. Numerical Methods, Algorithms and Software[END_REF][START_REF]Simulation Techniques in Applied Dynamics[END_REF]. In this paper, we compare several solvers for index-1, index-2 and index-3 DAEs in the specific context of unilateral contact and impact. These solvers will be compared in terms of drift of the constraints, implementation and the way they handle cases when the constraints are not sufficiently differentiable.

Figure 1: Some illustrations of systems with nonsmooth constraints. A valve system with a zoom on the cylindrical cam and a standard cam-follower system.

The dynamics of a rigid multibody system with m frictionless contact points can be written, using the Lagrangian formalism, as

     q(t) = v(t)
(1a) M (q(t)) v(t) = F (t, q(t), v(t)) + G (q(t))λ(t) (1b) 0 ≤ λ(t) ⊥ g(q(t)) ≥ 0,

where q(t) ∈ R n is the generalized coordinates vector, v(t) ∈ R n is the generalized velocities vector, M (q(t)) ∈ R n×n is the symmetric definite positive matrix of inertia, F (t, q(t), v(t)) comprises the external applied loads, the non linear interactions between bodies and the non linear inertial terms, g(q(t)) ∈ R m is the vector of constraints, G = ∇ q q(q) is the Jacobian of the constraints and λ(t) ∈ R m is the Lagrange multiplier which is related to the force at the contact point. The complementarity condition (1c) models the unilateral contact, also termed as the Signorini condition at the position level.

We suppose that impacts occur in infinitely short periods so that the displacements of the bodies during the contact period can be negligible and we use the Newton impact law with a coefficient of restitution e. If g(q(t i )) = 0, that is the contact is closed (we say also active) at time t i , we compute the impulse p i and the post-impact velocity v + (t i ) by solving the impact equations

       M (q(t i ))(v + (t i ) -v -(t i )) = G (q(t i ))p i U + N (t i ) = G(q(t i ))v + (t i ) U - N (t i ) = G(q(t i ))v -(t i ) 0 ≤ U + N (t i ) + e U - N (t i ) ⊥ p i ≥ 0, (2)
where U N (t) = ġ(q(t)) = G(q(t)) q(t) is the relative normal velocity at contact. The complementarity condition in (2) describes the Signorini condition written at the velocity level and augmented of the impact law.

When an event-driven strategy is chosen for the numerical time integration, a solver for differential algebraic equations has to be used between two events for solving the smooth system with equality constraints g α (q) = 0, α ∈ I 2 . The index set I 2 belongs to the family of index sets generally introduced in order to characterize the state of the contacts:

• the index set I 0 of all possible constraints in the system I 0 = {1, ..., m},

• the index set I 1 of contacts activated in position:

I 1 = {i ∈ I 0 |g i (q(t)) = 0}, • the index set I 2 of all closed contacts I 2 = {i ∈ I 1 | ġi (q(t)) = 0}.
From the numerical point of view and with the assumption that the constraints are smooth enough, the following smooth dynamics,

     q(t) = v(t) (3a) M (q(t)) v(t) = F (t, q(t), v(t)) + G (t, q(t))λ (3b) g α (q(t)) = 0, λ α ≥ 0, α ∈ I 2 . ( 3c
)
is integrated with a DAE solver, during periods of time over which the index set I 2 is constant. An event detection is performed to find the time where new contacts occur or some contacts need to be released. When an event is detected, the impact equations (2) are solved and the state of the dynamical system is updated. A complete algorithm for the event-driven schemes is detailed in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF].

It is well known that index-3 DAEs are difficult to numerically handle [START_REF] Eich-Soellner | Numerical methods in multibody dynamics[END_REF][START_REF] Brüls | Convergence of the generalized-α scheme for constrained mechanical systems[END_REF]. Therefore the dynamics is usually integrated as an Ordinary Differential Equation (ODE) by reducing the original index of the system from 3 to 1. It amounts to solving the problem at the acceleration level by differentiating twice the constraints. Roughly, index reduction consists in differentiating the constraint as many times as necessary to get a set of equations that may be solved using methods for lower index problems. Hence, if the constraint g is differentiated once with respect to time, one obtains an index-2 DAE, we solve in this case (3a) and (3b) together with

G(t, q(t))v(t) = 0. ( 4
)
If g is twice differentiated, one gets an index-1 DAE composed of (3a) and (3b) with

G(t, q(t)) v(t) + G qq (v, v) = 0. (5)
Numerically, the index reduction leads to the phenomenon of drift of the constraints, at different levels depending on the index. Indeed, as opposed to the continuous time case, enforcing (5) in discrete time does not imply that (3c) and (4) are satisfied.

2 Presentation of the chosen time-integration schemes 2.1 An index-3 DAE solver: the generalized-α-scheme

This scheme dedicated to the integration of index-3 DAEs (3) has been studied by O. Brüls and M. Arnold [START_REF] Brüls | Convergence of the generalized-α scheme for constrained mechanical systems[END_REF]. It is an adaptation of the standard α-scheme for the resolution of index-1 DAEs in the context of computational Mechanics of solids. Let us introduce the acceleration-like variable a defined by the recurrence relation

(1 -α m )a n+1 + α m a n = (1 -α f )q n+1 + α f qn . (6)
The discretization of equation ( 3) is given by

q n+1 = q n + h qn + h 2 ( 1 2 -β)a n + h 2 βa n+1 qn+1 = qn + h(1 -γ) + hγa n+1 , (7)
where the constants α, α m , β and γ should be chosen suitably so that the scheme is stable. In [START_REF] Brüls | Convergence of the generalized-α scheme for constrained mechanical systems[END_REF], it is said that the algorithm is unconditionally stable if the coefficients are chosen such that for ρ ∞ < 1

α m = 2ρ ∞ -1 ρ ∞ + 1 < 1 2 , α f = ρ ∞ ρ ∞ + 1 < 1 2 , γ = 1 2 + α f -α m > 1 2 , β = 1 4 (γ + 1 2 ) 2 . ( 8
)
The scheme is based on a prediction step and a correction step where some Newton iterations are performed in order to reduce the dynamical and the constraint residuals defined by

R q = M (q, t)q -F (q, q, t) -G (q)λ R λ = g(q). (9)
The Newton iterations amount to solving of the following linear system β M (q, t) + γ C t (q, q, t)

+ K t (q, q, t) -G T (q) G(q) 0 ∆q ∆λ = R q R λ , (10)
where

β = 1-αm h 2 β(1-α f ) , γ = γ hβ , M is the inertia matrix, K t = ∂(M q-F +G λ)
∂q is the stiffness matrix, and

C t = -∂F ∂ q
is the damping matrix. Note that this algorithm aims at maintaining the constraints at the position level, but it can also be reformulated to write the constraints at the velocity level or at the acceleration level [START_REF] Ch | Solving constrained mechanical systems by the family of Newmark and α-methods[END_REF].

For the control of time step size, we compute an estimation of the solution y = [q n+1 , qn+1 ] using a 1 st order Newmark scheme with γ > 1/2, and we define the truncation error as

err = y α-scheme -y Newmark . (11)
The optimal step size is then given with

h opt = safe tol h err 1/p h, ( 12
)
where safe is a safety factor, tol the user required tolerance and p = 2 is the order of the α-scheme.

An index-2 DAE solver: the Half-Explicit Method of order 5 (HEM5)

The HEM5 solver is based on a half-explicit 5 order Runge-Kutta(RK) method with 8 stages, described in [START_REF] Brasey | Hem5 user's guide[END_REF]. Half-Explicit RK methods have been deeply studied by E. Hairer in [START_REF] Hairer | Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems[END_REF]. The HEM5 solver is suited to the integration of DAEs, reduced to the index-2 formulation. The numerical approximation of the equations of motion is given by

   M (Q i , t n + c i h) Vi = F (t n + c i h, Q i , V i ) + G (t n + c i h, Q i )Λ i Qi = V i G(t n + c i h, Q i )V i = 0, (13)
where the stages are evaluated as follows :

Q i = q n + h j<i a ij Qj and V i = v n + h j<i a ij Vj .
At each stage, the estimations of position and velocity Q i and V i are explicitly computed, while Λ i and Vi are implicitly obtained by solving the linear system

M (Q i , t n + c i h) -G T (Q i , t n + c i h) G(Q i+1 , t n + c i+1 h) 0 Vi Λ i = F (Q i , V i , t n + c i h) r i , (14)
where

r i = - G(Q i+1 ) ha i+1,i (v n + h i-1 j=1
a i+1,j Vj ). At the end of the time step, the numerical solution is given by q n+1 =

Q 9 = q n + h 8 i=1 b i Qi and v n+1 = V 9 = v n + h 8 i=1 b i Vi , with b i = a 9i
. The computation of coefficients c i and a ij is detailed in [START_REF] Brasey | Hem5 user's guide[END_REF]. It is worthwhile to note that the HEM5 solver is constructed in a way that the constraint at the velocity level is satisfied at each stage: G(Q i )V i = 0 , ∀i = 1 . . . 8. In order to get the acceleration v and the Lagrange multiplier λ at the end of the time step, an additional linear system has to be solved

M -G T G 0 v λ = F r with r = -G qq (v, v). (15)
Concerning the HEM5 solver, V. Brasey and E. Hairer [START_REF] Brasey | Hem5 user's guide[END_REF] define an error based on the 7 th and 8 th estimations Q 7 and Q 8 such that

err 1 = q n+1 -Q 8 s = O(h 4 ) err 2 = q n+1 -q n -h( 5 2 Q 7 -3 2 Q 8 ) s = O(h 3 ). (16)
Finally, the practical error estimation is computed by

err = err 1 2 err 1 + 0.01 err 2 = O(h 5 ). (17)
The optimal step size is computed with the same rule ( 12) with p = 5.

An index-1 DAE solver: the Runge-Kutta-Fehlberg of order 4 (RKF45) method

The RKF45 is a 4 order method suited for the integration of ODEs. It belongs to the family of RK embedded methods and provides 2 approximations of the solution to evaluate the error. The scheme is based on 6 estimations of the derivatives defined as

Y i = y n + h f ( j<i a ij Ẏj ).
In our case, we define Y = col( q, q ) and

f = v -M -1 G (GM -1 G ) -1 ((Gv) q v + GM -1 F ) + M -1 F .
The 4 th order approximation of the solution at the end of the step is computed by y n+1 = y n + i≤6 b i Ẏj . A 5 th order estimation, used for the computation of the error, is given with ỹn+1 = y n + i≤6 bi Ẏj . The computational error is defined with err = ỹn+1 -y n+1 , and the computation of the optimal step size is given with [START_REF] Jean | The non smooth contact dynamics method[END_REF] with p = 4.

Numerical tests

In this section, numerical tests with different values of the tolerance on the truncation error will be performed on two mechanisms, using the time step control strategy presented in section 2. The aim of these tests is to compare the selected solvers in terms of computational effort and drift of the constraints. Table 1 presents the parameters of time-step control strategy. We set the tolerance on the drift of the constraints to a large value (10 -3 ) in order to prevent from a too severe requirement for some schemes when applying the time-step control strategy. The slider crank system is depicted in Figure 2a. Let m 1 , m 2 and l 1 , l 2 denote the masses and the lengths of the rods, and m 3 the mass of the slider. The system is described with the vector of generalized coordinates q = [α 1 , α 2 , x] as defined on Figure 2a. The revolute joint and the prismatic joint lead to two constraints which are written as:

g 1 (q) = l 1 sin(α 1 ) + l 2 sin(α 2 ) = 0 g 2 (q) = x -l 1 cos(α 1 ) -l 2 cos(α 2 ) = 0 (18)
The parameters of the simulation are: q 0 = [0, 0, 3] , v 0 = [8, -4, 0] , m 1 = 1 kg, m 2 = 1 kg, m 3 = 1 kg, l 1 = 1 m, l 1 = 2 m, gr = 9.81 m.s -2 (Earth gravity), t 0 = 0 s and T end = 2 s.

The rocking block

The aim is to study the motion of the center of mass of the block while maintaining an edge of the block in contact with the ground, as depicted in Figure 2b. The motion of the block can be described by the generalized coordinates q = [x, z, θ] . The characteristics of the simulation are: q 0 = [0.5L cos( π 3 ) -0.5l sin( π 3 ), 0.5L sin( π 3 ) + 0.5l cos( π 3 ), π 3 ] , v 0 = [0, 0, 0] , m 1 = 1 kg, l = 0.5, m, L = 1 m, gr = 9.81 m.s -2 (Earth gravity), t 0 = 0 s and T end = 0.5 s. During the simulation time, the point A of the block should stay in contact with the ground on a nontrivial time-interval. This constraint in A is formulated as g(q) = z -0.5 L sin(θ) -0.5 l cos(θ) = 0.

(19)

Results

Computational effort

In Figures 3a and3b, we present the precision-work diagrams for the two mechanisms. From the computational effort point of view, Figure 3 shows a big gap between the α-scheme on the one hand and the HEM5 solver and the RKF45 scheme on the other hand. The α-scheme solver is almost 9 times more expensive than HEM5 and 3 times more than the RKF45 scheme. We can say that overall the HEM5 solver is the most computationally efficient since the tolerances on the integration error and on the drift of the constraints are respected with smaller time step sizes.

Drift of the constraints

For each value of the tolerance in the range [10 -10 , 10 -2 ], we computed the maximum drift of the constraints at both position and velocity levels, for each mechanism and each scheme. Tables 3 and4 show the results for the slider-crank mechanism and for the rocking block.

Slider crank: With the HEM5 solver, the drift of the constraints at the velocity level vanishes (at the machine accuracy). The drift at the position level is acceptable even for large tolerances on the integration error. With the RKF45 method, the tolerance at the integration error is met with large time step sizes ( 10 -2 , 2 nd column of Table 3b) but the drift of the constraints is much higher than with the HEM5 solver. For the slider crank mechanism, the drift of g 1 and g 2 (columns 4 and 5 of Table 3c) using the α-scheme, vanishes (at the machine accuracy). The error at the velocity level is quite large (one order less than the time step size). If the magnitude of this drift exceeds the tolerance set for the index set I 2 , then the evaluation of this index set will be wrong or difficult to perform.

Rocking block: From Tables 4a, 4b and 4c, we can observe that for the large tolerances, the HEM5 solver and the RKF45 method are equivalent from the computational cost point of view, but this later is more costly for tight tolerances. But both solvers are doing better than the α-scheme. The drift of the constraints at the velocity level for the α-scheme is quite high for large tolerances (large time step sizes) by contrary to the two other solvers. This drift may lead to a wrong estimation of the closed contacts belonging to the set I 2 .

As said in the sections describing the HEM5 solver and the α-scheme, the first solver is constructed in a manner that enables to hold the constraints at the velocity level while the second one maintains the constraints at the position level, and both of them achieve very well these goals. The drift of the constraints using the RKF45 scheme is maintained at acceptable levels depending on the time step size being used. Overall, the HEM5 solver gives the best results with regards to the drift of the constraints at both position and velocity levels, even for large time step sizes. With regards to these results, it seems that discretizing the constraints at the velocity level is a good compromise to reduce their drift at both position and velocity levels.

Bodies with geometrical singularities

In this section, we will discuss the changes that occur when there are some edges or other discontinuities in a given geometry. These discontinuities make the integration of the equations of motion difficult, and could even lead to incoherent results if they are not correctly treated. In the case of the HEM5 solver and the RKF45 scheme, these singularities lead to a discontinuity in the term G qq of Equation ( 5), this leads to a jump in the acceleration and in the contact force. But these discontinuities also make the systems ( 14) and ( 10) ill-posed. Let us start with illustrating the problem with some examples. Then, we will propose some solutions to overcome this problem.

Illustration of the problem

Let us illustrate this with two problems presented in Fig. 4a and Fig. 4b. For the system of Fig. 4a, when the sphere reaches the connection, both constraints are active and with the HEM5 solver for example, the system ( 14) is not well-posed, here is a numerical example of such a system at the discontinuity:

    1 0 -1 -1 0 1 0 0 1 0 0 0 1 6.65031e -5 0 0     V1 Λ 1 =     -9.81 0 -3.62238e -2 4.41971e 1    
At x = 0, two constraints are active (g 1 : z = 0 and g 2 : z -x 2 = 0), but the Hessians of these constraints are not equal at the edge and this leads to a discontinuity in term G qq . In the next section, we propose two methods to solve this problem. The first solution consists in re-writing the linear systems ( 14) and (15) (in the case of HEM5, computation of the estimations of the accelerations and Lagrange multipliers) as a Linear Complementarity Problem (LCP)

     M (Q i ) Vi -G (Q i )Λ i = F (Q i , V i , t i ) G(Q i+1 ) Vi -r i ≥ 0 0 ≤ Vi ⊥ Λ i ≥ 0. (20)
This formulation transforms the unfeasible linear problem given by system [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF] with two equality constraints into a feasible LCP. This transformation enables to find a solution to the problem by relaxing one of the equality constraints at the edge. For our simulations, we used for the resolution of these LCPs Lemke's algorithm available in the open-source software SICONOS [START_REF] Acary | An introduction to Siconos[END_REF] 1 . When applied to the examples of Fig. 4a and 4b, Lemke's algorithm provides some estimations of the acceleration and the contact forces which correspond to the expected theoretical ones.

4.2.2 Solution 2: using the Moreau-Jean's time-stepping scheme

The 2 nd solution consists in using an event-capturing time-stepping scheme, for instance the Moreau-Jean's timestepping scheme [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF]. The event-capturing time-stepping scheme is used here to overcome the difficulty triggered by geometries with non-continuous curvatures. More precisely, it is used to compute the dynamics in the zones containing the geometrical discontinuities. The following pseudo-algorithm explains the way an event-capturing scheme method is used in the whole integration process:

• Detect the time t d when the discontinuity (edge, undetermined normal) occurs.

• Integrate the dynamics over the time-step [t d , t d + h] with the event-capturing scheme.

• Compute the set I d of active constraints at t d + h.

• Go back at time t d and solve the dynamics with a DAE solver with taking into account only the constraints belonging to I d .

We favor the Moreau-Jean method as a solution to the presented problems rather than the LCP (previous subsection) solution since it relies on a mechanical reasoning for finding the correct index set of active constraints on the right of a discontinuity. Thus, the results do depend only on the dynamics of the mechanical problem.

Conclusion

This article aimed at comparing index-3, index-2 and index-1 DAE solvers. We chose the generalized-α scheme for the first category, the HEM5 solver for the 2 nd one and the RKF45 method for the 3 rd one. We wanted to compare the solvers in terms of stabilization of the drift of the constraints, of computational effort and the problems they face in case of geometrical singularities. For this aim, we made simulations on two mechanisms, using fixed time steps and also with adapting the time step size to the user defined tolerance on the truncation error and the maximum drift. After analyzing the time step sizes used by each solver, the computational cost, the maximum drift of the constraints and the behavior in the zones containing singularities, we can draw the following conclusions:

• Computational effort. Even though the HEM5 solver contains more stages (8) than the RKF45 scheme (6), the numerical effort of the first solver is lower than that of the second one when a strategy of time step control is used. Indeed, the order (5) of the HEM5 solver in addition to its characteristics of reducing the violation of the constraints, enable to use larger step sizes than those used for the RKF45 scheme, and then to reduce the computational effort. Both HEM5 and RKF45 give better results than the generalized-α scheme that needs very small time step sizes mainly due to its lower order (2) and the drift at the velocity level. But we know that the generalized-α scheme is a low order scheme that is well suited for flexible multibody dynamics where high frequency non-linear dynamics can render the integration difficult. In this latter case, half-explicit methods can be in troubles.

• Drift of the constraints. The HEM5 solver enforces the constraints at the velocity level and it enables to reduce drastically the drift at the position level. This enables often to perform the integration without any procedure of projection of the constraints on the admissible manifold. This is not the case of the RKF45 scheme where a projection on the constraints is mandatory in some cases. The generalized-α scheme enforces the constraints at the position level but the drift at the velocity level may lead to numerically losing the contact depending on the tolerances that are chosen for the index set I 2 of closed contacts defined in Section 1. Discretizing the dynamics with a formulation of the constraint at the velocity level seems to be a good compromise to stabilize the drift at acceptable tolerances.

• Geometrical discontinuities. None of the solvers enable to avoid integration issues (systems ill-posed when constraints are not compatible,. . . ) that occur at some geometrical singularities. We presented two solutions to this problem: formulating the systems that become ill-posed at the singularities as LCPs and solve them with some dedicated LCP algorithm. To this aim, a range of LCP algorithms is proposed in the SICONOS software [START_REF] Acary | An introduction to Siconos[END_REF]. This solution gives correct results when the problem is feasible. The second solution is using en event-capturing timestepping scheme [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF] to compute the dynamics in the singularity zone. The last solution is promising since it is only driven by the dynamics and does not depend on any feasibility condition.

In the light of this preliminary study, some further work is planned. The PHEM56 solver of Murua (1995), which is a 6 stages partitioned half-explicit method of order 5 constructed in the same manner as the HEM5 solver, may be used instead of this latter to save the computational effort due to the 2 additional stages in HEM5. In this article, we worked mainly with semi-implicit schemes. Further works will deal with extrapolation solvers such as the MEXX solver by Lubich and al.(1992), and the use of hybrid mixed time stepping integrators in [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF] and projected integration schemes in [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF].

1

 1 The slider-crank system

Figure 2 :Figure 3 :

 23 Figure 2: Two academic mechanisms.

Figure 4 :

 4 Figure 4: Systems with geometrical discontinuities, illustrative examples

4. 2 Solutions to the problem 4 . 2 . 1 Solution 1 :

 24211 Formulating the system as a LCP

Table 1

 1

		: Parameters for time step control
	Absolute tolerance (tol) Minimum time step	Tolerance of	Maximum drift of	safety factor (safe)
			Newton's loop	g and	ġ
	[10 -10 , 10 -2 ](*)	10 -6 s	10 -10	10 -3		0.9

(*) We vary the value of tol to compare the computational effort and the drift of the constraints

Table 2

 2

summarizes the computational effort for a single time step for each solver.

Table 2 :

 2 Computational effort of the solvers during a single time step

	Solver	F	G	G q	K t	C t	# Linear equa-	Work per time
							tions / Matrices	step
							inversions
	HEM5	8	8	1	-	-	9	26
	RKF45	6	6	6	-	-	12	30

Table 3 :

 3 Slider-crank: average time step size, work and drift of the constraints for different tolerances

	(a) HEM5 solver

Table 4 :

 4 Rocking block: average time step size, work and maximum drift

			(a) HEM5 solver
	tolerance	average h	work	g max	ġmax
	10 -2	0.336	234	3, 47 10 -5	4, 81 10 -16
	10 -3	0.260	460	1, 13 10 -7	3, 34 10 -16
	10 -4	0.179	624	6, 34 10 -8	6, 72 10 -16
	10 -5	0.104	988	3, 16 10 -9	6, 97 10 -16
	10 -6	0.056	1924	4, 76 10 -10 4, 58 10 -16
	10 -7	0.054	2524	3, 78 10 -10 6, 45 10 -16
	10 -8	0.046	3120	6.00 10 -11	6, 13 10 -16
	10 -9	0.026	5824	8, 01 10 -11 6, 84 10 -16
	10 -10	0.018	8580	1, 29 10 -11 7, 14 10 -16
			(b) RKF45 method
	tolerance	average h	work	g max	ġmax
	10 -2	0, 065	390	1, 70 10 -4	9, 50 10 -5
	10 -3	0, 063	450	2, 21 10 -5	8, 07 10 -5
	10 -4	0, 051	630	1, 21 10 -5	8, 07 10 -5
	10 -5	0, 039	840	2, 02 10 -6	2, 38 10 -6
	10 -6	0, 0009	3036	5, 62 10 -7	9, 98 10 -7
	10 -7	0, 017	2700	4, 93 10 -8	9, 75 10 -8
	10 -8	0, 0029	13560	7, 13 10 -9	1.00 10 -8
	10 -9	0, 0008	5910	4, 17 10 -10 9, 96 10 -10
	10 -10	7, 6 10 -4	48510	5, 30 10 -11 1.00 10 -10
			(c) α-scheme
	tolerance	average h	work	g max	ġmax
	10 -2	7.11 10 -3	2180	8.33 10 -17	9.33 10 -5
	10 -3	7.11 10 -3	2180	8.33 10 -17	9.33 10 -5
	10 -4	1.26 10 -3	19815	8.33 10 -17	1.43 10 -3
	10 -5	6.52 10 -4	30870	8.33 10 -17	2.28 10 -4
	10 -6	3.12 10 -4	55725	8.33 10 -17	2.79 10 -5
	10 -7	1.28 10 -4	129555	8.33 10 -17	5.77 10 -6
	10 -8	4.86 10 -5	303435	8.33 10 -17	4.12 10 -7
	10 -9	1.97 10 -5	667055	8.33 10 -17	5.30 10 -8
	10 -10	6.48 10 -6	1241790	1.11 10 -16	7.50 10 -9

http://siconos.gforge.inria.fr


Mounia Haddouni
email: mounia.haddouni@inria.fr

Vincent Acary
email: vincent.acary@inria.fr

Jean-Daniel Beley
email: jean-daniel.beley@ansys.com

Stéphane Garreau
email: stephane.garreau@ansys.com

Comparison of index-2 and index-1 DAE solvers for nonsmooth multibody systems with unilateral or bilateral constraints

Keywords: Nonsmooth mechanics, event-driven schemes, constraints, DAEs, ODEs, numerical solvers, LCP, NSCD

come

Introduction

Using the Lagrangian formalism, the dynamics of a rigid multibody system with m frictionless contact points can be described by

     q(t) = v(t) (1a)
M(q(t)) v(t) = F(t, q(t), v(t)) + G T (t, q(t))λ(t) (1b) 0 ≤ λ(t) ⊥ g(q(t)) ≥ 0,

where q(t) ∈ R n is the generalized coordinates vector, v(t) ∈ R n is the generalized velocities vector, M(q(t)) ∈ R n×n is the symmetric definite positive matrix of inertia, F(t, q(t), v(t)) comprises the external applied loads, the non linear interactions between bodies and the non linear inertial terms, g(q(t)) ∈ R m is the vector of constraints, G = ∂ ∂q g(q(t)) is the Jacobian of the constraints and λ(t) ∈ R m is the Lagrange multiplier which is related to the force at the contact point. The complementarity condition (1c) models the unilateral contact, or the Signorini condition at the position level. We suppose that impacts occur in an infinitely short period so that the displacements of the bodies during the contact period can be negligible and we use the Newton impact law with a coefficient of restitution e. If g(q(t)) = 0, we compute the impulse p i and the post-impact velocity v + (t i ) by solving the impact equations

       M(q(t i ))(v + (t i ) -v -(t i )) = G T (t i , q(t i ))p i U + N (t i ) = G(t i , q(t i ))v + (t i ) U - N (t i ) = G(t i , q(t i ))v -(t i ) 0 ≤ U + N (t i ) + eU - N (t i ) ⊥ p i ≥ 0, (2)
where U N (t) = dg dt (q(t)) = G(q(t)) q(t) is the relative normal velocity at contact. The complementarity condition in (2) describes the Signorini condition written at the velocity level, augmented of the impact law.

When an event-driven strategy is chosen for the time integration strategy of nonsmooth multibody systems with unilateral or bilateral constraints, any solver for differential algebraic equations (DAE) can be used between two events for the time integration of the smooth systems with equality constraints g(q) = 0. Three index sets are generally introduced in order to characterize the state of the contacts:

• The index set I 0 of all possible constraints in the system: I 0 = {1, ..., m}

• The index set I 1 of all closed contacts:

I 1 = {i ∈ I 0 |g i (q(t)) = 0}
• The index set I 2 of all persistent contacts:

I 2 = {i ∈ I 1 | ġi (q(t)) = 0}
From the numerical point of view, the smooth dynamics is integrated with some DAE solver, during periods of time over which the index set I 2 is constant. An event detection is performed to find the time where a new contact occurs or a contact needs to be released. When an event is detected, the impact equations ( 2) are solved and the state of the dynamical system is updated. A complete algorithm for the event-driven schemes is detailed in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]. As we said above, between two events, we have to solve the index-3 DAE (3) over the closed contact set

I 2      q(t) = v(t) (3a) M(q(t)) v(t) = F(t, q(t), v(t)) + G T (t, q(t))λ (3b) g α (q(t)) = 0, λ α ≥ 0, , α ∈ I 2 . ( 3c
)
It is well known that index-3 DAEs are difficult to numerically handle. Therefore, the dynamics is usually integrated as an ODE by reducing the original index 3 of the system to 1. It amounts to solving the problem at the acceleration level by differentiating twice the constraints. Index reduction consists in differentiating the constraint as many times as necessary to get a set of equations that may be solved using methods for lower index problems. Hence, if the constraint g is differentiated once with respect to time, one obtains an index-2 DAE, we solve in this case (3a) and (3b) together with

G(t, q(t))v(t) = 0. ( 4
)
If g is twice differentiated, one gets an index-1 DAE composed of (3a) and (3b) with

G(t, q(t)) v(t) + G qq (v, v) = 0. (5)
Numerically, the index reduction leads to the phenomenon of drift of the constraints, at different levels depending on the index. Indeed, as opposed to the continuous time case, enforcing (5) in discrete time does not imply that (3c) and (4) are satisfied.

2 Quick recall on the numerical time-integration schemes 2.1 An index-2 DAE solver. The HEM5 solver

The HEM5 solver is based on a half-explicit 5 order Runge-Kutta method with 8 stages, described in [START_REF] Brasey | Hem5 user's guide[END_REF]. Half-Explicit Runge-Kutta methods have been introduced by Ernst Hairer in [START_REF] Hairer | Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems[END_REF]. The HEM5 solver is suited to the integration of DAEs, reduced to the index 2 formulation. Let us denote q n , v n and t n the values of position, velocity and time at the beginning of a time step, and let h be the length of the time step. The numerical approximation is given by

   M(Q i ,t n + c i h) Vi = F(t n + c i h, Q i ,V i ) + G T (t n + c i h, Q i )Λ i Qi = V i G(t n + c i h, Q i )V i = 0, (6)
where the stages are evaluated as follows :

Q i = q n + h ∑ j<i a i j Q j and V i = v n + h ∑ j<i a i j Vj .
At each stage, the estimations of position and velocity Q i and V i are explicitly computed, while Λ i and Vi are implicitly obtained by solving the following linear system

M i -G T i G i+1 0 Vi Λ i = F i r i , (7)
where

r i = - G(Q i+1 ) ha i+1,i (v n + h i-1 ∑ j=1 a i+1, j Vj ).
At the end of the time step, the numerical solution is given by

q n+1 = Q 9 = q n + h ∑ 8 i=1 b i Qi and v n+1 = V 9 = v n + h ∑ 8 i=1 b i Vi , with b i = a 9i
. The computation of coefficients c i and a i j is detailed in [START_REF] Brasey | Hem5 user's guide[END_REF]. It is worthwhile to note that the HEM5 solver is constructed in a way that the constraint at the velocity level is satisfied at each stage: G(Q i )V i = 0 , ∀i = 1..8. In order to get the acceleration and the Lagrange multiplier at the end of the time step, an additional linear system has to be solved:

M -G T G 0 v λ = F r , (8)
where r = -G qq (v, v).

The embedded schemes technique is the most suited way to compute the optimal time step when we deal with Runge-Kutta schemes. Concerning the HEM5 solver, V.Brasey and E.Hairer [START_REF] Brasey | Hem5 user's guide[END_REF] define an error based on the 7 th and 8 th estimations such that

err = err 1 2 err 1 + 0.01err 2 = O(h 5 ), (9)
with

err 1 = q n+1 -Q 8 s = O(h 4 ) and err 2 = q n+1 -q n -h( 5 2 Q 7 - 3 2 Q 8 ) s = O(h 3 ).
Finally, The optimal step size is computed with (10):

h opt = s h ( tol h err ) 1 p , (10)
with p = 5 is the local order of the method, s a safety factor, and tol the user defined tolerance.

A standard ODE solver. The Runge-Kutta-Fehlberg scheme

This method is suited for the integration of ODEs. It belongs to the family of Runge-Kutta-Embedded methods and provides 2 approximations of the solution to evaluate the error. The scheme is based on 6 estimations of the derivatives defined as:

Y i = y n + h F( ∑ j<i a i j Ẏj ).
In our case, we have:

Y = q q and F = v -M -1 G T (GM -1 G T ) -1 ((Gv) q v + GM -1 f ) + M -1 f .
A 4 th and 5 th order approximations of the solution at the end of the step are computed by y n+1 = y n + ∑ i≤6 b i Ẏj and ỹn+1 = y n + ∑ i≤6 bi Ẏj . The computational error is defined with:

err = ỹn+1 -y n+1 . (11)
Finally, the computation of the optimal step size is given by (10) with p = 4.

3 Numerical study of the control of the drift of the constraints

The four-bar linkage

In this section we consider the four-bar linkage system, with 4 bars linked with revolute joints, and one of them being grounded, as described in Fig. 1: Fig. 1: Four-bar linkage

This system is driven by a torque of value τ = 6 N.m and is described by the vector of generalized coordinates q = [ϕ 1 , ϕ 2 , ϕ 3 ] as depicted on the figure. During the simulation, the lengths of the bars must be kept constant. To this aim, two constraints are defined:

g 1 (q) := l 1 cos(ϕ 1 ) + l 2 cos(ϕ 2 ) -l 3 cos(ϕ 3 ) -d 1 = 0 g 2 (q) := l 1 sin(ϕ 1 ) + l 2 sin(ϕ 2 ) -l 3 sin(ϕ 3 ) = 0. ( 12
)
The initial conditions of position and velocity are given by q 0 = [ π 2 , 0, π 6 ] T and v 0 = [0, 0, 0] T . They are compatible with the constraints at both position and velocity levels. The time step size is controlled in order to meet a tolerance of 10 -10 on the integration error but also to maintain the drift of the constraints (at position and velocity levels) at an acceptable level (10 -6 ). We set the minimum step size at 10 -6 s and the maximum step size at 10 -2 s. The drift of the constraints at both position and velocity levels, for both solvers, are plotted in Fig. 2 and Fig. 3. For the levels of accuracy that we set, the Runge-Kutta-Fehlberg scheme uses the minimum step size (10 -6 s) at every time step. This means that the prescribed tolerance is hardly achieved by the scheme in this range of time-steps. The time step sizes needed for the simulation with the HEM5 and the RKFehlberg schemes are presented in Fig. 4. Thanks to its intrinsic characteristic of maintaining the constraint at the velocity level, the HEM5 solver reduces the drift of the constraint drastically at both position and velocity levels, and for a much lower computational effort compared to the Runge-Kutta-Fehlberg scheme.

The rocking block

The aim is to study the motion of the center of mass of the block depicted in Fig. 5.

Fig. 5: Rocking block

The motion of the block can be described by the generalized coordinates q = [x, z, θ]. During the simulation time (0.5s), the point A of the block should be in contact with the ground. We varied the tolerances on the truncation error and on the drift of the constraints, and then we varied the time step size. The Work/Precision and the error/h diagrams are given in Fig. 6a and Fig. 6b. We can see that holding the drift of constraints under very low thresholds enables the HEM5 solver to reduce the computational effort for high user-defined tolerances. The Runge-Kutta-Fehlberg needs lower step sizes to meet the tolerance on the drift of the constraints. We varied the tolerance for both HEM5 and RKfehlberg schemes, and we computed the average time step size, the computational effort1 and the maximum drift of g and ġ. Results are detailed in

Problems with geometrical discontinuities

In this section, we discuss the problems that occur when in a given geometry, there are some edges or other discontinuities in the contact interface that yield nonsmoothness in the constraint. These discontinuities make the integration of the equations of motion more difficult, and could even lead to incoherent results if they are not correctly treated. In the case of problems with geometrical discontinuities (edges for example), the term G qq of Equation ( 5) is discontinuous, this leads to a jump in the acceleration and in the contact force. Numerically, these discontinuities result in a difficulty when computing the estimations of acceleration or Lagrange multipliers because the system ( 7) is not well-posed. Let us illustrate this with two problems presented in Fig. 7a and Fig. For the system of Fig. 7a, when the sphere arrives at the edge, both constraints are active and with the HEM5 solver for example, the system ( 7) is not well-posed, here is an example of such a system:

    1 0 -1 -1 0 1 0 0 1 0 0 0 1 6.65031e -5 0 0     V1 Λ 1 =     -9.81 0 -3.62238e -2 4.41971e 1    
At x = 0, two constraints are active (g1 : z = 0 and g2 : zx2 = 0), but the Hessian of these constraints are not equal at the edge and leads to a discontinuity in term G qq . We propose two methods to solve this problem.

1 st method: Formulating the system as an LCP The first solution consists in re-writing the systems for the computation for the estimations of accelerations and Lagrange multipliers (for example system [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF] for the HEM5 solver) as a Linear Complementarity Problem:

     M(Q i ) Vi -G T (Q i )Λ i = F i G(Q i+1 ) Vi -r i ≥ 0 0 ≤ Vi ⊥ Λ i ≥ 0. ( 13
)
This formulation transforms the unfeasible linear problem given by system [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF] with two equality constraints into a feasible complementarity problem, which enables to find a solution to the problem. For our simulations, we used for the resolution of these LCPs Lemke's algorithm available in the open-source software SICONOS 2 . When applied to the examples of Fig. 7a and 7b, Lemke's algorithm provides some estimations of the acceleration and the contact forces which correspond to the expected theoretical ones.

2 nd method: Using the NSCD method The other solution that we propose is to compute the dynamics in the zone that presents discontinuities with a time-stepping scheme, for instance the Non Smooth Contact Dynamics (NSCD) method (Moreau-Jean's time-stepping scheme [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF]). The NSCD method is used here to overcome the difficulty triggered by geometries with non-continuous curvatures. More precisely, it is used as a way to sort the constraints to be activated, that is decide the 'direction' to go through when crossing an edge for example. The following pseudo-algorithm explains the way the NSCD method is used in the whole integration process: We favor the NSCD method as a solution to the presented problems rather than the LCP solution (previous subsection) since it is a physical solution, so the results are only depending on the dynamics of the given problem.

Conclusion

Even though the HEM5 solver contains more stages (8) than the Rune-Kutta-Fehlberg scheme (6), the numerical effort of the first solver is lower than that of the second one when they are used with an adaptive time step. Indeed, the order (5) of the HEM5 solver in addition to its characteristics of reducing the violation of the constraints, enables using smaller time step sizes than those used for the RKFehlberg scheme, and thus reduces the computational effort. The HEM5 solver enforces the constraints at the velocity level and it enables to reduce drastically the drift at the position level. This enables often to perform the integration without any procedure of projection of the constraints on the admissible manifold. This is not the case of the RKFehlberg scheme where a projection on the constraints is mandatory. Discretizing the dynamics with a formulation of the constraint at the velocity level seems to be a good compromise.

None of the solvers enable to avoid integration issues (systems not well-posed when constraints are not compatible) that occur at some geometrical discontinuities. We presented two solutions to this problem: using Lemke's algorithm which gives correct results when the problem is feasible, or using the NSCD method. The last solution is promising since it is only driven by the dynamics and does not depend on any feasibility condition.

Further works will deal with index-3 DAE resolution as in the work of O.Bruls and M.Arnold in [START_REF] Brüls | Convergence of the generalized-α scheme for constrained mechanical systems[END_REF], and the use of hybrid mixed time stepping integrators in [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF] and projected integration schemes in [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF].

  (a) g 1 versus time. (b) g 2 versus time.

Fig. 2 :

 2 Fig. 2: Four bar linkage example. Drift of the constraints at the position level

Fig. 3 :

 3 Fig. 3: Four bar linkage example. Drift of the constraints at the velocity level

Fig. 4 :

 4 Fig. 4: Four-bar linkage. Evolution of the time step with the simulation time.

  Fig. 6: error and precision diagrams for the rocking block

  7b

  .

Fig. 7 :

 7 Fig. 7: Systems with geometrical discontinuities, illustrative examples

Table 1

 1

	and 2.
	tolerance average h work	g max	ġmax
	10 -2	0.336	234	3, 47.10 -5 4, 81.10 -16
	10 -3	0.260	460	1, 1310 -7	3, 3410 -16
	10 -4	0.179	624	6, 3410 -8	6, 7210 -16
	10 -5	0.104	988	3, 1610 -9	6, 9710 -16
	10 -6	0.056	1924 4, 7610 -10 4, 5810 -16
	10 -7	0.054	2524 3, 7810 -10 6, 4510 -16
	10 -8	0.046	3120	610 -11	6, 1310 -16
	10 -9	0.026	5824 8, 0110 -11 6, 8410 -16
	10 -10	0.018	8580 1, 2910 -11 7, 1410 -16
	10 -11	0.010	15548 4, 1710 -11 9, 1710 -16

Table 1 :

 1 Rocking block: average h, work and maximum drift with HEM5

	tolerance average h	work	g max	ġmax
	10 -2	0, 065	390	1, 710 -4	9, 510 -5
	10 -3	0, 063	450	2, 2110 -5	8, 0710 -5
	10 -4	0, 051	630	1, 2110 -5	8, 0710 -5
	10 -5	0, 039	840	2, 0210 -6	2, 3810 -6
	10 -6	0, 0009	3036	5, 6210 -7	9, 9810 -7
	10 -7	0, 017	2700	4, 9310 -8	9, 7510 -8
	10 -8	0, 0029	13560	7, 1310 -9	10 -8
	10 -9	0, 0008	5910	4, 1710 -10 9, 9610 -10
	10 -10	7, 6 -4	48510	5, 310 -11	10 -10
	10 -11	1, 9710 -6 9585990 1.2310 -11 9, 7910 -12

Table 2 :

 2 Rocking

block: average h, work and maximum drift with Runge-Kutta-Fehlberg

  1. Detect time t d when the discontinuity (edge, undetermined normal) occurs 2. Compute the dynamics with the NSCD method 3. Compute the set I d of active constraints at the end of the step 4. Go back at time t d and solve the dynamics with a DAE solver with taking into account only the constraints belonging to I d

The computational effort is computed by counting the number of function evaluations and the number of linear equations that are solved.

http://siconos.gforge.inria.fr