Van-Vuong Lai
email: van-vuong.lai@polytech-lille.fr

Olivier Chiello
email: olivier.chiello@ifsttar.fr

Jean-François Brunel

Philippe Dufrenoy

V-V Lai

J.-F Brunel

P Dufrénoy

Dynamic FEM simulation of wheel-rail rolling contact with friction in an Eulerian frame -Application to curve squeal

Introduction

In the case of a tight radius curve, squeal can be generated due to the unstable slip of the wheel on the head of the rail. The mechanism causing the instability is however still controversial. A negative slope of the friction coefficient depending on creepage velocity is introduced in some models as a source of the instability whereas squeal can also occur in the case of constant friction due to mode coupling instability in other models.

There are actually wheel/rail rolling contact models with some simplifications. Equivalent point-contact models are based on analytic formulas or heuristic laws of frictional rolling contact [START_REF] Ayasse | 4 Wheel-Rail Contact, Handbook of railway vehicle dynamics[END_REF][START_REF] Jj | Wheel-rail rolling contact theory[END_REF]. They give good results in slightly non linear dynamic cases close to pure rolling [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]. However, in the cases of large slips and high-frequency dynamics, surface contact models are more adapted. They need a discretization of the contact zone. The variational theory of Kalker [START_REF] Jj | Wheel-rail rolling contact theory[END_REF] which is used in Pieringer's model [START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF] performs some simplifications as the semi-analytical computation of the local contact flexibilities or influence functions (Boussinesq and Cerruti elastic half-space assumptions, contact/friction decoupling). In the case of frictioninduced vibrations in curves characterized by large slips, high frequency dynamics and very fast evolutions in the contact zone, the impact of these simplifications is unknown.

Finite element wheel/rail frictional rolling contact models use commercial softwares [START_REF] Dahlberg | Transient rolling of cylindrical contacts with constant and linearly increasing applied slip[END_REF][START_REF] Toumi | Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods[END_REF] but the calculation is generally performed in a Lagrangian reference frame, which greatly limits its performance and is not adapted to high frequency dynamics.

The aim of the paper is to propose a full finite Element method for modeling friction induced vibration of the wheel-rail contact in curves. The wheel-rail contact is simulated by surface contact models with a discretization of the contact zone and non-smooth frictional contact laws. To reduce computational costs, simulations are performed in a Eulerian reference frame and modal reduction strategies are used. Instabilities are analyzed through 2 types of analysis: stability analysis to predict unstable modes assuming a full steady steading in the contact zone, and a transient analysis which allows to introduce non linearities ans to determine amplitude of vibrations. By using finite element method, discretization of the contact surface, non-smooth frictional contact laws and model reduction techniques in an Eulerian frame. Reduction base including free-interface modes and static attachment modes is used to reduce computing duration. The proposed approach is applied for contact between the rail and the wheel during curving.

2 Finite Element (FE) formulation in the time-domain

Contact problem in the Eulerian frame

A wheel/rail rolling contact is considered, as shown in Fig. 1. The z-axis is chosen to coincide with the common normal to the two surfaces in contact, the longitudinal x-axis corresponds to the rolling direction and the y-axis refers to the lateral direction. The wheel rolling speed is V in a direction parallel to x-axis . In case of curving, wheel relative lateral velocities ∆V y and longitudinal velocities ∆V x have to be considered. The wheel may also have a relative angular velocity ∆ω z around their common normal (or spin). According to [START_REF] Johnson | Contact mechanics[END_REF], in the Eulerian frame which moves with the point of contact, the relative sliding instantaneous velocities (or creep velocities) between the wheel and the rail at a fixed point of the potential contact interface are given by:

ṡx = v W x -v R x = ∆V x -∆ω z y + V ( ∂u W x ∂x - ∂u R x ∂x ) + ( uW x -uR x ) ṡy = v W y -v R y = ∆V y + ∆ω z x + V ( ∂u W y ∂x - ∂u R y ∂x ) + ( uW y -uR y ) (1)
where R, W denote respectively the rail and the wheel. u(x, y, z, t) and v(x, y, z, t) denote respectively the displacement and velocity fields of the structure in the Eulerian frame and the notation u = ∂u ∂t refers to the time partial derivative. The terms involving V ∂u ∂x represent the deformation contributions due to rolling in the Eulerian frame whereas the terms involving u simply represent the dynamic contributions.

Contact laws

To deal with unilateral contact on the interface, a non-regularized Signorini law is chosen:

∆u n -g ≤ 0 r n ≤ 0 (∆u n -g)r n = 0 (2)
where ∆u n = u W n -u R n is the normal relative displacement on the potential contact interface, r n is the normal contact stress and g is the initial gap. This law simply conveys that (i) there is no interpenetration between the two bodies, (ii) the interface only undergoes compression and (iii) it respects the condition of complementarity. A equivalent semi-regularized form of the law uses the projection on the negative real set (cf. for instance [START_REF] Jean | The non-smooth contact dynamics method[END_REF]):

r n = Proj R -(r n -ρ n ∆u n ) ∀ρ n > 0 (3)
where Proj R -(x) = min(x, 0) and ρ n is a positive real number called normal augmentation parameter.

To deal with frictional contact, a non-regularized Coulomb law with a constant friction coefficient µ is used:

r t ≤ -µr n r t = -µr n ⇒ ∃λ > 0, ṡt = -λr t r t < -µr n ⇒ ṡt = 0 ( 4
)
where ṡt is the vectorial creep velocity with components ṡx and ṡy defined in equation ( 1) and r t is the tangential stress due to friction. This law states that the frictional reaction cannot be greater than a limit -µr n . If the limit is reached, the material particle slides and has a direction opposed to the direction of the reaction (slip state). In other cases, the relative velocity is null (stick state).

A equivalent semi-regularized form of the law uses the projection on the Coulomb cone (cf. for instance [START_REF] Jean | The non-smooth contact dynamics method[END_REF]):

r t = P roj C (r t -ρ t ṡt ) ∀ρ t > 0 (5)
where P roj C (x) = min( µ|rn| x , 1)x and ρ t is a positive real number called tangential augmentation parameter.

Dynamics equations and FE discretization

Convective terms in the equations of motion induced by the Eulerian frame are only significant in the case where the rotational speed is in the order of magnitude of the eigenfrequencies ( [START_REF] Chambrette | Stability of a beam rubbed against a rotating disc[END_REF]). In the following application, the effects of the convective terms are neglected because the wheel rotational speed is assumed low in comparison to the wheel/rail eigenfrequencies.

The displacement field u must verify the equations of continuum mechanics [START_REF] Raous | Numerical characterization and computation of dynamic instabilities for frictional contact problems, Friction and Instabilities[END_REF] together with the local formulations of contact equations Eqs. (3) and [START_REF] Dahlberg | Transient rolling of cylindrical contacts with constant and linearly increasing applied slip[END_REF]. In order to introduce FE approximations, the principle of virtual power is used for the system dynamics and the contact laws are written in weak forms [START_REF] Kudawoo | Problèmes industriels de grande dimension en mécanique numérique du contact: performance, fiabilité et robustesse[END_REF]:

Find u ∈ U and r such as ∀u * ∈ U o and ∀r *

Ω ρu * .üdΩ + Ω (u * ) : σ(u)dΩ = Ω u * .f s dΩ + ∂Ω F u * .f d ds + Sc u * .rds Sc r * n r n ds = Sc r * n Proj R -(r n -ρ n ∆u n ) ds Sc r * t r t ds = Sc r * t Proj C (r t -ρ t ṡt ) ds ( 6
)
where Ω is the whole structure. f s and f d are respectively a volume load in Ω and a surface load on a surface δΩ F . (u) and σ(u) stand respectively for the symmetric gradient and the stress tensor corresponding to any displacement field u, the notation ü = ∂ 2 u ∂t 2 refers to the double time partial derivative, U = {u|u = u d on ∂Ω u } and U 0 = {u|u = 0 on ∂Ω u }. δΩ u is the part of the boundary of the domain Ω where a displacement u d is prescribed on. S c is the potential contact zone. By using appropriate shape interpolation functions for unknown and test displacement and reaction fields (u, r, u * , r * ), finite element discretization of Eq. ( 6) directly gives:

M Ü + C U + KU = F + P n T R n + P t T R t R n = Sc N n T Proj R -(N n (H -1 n R n -ρ n (P n U -G))ds R t = Sc N t T Proj C (N t (H -1 t R t -ρ t Ṡt (P t U, P t U)))ds (7)
where U, R n and R t denote respectively the vectors of nodal displacements, equivalent normal reactions and equivalent friction forces, M, C, K are respectively the mass, damping and stiffness matrices of the structure without contact, G is the vector of nodal initial gaps and N n , N t are the shape function vectors on the contact interface. In addition, P n , P n are matrices allowing to pass the contact reactions from the local relative frame to the global frame whereas H n = Sc N n T N n ds and H t = Sc N t T N t ds are transformation matrices from nodal to equivalent forces. Finally, Ṡt (P t U, P t U) denotes the vector of nodal creep velocities which can be determined linearly from local displacement and velocity vectors, taking into account quasistatic creep velocities (cf. Eq. ( 1)).

Reduction strategies

The principle is to search an approximated solution U = Bq r of the problem spanned by a reduced basis B which leads to a reduced dynamics equation:

M r qr + C r qr + K r q r = B T (F + P n T R n + P t T R t ) (8)
where M r , C r , K r = B T (M r , C r , K r )B and the size of the system is reduced to the number of modes in basis B.

In this paper, the proposed basis B includes free-interface normal modes of the structure and static attachment modes B = [Φ Φ s ]. Matrix Φ contains the real solutions of the free and undamped system:

(K -ω 2 M)U = 0 (9)
whereas Φ s make the base statically complete by adding static solutions to unitary forces (normal and tangential) on the contact interface:

KΦ s = P c T (10)
An approximation consists in neglecting the dynamic terms in the reduced equations corresponding to attachment modes. However, bases Φ and Φ s are not orthogonal which is a necessary condition in order to separate the contributions of normal modes and attachments modes. This can be fixed by using residual attachment modes defined by:

Φs = Φ s -Φ(Φ T KΦ) -1 Φ T P c T ( 11
)
These modes are the static displacement responses to unit contact reactions after the elimination of the contribution of normal modes. It can be easily verified that ΦT s KΦ = 0. The use of basis B = [Φ Φs ] in Eq. ( 8) together with the elimination of the dynamic terms relating to attachments modes gives:

Φ T MΦq + Φ T CΦ q + Φ T KΦq = Φ T (F + P n T R n + P t T R t ) ΦT s K Φs q s = ΦT s (F + P n T R n + P t T R t ) (12)
where q and q s are the generalized coordinate vectors corresponding respectively to normal modes and residual attachment modes such that U = Φq + Φs q s . A total decoupling between the generalized equations corresponding to free-interface normal modes and residual static attachment modes is obtained.

In addition, noticing that:

ΦT s K Φs q s = (Φ T s -P c Φ(Φ T KΦ) -1 Φ T )K(Φ s -Φ(Φ T KΦ) -1 Φ T P c T )q s = (P c Φ s -P c Φ(Φ T KΦ) -1 Φ T P c T )q s = P c Φs q s = P c (U -Φq) (13)
and including the above expression in Eq. ( 8) gives:

Φ T MΦq + Φ T CΦ q + Φ T KΦq = Φ T (F + P n T R n + P t T R t ) P c U = ΦT s (F + P n T R n + P t T R t ) + P c Φq (14)
It must be emphasized that, with this formulation, contact displacements P c U can be directly calculated from expanded forces and normal modes coordinates.

This reduction strategy consists in solving the global dynamics using free-interface normal modes and adding a local static residual flexibility controlled by matrix ΦT s in the expression of the contact displacements. This method gets very close to models such Pieringer's one [START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF] combining modal dynamics and static local resolutions of the kind of CONTACT software. However, in the proposed method, the local flexibility is calculated numerically by finite elements instead of using semi-analytical influence functions derived from Boussinesq's approximations. In addition, it takes into account the static normal/tangential coupling in the contact zone.

Stability analysis

The aim of the stability analysis is to address the mechanism of curve squeal due to frictional contact through the determination of the evolution of small perturbations around the steady sliding equilibrium.

Quasi-static equilibrium

The quasi-static or steady sliding equilibrium is first obtained by neglecting the dynamic terms in Eqs. ( 1) and [START_REF] Johnson | Contact mechanics[END_REF]:

KU e = F + P T n R e n + P T t R e t R e n = Sc N T n Proj R -(N n (H -1 n R e n -ρ n (P n U e -G)))ds R e t = Sc N T t Proj C (N t (H -1 t R t -ρ t Ṡe t )ds (15)
where Ṡe t = Ṡt (P t U e , 0) is the vector of nodal quasi-static creepage velocities. Assuming that some solutions of Eq. ( 15) exist and can be calculated, it notably provides the status of the nodes on the contact interface as a function of the equivalent normal reactions r e n and friction forces r e t . The facing nodes are in contact if r e n < 0. These contact nodes are then sliding if r e t = -µr e n and sticking if r e t < -µr e n .

Complex Eigenvalue Analysis (CEA) in case of full steady sliding

Stability analysis is only carried out in the case of full steady sliding (no sticking region) and maintained contact configuration: it is thus assumed that for each node in contact at equilibrium, bilateral contact and sliding Coulomb friction laws apply. In order to perform the stability analysis, these laws have to be linearized. On effective contact region, the linearized forms of Eqs. ( 3) and ( 5) with the above assumptions can be written (cf. for instance [START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF]):

∆u n = 0 r t = -µr n t -c b ub b (16)
where c b = -µr e n / ṡe t is a damping term due to the linearisation of the sliding direction of the friction force and b is the tangential direction orthogonal to t.

Searching a discrete solution of the form U e + Ũ exp(λt) where Ũ stands for the complex displacement vector corresponding to small harmonic perturbations around the equilibrium, the linearized form of equation ( 7) leads to a constrained non symmetric eigenvalues problem:

(λ 2 M + λ(C + C b ) + K) Ũ = ( PT n + µ PT t ) Rn Pn Ũ = 0 ( 17
)
where C b is the damping matrix provided by the linearisation of the sliding direction of friction force [START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF][START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. Application to an elastic layer[END_REF][START_REF] Lorang | TGV disc brake squeal[END_REF] and Pn , Pt are new projections matrices such that Pn is the restriction of P n on nodes in the effective contact region at equilibrium whereas Pt is the restriction of P t on components in direction t on nodes in the effective contact region. In stable cases, the perturbations vanish and no vibration occur. In unstable cases, Complex modes and eigenvalues of the problem are then calculated. Modes corresponding to eigenvalues with negative real part are stable. In this case, the perturbations corresponding to these modes vanish and no vibration occur. On the other hand, modes corresponding to eigenvalues with positive real part are unstable. Some perturbations corresponding to these modes tend to diverge which can lead to self-sustained vibrations.

The divergence rate of a mode is notably defined as Re(λ)/Im(λ) where (Re(λ), Im(λ)) are respectively the real and imaginary parts (pulsation) of the mode. This rate corresponds to a negative damping rate.

Reduced CEA formulations

Solving such a large non-symmetric eigenvalue problem needs model reduction. The approximation consisting in neglecting the dynamic part of the attachments modes is also tested for the stability analysis. The same technique used in section 3 is applied. The reduction is first performed with basis B = [Φ Φsn ] composed of free-interface component modes Φ and static attachment residual modes Φsn defined on effective contact region by:

Φsn = Φ sn -Φ(Φ T KΦ) -1 Φ T PT n (18)
The dynamic effects corresponding the residual attachment modes are then neglected leading to the following constrained eigenvalue problem:

(λ 2 Φ T MΦ + λΦ T (C + C b )Φ + Φ T KΦ)q = Φ T ( PT n + µ PT t ) Rn Pn Ũ = ΦT sn ( PT n + µ PT t ) Rn + Pn Φq Pn Ũ = 0 (19)
where q is the generalized vector corresponding to free-interface normal modes.

From the two last lines of Eq. ( 19), normal reactions Rn can be expressed as a function of the generalized vector q:

Rn = -ΦT sn ( PT n + µ PT t )) -1
Pn Φq (20)

Finally, using the above expression in the reduced eigenvalue problem gives:

(λ 2 Φ T MΦ + λΦ T (C + C b )Φ + Φ T (K + K c )Φ)q = 0 (21)
where K c is a non-symmetric stiffness matrix taking account the effect of the local residual flexibility of the structure due to normal and friction forces:

K c = ( PT n + µ PT t )(Φ T sn ( PT n + µ PT t )) -1 Pn (22)
This reduction strategy consists in solving the global dynamics using free-interface normal modes and adding a local static residual stiffness controlled by matrix K c . This matrix plays the same role than 2-DOF Hertzian stiffness/Coulomb friction relations in the case on equivalent point contact models but is calculated by finite elements.

Numerical methods for non linear simulations

For the computation of the quasi-static solution (Eq. ( 15)) and dynamics solutions at each time step (Eq. ( 7)), an iterative fixed point algorithm on equivalent contact reactions and friction forces is used with a stop criterion based on forces convergence [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. Application to an elastic layer[END_REF]. This algorithm is appropriate to the formulation of the frictional contact laws as non linear projections. The main advantage of the fixed point algorithm is that the integrator matrix remains constant at each iteration. In the simulations presented in this paper, the augmented parameters ρ n , ρ t are chosen as the smallest eigenvalue of the integrator matrix condensed on the contact degrees of freedom [START_REF] Khenous | Hybrid discretization of the Signorini problem with Coulomb friction[END_REF].

For the computation of the transient solution, the chosen time integration method is a modified θ-method. This is a first-order scheme developed by Jean [START_REF] Jean | The non-smooth contact dynamics method[END_REF] and appropriate to unilateral contact dynamics. It notably allows to compute quasi-inelastic shocks. The detail of the scheme with an application to friction-induced vibrations can be found for instance in Loyer's work [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. Application to an elastic layer[END_REF].

Results

In this section, an application of the methodology presented in the previous sections is proposed for the case of rolling contact of a wheel and a rail with lateral creepages. The model is first described. Secondly, the model is validated in quasi-static conditions with CONTACT software. A stability analysis is then carried out. Transient calculations are finally performed with reduction and results are discussed.

Description of the model

The material behavior is assumed to be linearly elastic, isotropic and undergoing small deformations. The material data of the wheel and the rail are listed in Tab. 1.

Data

Wheel rail Young's modulus (GPa) 206.9 205 Poisson's ratio 0.288 0.3 Density (kg/m 3 ) 7800 7800

Table 1: Material data of the wheel and the rail A vertical displacement u z0 is applied at the hub of the wheel. For the discretization by finite elements, compatible meshes on the interface are considered: facing nodes of the wheel and the rail in potential contact have identical tangential coordinates x and y. The potential contact area is meshed with elements of length 1 mm as shown in Fig. 2. The dynamic behavior of the wheel is assumed to be linear and lightly damped so that a real modal basis can be used. This basis is presumed to be known in the frequency range of interest. A mode is characterized by its natural circular frequency ω i , its vertical and lateral modal amplitudes at the contact point Φ z and Φ y (assumed to be normalized with respect to mass) and its structural damping factor Ξ.

The 100 first natural frequencies and corresponding free-interface modes have been calculated up to 8000 Hz. Three kinds of mode may be distinguished: the radial modes (r,m,n), the axial modes (a,m,n) and the circumferential modes (c,m,n) where n is the number of nodal diameters and m is the number of nodal circles. This classification is necessary to obtain a convergence of eigenvalues depending on the model mesh. Chosen modal damping factors η depend on the nodal diameters [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]:

   η = 10 -3 if n = 0 η = 10 -2 if n = 1 η = 10 -4 if n ≥ 2 (23)

Rail dynamics

The rail model consists of one periodically supported rail of type UIC60 provided by ESI-Group (Fig. 4). This rail has 48m length. The space between sleepers is 60cm. However the dynamic of sleepers and ballast is neglected because of its low frequency domain.

The rail hysteresis damping is η = 0.02. The rail ends with an anechoic termination of 6 meters composed of 5 lengths L = (0.6, 0.6, 0.6, 1.2, 3.0)m which are associated variable structural damping η = 0.1|1 to avoid the return of waves. Pads that connect the rail and the sleepers are modeled by springs (Kx = Ky = 36e6, Kz = 180e6)N/m. The pad hysteresis damping is η = 1.

The dynamic response obtained by FEM is compared with the dynamic response obtained with the analytical model as Timoshenko beam [START_REF] Hamet | Railway noise: use of the Timoshenko model in rail vibration studies[END_REF]. The vertical mobilities at the center point between two sleepers and on

Quasi-static results

In this section the wheel/rail quasi-static rolling contact with lateral creepage is considered. The rolling is performed in the -x direction with V = 10 m/s. An vertical displacement u z0 = 0.35 mm is applied leading For an imposed lateral creepage ∆V y /V = 0.2%, the longitudinal distributions of the normal contact reactions and frictional reaction on the center line of the contact zone (y = 0) are presented in Fig. 7. As expected, a stick zone occurs at the leading edge of the contact and a slip zone occurs at the trailing edge of the contact. The comparison of the results obtained with the proposed full FE method and the results provided by CONTACT software shows a good agreement in the slip zone and some differences for the tangential stresses in the stick zone.

Stability results

Stability analysis is carried out in case of full sliding using a higher creepage ∆V y /V = 1%. In order to solve the non symmetric eigenvalue problem, the reduction strategy presented in section 4.3 ( contact static approximation) is performed. 4 unstable modes whose real parts are positive are obtained. The frequency of these unstable complex modes are 334, 918, 1670 and 3417 Hz respectively. These frequencies obtained get close to those observed experimentally. Figs. 9 to 12 show the unstable modes shape. There are 3 axial wheel modes of (2,3,4) nodal diameters and one wheel mode (a,1,6).

Transient results

Transient results corresponding to the unstable case founded in the previous section are determined using a numerical time integration from given initial conditions. In all the following results, the integration starts from the equilibrium i.e. the initial displacements are the displacements obtained from the quasi-static solution and the initial velocities are null. The parameters V x = 10 m/s, µ = 0.3 and ∆V y /V = 1% are used. The time step for the integration is ∆t = 1µs. The time series of the lateral velocity of a point outside the contact zone is presented in the Figure 15a. The spectrogram in Fig. 15b allows to observe which frequencies are present in the solution and when they appear. At the beginning, there are two major frequencies: 334 and 919 Hz which get very close to the natural frequencies of to the two premier unstable modes provided by stability analysis. In the stationary step where nonlinear events happen, the dominant frequency is f 0 = 919 Hz and its harmonic frequencies f k = kf 0 appear.

Conclusion

In this paper, a method is presented for the modeling and analysis of the high-frequency friction-induced instabilities of wheel/rail frictional rolling contact in time and frequency domains. A full finite element formulation around the stationary position in an Eulerian reference frame is derived with a fine discretization of the contact surface combined with unilateral and Coulomb friction laws with constant friction coefficient. Both stability analysis and transient calculation are performed not only to determine unstable modes and frequencies but also to determine the amplitude and the full spectrum of vibrations.

In order to reduce the computational effort, a reduction strategy is proposed for both domains. The technique consists in simply adding a residual static contact flexibility to the free-interface normal modes when solving the frictional contact equations (contact static approximation).

The method is tested in the case of wheel/rail frictional rolling contact. The obtained results are first compared to analytical results and Kalkers theory in the quasi-static case. Application to narrow curves configuration shows that self-excited vibrations of the wheel / rail curve contact could occur due to mode coupling instability (with constant friction coefficient). The results of transient dynamics exhibit are consistent with the stability analysis and exhibit the corresponding localized stick/slip oscillations in the contact zone.

Figure 1 :

 1 Figure 1: Coordinate system, kinetic variables and applied forces

Figure 2 : 6 . 1 . 1

 2611 Figure 2: FE model of the wheel with fine mesh on the contact zone

Figure 3 :

 3 Figure 3: a: Mode axial of 136.4 Hz (m = 0, n = 1), b:Mode radial of 3281 Hz (m = 1, n = 5) et c: Mode axial of 3417 Hz (m = 0, n = 6)

Figure 4 :Figure 5 :

 45 Figure 4: Rail model

Figure 6 :

 6 Figure 6: Vertical mobility

Figure 7 :

 7 Figure 7: Quasi-static frictional stresses on a line y = 0 for ∆V y /V = 0.2% and µ = 0.3

Figure 8 :

 8 Figure 8: Divergence rate of the complex modes

Figure 9 :

 9 Figure 9: Unstable mode shape 1 (334 Hz)

Figure 10 :

 10 Figure 10: Unstable mode shape 2 (918 Hz)

Figure 11 :

 11 Figure 11: Unstable mode shape 3 (1670 Hz)

Figure 12 :

 12 Figure 12: Unstable mode shape 4 (3417 Hz)

Fig. 13a showFigure 13 :

 13 Fig.13ashow the time series of the lateral contact resultant forces F t . The tangential resultant force increases until a pronounced stick/slip oscillation builds up as shown in Fig.13b. When the tangential contact resultant force is smaller than the traction bound µF n , a transient stick zone appears at the leading edge of the effective contact region as shown in Fig.14.

Figure 14 :Figure 15 :

 1415 Figure 14: Status of nodes in the potential contact zone: no contact zone (red), slip zone (blue) and stick zone (green)

Acknowledgements

This work was carried out within the CERVIFER project backed by the ADEME organization and the Hautde-France region. It was also supported by the LabEx CeLyA of University of Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007). The authors gratefully acknowledge the support of these institutions.


Nicolas Coste
email: nicolas.coste@st.com

Hubert Garavel
email: hubert.garavel@inria.fr

Holger Hermanns
email: hermanns@cs.uni-sb.de

Richard Hersemeule
email: richard.hersemeule@st.com

Quantitative Evaluation in Embedded System Design: Validation of Multiprocessor Multithreaded Architectures *

As levels of parallelism are becoming increasingly complex in multiprocessor architectures, GALS, and asynchronous circuits, methodologies and software tools are needed to verify their functional behavior (qualitative properties) and to predict their performance (quantitative properties). This paper presents the work currently done in the Multival project (pôle de compétitivité mondial Minalogic), in which verification and performance evaluation tools developed at INRIA and Saarland University are applied to three industrial architectures designed by Bull, CEA/Leti and STMicroelectronics.

Introduction

Asynchronous computing tends to be pervasive at many levels: grids, clusters, multiprocessor architectures, multicore processors, GALS, asynchronous logic. The absence of global synchronization (no clock or several clocks) is a major departure from the classical synchronous design approach. Due to their functional complexity, there is a need to prove correctness of asynchronous architectures and circuits, but no industrial methodology is available to date. Moreover, a high degree of concurrency may introduce latencies, thus decreasing the expected performance. Consequently, correctness of functional behavior may not be sufficient for the targeted applications: performance properties * This work has been partially funded by the French government and by Conseil Général de l'Isère as part of the pôle de compétitivité Minalogic (Multival project).

(such as latencies, throughputs, resource usage, etc.) should also be considered at design time.

Launched in December 2006 after successful collaborations between its partners [START_REF] Garavel | On combining functional verification and performance evaluation using CADP[END_REF][START_REF] Hermanns | A set of performance and dependability analysis components for CADP[END_REF][START_REF] Wodey | LOTOS code generation for model checking of STbus-based SoC: The STbus interconnect[END_REF]8], the Multival project1 aims at using formal methods for functional verification and performance evaluation of asynchronous hardware. This project focuses on three different architectures that share massive asynchronous parallelism and complex communication protocols, together with strong industrial potential:

• FAME2 (Flexible Architecture for Multiple Environments) is a CC-NUMA multiprocessor architecture developed at Bull for teraflops mainframes;

• FAUST (Flexible Architecture of Unified System for Telecom) is a Network-on-Chip (NoC) based platform designed at CEA/Leti for wireless telecom applications (4G, MIMO, etc.);

• xSTream is a multiprocessor dataflow architecture for high performance embedded multimedia streaming applications designed at STMicroelectronics.

Formal Modeling Flow in Multival

A first step in the Multival approach is to model the architectures under study using formal languages suitable for representing asynchronous behaviors. We use either the LOTOS language (ISO standard 8807) or the CHP (Communicating Hardware Processes) language that can be translated automatically into LOTOS [START_REF] Salaün | Translating hardware process algebras into standard process algebras -illustration with CHP and LOTOS[END_REF].

Multival already produced a CHP model of the asynchronous network router of FAUST, and LOTOS models for the network queues of xSTream, the higher level protocols of FAUST, the FAME2 circuits implementing routing/cache-coherency protocols, and the MPI software layer and MPI benchmark applications to be run over FAME2 mainframes. These models are structural (in a bottom-up way using composition of sub-modules) or functional (in a top-down way using successive refinements).

Functional Verification Flow in Multival

The second step is to use the CADP toolbox [START_REF] Garavel | CADP 2006: A toolbox for the construction and analysis of distributed processes[END_REF] to translate these LOTOS models into Labeled Transition Systems (LTSs), whose vertices and edges correspond respectively to the states and the internal/external communications in the model. Those LTSs, which enumerate the state space of the model, can be verified using the model checking tools (based on temporal logics) and/or the equivalence checking tools (based on bisimulations) included in CADP. To avoid state space explosion, refined approaches based on compositional verification and/or co-simulation are also used.

Up to now, several significant results have been obtained: two functional issues in xSTream have been highlighted [START_REF] Lantreibecq | Verification of the xSTream communication model using CADP[END_REF]; the FAUST NoC router has been verified formally [START_REF] Salaün | Formal verification of CHP specifications with CADP -illustration on an asynchronous Network-on-Chip[END_REF], and theoretical results on isochronous forks in asynchronous circuits have been demonstrated automatically.

Performance Evaluation Flow in Multival

Performance evaluation in the Multival project is based on the IMC (Interactive Markov Chains) formalism [START_REF] Hermanns | Interactive Markov Chains and the Quest for Quantified Quality[END_REF], which combines the concepts of concurrency theory (LO-TOS, LTSs, bisimulations, etc.) with the theory of Markov chains, and which is supported by the CADP toolset [START_REF] Garavel | On combining functional verification and performance evaluation using CADP[END_REF][START_REF] Hermanns | A set of performance and dependability analysis components for CADP[END_REF].

In this approach, the LOTOS functional models for the architectures are decorated with timing information describing the stochastic behavior. This can be done either directly, by inserting stochastic transitions into the LOTOS model, or compositionally, in the following steps: (1) localization of the relevant delays in the functional model; (2) exposition (using LOTOS gates) of the start and end of each such delay; (3) instantiation of each delay by synchronizing these LOTOS gates with an auxiliary LOTOS process expressing the delay as a phase-type distribution. The decorated model is then turned into an IMC using a compositional approach (which alternates state space generation and stochastic state space minimization). This IMC is then transformed into a Markov chain, for which the Markov solvers included in CADP can compute steady-state or time-dependent state probabilities and transition throughputs.

Using this flow, early results have been obtained in the Multival project. Bull was able to predict the latency of an MPI benchmark in different topologies, different software implementations of the MPI primitives, and different cache coherency protocols. STMicroelectronics explores this flow to predict latency, throughputs in the communication architecture, and occupancy within xSTream queues.

Conclusion

The joint study of functional verification and performance evaluation is a key point in the Multival project; so far, these two aspects have been addressed independently, using different modeling languages and software tools, thus increasing costs and delays. Although formal modeling of asynchronous hardware differs significantly from the standard synchronous design methodologies, and despite its steep learning curve, at the current development stage of Multival, LOTOS was sufficient to model the three architectures. These functional models are advanced enough to enable formal verification using CADP. Compositional verification (which requires expertise) is used to address state space explosion issues. Work on performance evaluation is currently exploring two issues: new algorithms to handle nondeterminism (currently not accepted by the Markov solvers of CADP), and representations of fixed-time delays, for which there is a space-accuracy tradeoff when approximating them in the IMC formalism.

http://www.inrialpes.fr/vasy/multival


Khadidja Meguelati
email: khadidja.meguelati@inria.fr

Benedicte Fontez
email: benedicte.fontez@supagro.fr

Nadine Hilgert
email: nadine.hilgert@inra.fr

Florent Masseglia
email: florent.masseglia@inria.fr

Dirichlet Process Mixture Models made Scalable and Effective by means of Massive Distribution

Keywords: Dirichlet Process Mixture Model, Clustering, Parallelism ACM Reference Format:

Clustering with accurate results have become a topic of high interest. Dirichlet Process Mixture (DPM) is a model used for clustering with the advantage of discovering the number of clusters automatically and offering nice properties like, e.g., its potential convergence to the actual clusters in the data. These advantages come at the price of prohibitive response times, which impairs its adoption and makes centralized DPM approaches inefficient. We propose DC-DPM, a parallel clustering solution that gracefully scales to millions of data points while remaining DPM compliant, which is the challenge of distributing this process. Our experiments, on both synthetic and real world data, illustrate the high performance of our approach on millions of data points. The centralized algorithm does not scale and has its limit on 100K data points, where it needs more than 7 hours. In this case, our approach needs less than 30 seconds.

In this case, data to be considered include data on plants and crop images, like the one illustrated by Figure 1, showing a view of a Durum crop. Automatic identification, from such images, of leaves, soil, and distinguishing plants from foreground, are of high value for experts since they provide the fundamental information used for popular supervised methods in the domain [START_REF] Ma | Machine learning for Big Data analytics in plants[END_REF]. One of the main difficulties, for clustering, is the fact that we don't know, in advance, the number of clusters to be discovered. To help performing cluster analysis, despite the unknown tackled number of clusters, statistics advocate for:

(1) Setting a number of clustering runs, with varying value of K, and selecting the one that minimizes a goodness of fit criteria. It may be a quadratic risk or the Residual Mean Squared Error of Prediction (RMSEP) [START_REF] James | An introduction to statistical learning[END_REF]. This approach needs the implementation of a cross-validation algorithm [START_REF] James | An introduction to statistical learning[END_REF]. The clustering approach in this case, may be a mixture model with an Expectation-Maximization (EM) algorithm [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], or K-means [START_REF] James | An introduction to statistical learning[END_REF], for instance. (2) Making a hierarchical clustering and then cut off the tree at a given depth, usually decided by the end-user. Different approaches for pruning with advantages and drawbacks exist, see [START_REF] James | An introduction to statistical learning[END_REF]. (3) Using a Dirichlet Process Mixture (DPM) which automatically detects the number of clusters [START_REF] Michael | Estimating normal means with a Dirichlet process prior[END_REF]. In this work, we focus on the DPM approach since it allows estimating the number of clusters and assigning observations to clusters, in the same process. Furthermore, its implementation is quite straightforward in a Bayesian framework. Such properties of DPM make it a very appealing solution for many use-cases. Unfortunately, DPM is highly time consuming. Consequently, several attempts have been done to make it distributed [START_REF] Lovell | Parallel markov chain monte carlo for dirichlet process mixtures[END_REF][START_REF] Wang | Scalable Estimation of Dirichlet Process Mixture Models on Distributed Data[END_REF][START_REF] Williamson | Parallel Markov chain Monte Carlo for nonparametric mixture models[END_REF]. However, while being effectively distributed, these approaches usually suffer from convergence issues (imbalanced data distribution on computing nodes) [START_REF] Gal | Pitfalls in the use of parallel inference for the Dirichlet process[END_REF][START_REF] Lovell | Parallel markov chain monte carlo for dirichlet process mixtures[END_REF][START_REF] Williamson | Parallel Markov chain Monte Carlo for nonparametric mixture models[END_REF] or do not fully benefit from DPM properties [START_REF] Wang | Scalable Estimation of Dirichlet Process Mixture Models on Distributed Data[END_REF] (see our discussion in Section 3). Furthermore, making DPM parallel is not straightforward since it must compare each record to the set of existing clusters, a highly repeated number of times. That impairs the global performances of the approach in parallel, since comparing all the records to all the clusters would call for a high number of communications and make the process impracticable. Our goal is to propose a parallel DPM approach that fully exploits parallel architectures for better performances and offers meaningful results. Our main contribution is to keep consistency of clusters among worker nodes, and between the worker and the master nodes with regards to DPM properties. Our motivating example comes from the biology use-case described above, where the processing time on one image may go up to several days in a centralized environment. These performances in response time are the main reason why DPM is not used in the domain. The results are very informative for experts, but the processing times are prohibitive. In this case, parallelization is an appealing solution but it has to guarantee that results remain as informative as the ones targeted by a centralized run. We propose DC-DPM (Distributed Clustering by Dirichlet Process Mixture), a distributed DPM algorithm that allows each node to have a view on the local results of all the other nodes, while avoiding exhaustive data exchanges. The main novelty of our work is to propose a model and its estimation at the master level by exploiting the sufficient statistics from the workers, in a DPM compliant approach. Our solution takes advantage of the computing power of distributed systems by using parallel frameworks such as MapReduce or Spark [START_REF] Zaharia | Spark: Cluster Computing with Working Sets[END_REF]. Our DC-DPM solution distributes the Dirichlet Process by identifying local clusters on the workers and synchronizing these clusters on the master. These clusters are then communicated as a basis among workers for local clustering consistency. We modified the Dirichlet Process to consider this basis in each worker. By iterating this process we seek global consistency of DPM in a distributed environment. Our experiments, using real and synthetic datasets, illustrate both the high efficiency and linear scalability of our approach. We report significant gains in response time, compared to centralized DPM approaches, with processing times of a few minutes, compared to several days in the centralized case. The paper is organized as follows. In Section 2 we state the problem and give the necessary background on Dirichlet Process Mixture. In Section 3 we discuss related work and in Section 4, we describe the details of our distributed solution for clutering by means of Dirichlet Process Mixture. Section 5 reports the results of our experimental evaluation to verify the efficiency and effectiveness of our approach, and Section 6 concludes.

PROBLEM DEFINITION AND BACKGROUND 2.1 Dirichlet Process Mixture Models

A Dirichlet Process (DP) is a probability distribution over distributions. In our use-case, a distribution over the image pixels could be "plant" with probability p 1 , and "not plant" with probability p 2 , with the property that p 1 +p 2 = 1. A DP generates a probability distribution G. We observe a sample θ 1 ,...,θ N from G. In our use-case, each θ i is the vector of possible pixel color values.

θ n | G iid ∼ G,n = 1,...,N G ∼ DP(α,G 0 ) G θ n

N

Where G is by construction a discrete probability distribution [START_REF] Sethuraman | A constructive definition of Dirichlet priors[END_REF] .

G = ∞ i=1 π i (v)δ ϕ i
Therefore, observed variables θ n have a non null probability of having the same value ϕ i and this allows for clustering. In our use-case, "plant" pixels will have the same color vector ϕ expressing the green value. Clustering is very sensitive to the DP parameters given by the end user. G 0 is a continuous probability distribution from which the (ϕ i ) i ∈N are initially drawn. In our use-case, G 0 gives the color probability of all possible clusters in the image.

ϕ 1 ,...,ϕ i ,... ∼G 0 α is a scale parameter (α > 0) which tunes the probability weights

π i (v). V 1 ,...,V i ,... ∼ Beta(1,α) π i (v) =v i i-1 j=1 (1-v j )
α tunes indirectly the probability mass function for k N , the number of unique values (namely ϕ i ) in a sample of size N [START_REF] Antoniak | Mixtures of Dirichlet Processes with Applications to Bayesian Nonparametric Problems[END_REF].

p(k N ) = S N ,k N N !α k N Γ(α) Γ(α +N ) (1)
where S N ,k N is the unsigned Stirling number of the first kind.

With a Dirichlet Process Mixture we observe the sample y 1 ,...,y N from a mixture of distributions F (θ n ). In our use-case, we assume that colors are observed with a noise distributed according to F . The mixture is controlled by a DP on the parameters θ n .

y n ∼ F (θ n ),n = 1,...,N θ n ∼ G G ∼ DP(α 0 ,G 0 )
In a Bayesian framework, the estimation of θ n is done on the posterior: P(θ 1 ,...,θ N |y 1 ,...,y N ). Instead of this representation, another parameterization is used to speed up computation of the posterior:

P(ϕ c 1 ,...,ϕ c N |y 1 ,...,y N )
Where θ n = ϕ c n , c n is the cluster label of observation n, and ϕ c n is the unique value of the θ n belonging to the same cluster. The Gibbs algorithm [START_REF] Gelman | Bayesian Data Analysis[END_REF] samples the cluster labels c 1 ,...,c N and next the cluster parameters (here ϕ c , for all c ∈ {1,...,K } where K designs the number of cluster label values) instead of θ 1 ,...,θ N . Several versions of a Gibbs sampling are proposed by Neal in [START_REF] Radford | Markov chain sampling methods for Dirichlet process mixture models[END_REF] to simulate values from the posterior. The principle is to repeat the following loops at least until convergence to the posterior:

α π c n G 0 ϕ k y n N K N Dirichlet Process
(1) Cluster assignment, for n = 1,...,N

• Remove observation n from its cluster. Check if the cluster is empty, if yes then remove the cluster and ϕ c n from the list {ϕ} of all possible values. • Draw c n from

P (c n =c | {c j } j n ,y n ,{ϕ}) ∝ #(c) N -1+α F (y n | ϕ c ) existing cluster α N -1+α ∫ F (y n | ϕ)dG 0 (ϕ) new cluster
Where #(c) designs the number of observations affected to cluster c (after removing observation n from the sample).

• If c designs a new cluster, need to draw ϕ c from P(ϕ |y n ) ∝ F (y n | ϕ)G 0 (ϕ) • End loops (2) Update of {ϕ},
• draw ϕ c from the posterior distribution of cluster c, P(ϕ | {y} c ) (which is proportional to the product of the prior G 0 and the likelihood of all observations affected to cluster c).

When distribution F and G 0 are conjugates, ϕ can be integrated out from the Gibbs sampling which becomes time-efficient (no need to update {ϕ}). Then

P (c n =c | {c j } j n ,y n ,{ϕ}) ∝ #(c) N -1+α ∫ F (y n | ϕ)dP(ϕ | {y} c ) existing cluster α N -1+α ∫ F (y n | ϕ)dG 0 (ϕ) new cluster

Massive Distribution and Spark

Clustering via Dirichlet Process Mixture based on Gibbs Sampling is unable to scale to large datasets due to its high computational costs associated with Bayesian inference. For this reason, we aim to implement a parallel algorithm for DPM clustering in a massively distributed environment called Spark which is a parallel programming framework aiming to efficiently process large datasets. This programming model can perform analytics with in-memory techniques to overcome disk bottlenecks. Similar to MapReduce [START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF], Spark can be deployed on the Hadoop Distributed File System (HDFS) [START_REF] Shafer | The Hadoop distributed filesystem: Balancing portability and performance[END_REF].

Unlike traditional in-memory systems, the main feature of Spark is its distributed memory abstraction, called resilient distributed datasets (RDD), that is an efficient and fault-tolerant abstraction for distributing data in a cluster. With RDD, the data can be easily persisted in main memory as well as on the hard drive. Spark is designed to support the execution of iterative algorithms.

To execute a Spark job, we need a master node to coordinate job execution, and some worker nodes to execute a parallel operation. These parallel operations are summarized to two types: (i) Transformations: to create a new RDD from an existing one (e.g., Map, MapToPair, MapPartition, FlatMap); and (ii) Actions: to return a final value to the user (e.g., Reduce, Aggregate or Count).

Problem Definition

The problem we address is as follows. Given a (potentially big) dataset of records find, by means of process performed in parallel, a partition of the dataset into disjoint subsets called clusters, such that:

• Similar records are assigned to the same cluster.

• Dissimilar records are assigned to different clusters.

• The union of the clusters is the original dataset.

RELATED WORK

We set our work in the context of parallel clustering. Previous works for distributed algorithms of unsupervised clustering already exist. Ene et al. [START_REF] Ene | Fast clustering using MapReduce[END_REF] gave a MapReduce algorithm for the k-center and k-median problems. Both algorithms use Iterative-Sample as a subprocedure to get a substantially smaller subset of points that represents all of the points well. To achieve this, they perform an iterative-Sample. However these algorithms require the number of clusters k to be specified in advance, which is considered as one of the most dificult problems to solve in data clustering. Debatty et al. [START_REF] Debatty | Determining the k in k-means with MapReduce[END_REF] proposed a MapReduce implementation of G-means which is an iterative algorithm that uses Anderson Darling test to verify if a subset of data follows a Gaussian distribution. However this algorithm overestimates the number of clusters, and then requires a post-processing step to merge clusters.

In our work, we focused on algorithms inspired by the DPM. Lovell et al. [START_REF] Lovell | Parallel markov chain monte carlo for dirichlet process mixtures[END_REF] and Williamson et al. [START_REF] Williamson | Parallel Markov chain Monte Carlo for nonparametric mixture models[END_REF] has suggested an alternative parametrisation for the Dirichlet process in order to derive non-approximate parallel MCMC inference for it, these approaches are criticized by Gal and Ghahramani in [START_REF] Gal | Pitfalls in the use of parallel inference for the Dirichlet process[END_REF], these latter showed that the approaches suggested are impractical due to an extremely imbalanced distribution of the data, and gave directions for future research like the development of better approximate parallel inference. The main idea when data is distributed is to perform a DPM in each worker. The issues are then to share information between workers, and to synchronize and update clusters arising from workers at the master level. For synchronization, the main challenge is a problem of identification and of label switching of clusters. In this context we can use a relabelling algorithm like for example the one proposed by Stephens [START_REF] Jasra | Markov chain Monte Carlo methods and the label switching problem in Bayesian mixture modeling[END_REF] for mixture models. For parallel Latent Dirichlet Allocation (LDA) and Hierarchical Dirichlet Process (HDP), Newman et al. [START_REF] Newman | Distributed algorithms for topic models[END_REF] suggested to measure distance between clusters and then proposed a greedy matching. Wang and Lin [START_REF] Wang | Scalable Estimation of Dirichlet Process Mixture Models on Distributed Data[END_REF] gave a detailed review of literature and recent advanced in this topic before giving a new proposal. They proposed to use a stepwise hierarchical classification at the master level with half chance for split or merge at each step. They began with a full model considering all clusters from all workers as different components of the model. Their algorithm uses the standard Bayes Factor to compare nested models and choose the best split or merge. In conclusion, at the master level, the proposed algorithms diverge from a DPM-classifier and are not a scalable estimations of a DPM. Moreover, Wang and Lin [START_REF] Wang | Scalable Estimation of Dirichlet Process Mixture Models on Distributed Data[END_REF] used a fixed value for the scale parameter (α) in their implementation of the DPM at the workers level. The number of final clusters is related to this value (see equation 1). Authors like Miller and Harrison [START_REF] Jeffrey | Inconsistency of Pitman-Yor process mixtures for the number of components[END_REF] have demonstrated the inconsistency for the number of components of a DPM model with fixed α value. If the number of components identified at the worker level is underestimated, then the number of clusters at the master level might be underestimated. The reverse will increase considerably the running time at the master level.

In our approach, we suggest to keep to a DPM algorithm as much as possible, even at the master level, to be close to the good properties of a DPM-classifier, despite the fact that data is distributed. We also suggest a modification of the DPM model to share information among workers. In this way we expect to improve our clustering (better estimation) and suppress label switching. Finally, we do not fix a value to α but allow a different estimation in each worker to add flexibility to our model.

Furthermore [START_REF] Wang | Scalable Estimation of Dirichlet Process Mixture Models on Distributed Data[END_REF] is restricted to specific cases where noise in the observations follows the conjugate distribution of the cluster centers distribution. For example, a Gaussian noise imposes a Gaussian distribution of the centers. Therefore, this method is not suited for centers having positive values only.

Our goal is to work on any data, even with exclusively positive centers.

DC-DPM: DISTRIBUTED CLUSTERING VIA DPM

In this section, we present a novel parallel clustering approach called DC-DPM, adapted for independent data. Parallelization calls for particular attention to two main issues. The first one is the load balance between computing nodes. In our approach we distribute data evenly across the different nodes, and there is no data exchange between nodes during the processing. The second issue is the cost of communications. In order to be efficient, nodes send and receive as few information as possible by performing many iterations of Gibbs Sampling independently in each worker before synchronizing the global state at the master level and only communicating sufficient statistics between workers and master. The challenge of using sufficient statistics, in a distributed environment, is to remain in the DPM approach at all steps, including the synchronization between the worker and master nodes. The novelty of our approach is to approximate the DPM model even at the master level when local data is replaced by sufficient statistics between iterations.

Architecture and Distributed Algorithm Data Distribution

Data is evenly distributed on the computing nodes when the process starts. This is a mere, sequential, distribution, that splits the dataset into equal sized partitions.

Worker

This level handles the innovation parts of DPM (detection of new clusters) and the individual cluster assignment in each worker. The updates of the cluster labels in worker j depend on sample size proportions of the distributed data:

P(c n,j =c | c n,j ,y n,j ,{ϕ}) ∝ #(c) j N j -1+α F (y n,j ,ϕ c ),c = 1,...,K α N j -1+α ∫ F (y n,j ,ϕ)dG 0 (ϕ) new
As the clusters are not known at the beginning, we cannot ensure that the sample size proportions of each cluster will be respected in each worker. If the data were uniformly distributed, each cluster would have only, in average, the same weight/proportion in all workers. Therefore we added a modification of the update, as can be seen in Algorithm 1.

Algorithm 1 DPM at worker level for each data y n do Draw c n,j from P(c n,j =c | {c l ,j } l n ,y n,j ,{ϕ},{w },α j ) ∝ #(c)+α j w c N j -1+α j F (y n,j ,ϕ c ),c = 1,...,K

α j w u N j -1+α j ∫ F (y n,j ,ϕ)dG 0 (ϕ) new Update of α j Draw ϕ c for new clusters
where the weight w c is the proportion of observations from cluster c evaluated on the whole dataset and w u the proportion of non affected observations (awaiting the creation, innovation, discover of their real clusters). Therefore, these parameters are updated at the master level during the synchronization. Now, the scale parameter α j can be viewed as a tuning parameter between local (worker) and global (master) proportions. Following [START_REF] Michael | Bayesian density estimation and inference using mixtures[END_REF] we use an inverse gamma prior to infer this parameter.

This modification of the update implies a slightly modified DPM in each worker j :

y n,j ∼ F (θ n,j ) θ n,j ∼ G j G j ∼ DP(α j ,G) α j ∼ IG(a,b) G = K c=1 w c δ ϕ c +w u G 0 , with w u + K c=1 w c = 1

Master

This level handles the final individual assignment in the master node and therefore the final common number of clusters K. The master gets from each worker the following input: sample size of cluster k in worker j (n j,k ), cluster parameter values-sufficient statistics, individual predictive value (traditionally/usually the cluster mean value in the worker: ŷn,j = ȳj,c n, j =k ). At the master level, the observations are assigned by clusters. A cluster corresponds to a set of individuals belonging to the same cluster of the same worker. Each cluster has a representative or individual predictive value which is used to perform the end of the Gibbs sampling at the master level:

P(c n,j =c | c n,j ) ∝ #(c) N -1+γ F ( ŷn,j ,ϕ c ),c = 1,...,K γ N -1+γ ∫ F ( ŷn,j ,ϕ)dG 0 (ϕ) new
Working at an individual level implies a slow Gibbs sampling with poor mixing [START_REF] Gonzalez | Parallel gibbs sampling: From colored fields to thin junction trees[END_REF]. So, we suggest an update by clusters. In this view, we denote z j,k the master label of the cluster k in worker j. To take into account the worker information ({ϕ

P(z j,k =c | z j,k ) ∝ #(c) N -1+γ F ( ȳj,k ,ϕ c ),c = 1,...,K γ N -1+γ ∫ F ( ȳj,k ,ϕ)dG(ϕ | ϕ workerj k )
The labels {c n,j } of all the observations in the cluster k of worker j are then assigned to the master label z j,k . Next, the cluster parameters ({ϕ c } c=1,...,K ) are updated from the posterior computed on the whole dataset. We assume that we don't need all the data but only sufficient statistics from all clusters from all workers to compute the posterior. This assumption is straightforward for many distributions, as the exponential family [START_REF]Exponential Family[END_REF]. Last, the synchronization of the workers is done through the definition of G using the updated parameters ({ϕ c } c=1,...,K ) and with weights drawn from a Dirichlet distribution Dir (n 1 ,...,n K ,γ ). The end user parameters of this Dirichlet distribution are updated at the master level from the whole dataset. The size n k is the sum of all observations having label k at the end of the master Gibbs sampling.

(w 1 ,...,w K ,w u ) ∼ Dir (n 1 ,...,n K ,γ ) γ ∼ IG(c,d)

By doing so, we do not have to consider label switching. Clusters are explicitly defined at the master level and parameter values are not updated in the worker. At the worker level, only innovation (creation of new clusters) is implemented. This is summarized by Algorithm 2.

Algorithm 2 DPM at master level for each (j,k) do Draw z j,k from P(z j,k =c | {c} j,k , ȳj,k ,{ϕ},γ ) ∝

#(c) N -1+γ F ( ȳj,k ,ϕ c ),c = 1,...,K γ N -1+γ ∫ F ( ȳj,k ,ϕ)dG(ϕ | ϕ workerj k )
Update of ϕ and (w 1 ,...,w K ,w u )

The workflow of our DC-DPM approach is illustrated by Figure 2. It consists in 4 steps:

(1) Identify local new clusters in the workers (2) Compute and send sufficient statistics and cluster sizes from each worker to the master (3) Synchronize and estimate cluster labels from sufficient statistics (4) Send updated cluster parameters and cluster sizes from master to workers Our first proposition concerns the synchronization and estimation of the DPM. It is done with a Gibbs sampling conditionally on the sufficient statistics instead of the whole dataset/individual observations. Our second proposition is a construction of a shared prior distribution updated at the master level and send to the workers' DPM. This distribution reflects the information/results collected from all workers and synchronized at the master.

The likelihood for one observationy i from the Exponential family is:

F (y i | η) =h(y i )exp η T ψ (y i )-a(η)
where :

F (y 1 ,...,y n | η) = N i=1 h(y i ) exp η T N i=1 ψ (y i ) -N a(η)
Among all the distributions included in the Exponential Family, we implemented the Normal case for the experiments:

F (. | ϕ c ) = N (ϕ c ,Σ 1 )
. This choice corresponds to the simple linear model y n =ϕ c +ε n and ε n is Normally distributed N (0,Σ 1 ). In this case, the sample mean of cluster c, namely ȳc is a sufficient statistic and the posterior distribution can be conditioned only on its value:

P(ϕ | {y n } c n =c ) =P(ϕ | ȳc ) ∝ F ( ȳc | ϕ)G 0 (ϕ)
When G 0 is not a conjugate prior (e.g., a normal distribution), the posterior distribution is not a usual one but a value from this posterior can be simulated with a Metropolis Hasting (MH) within Gibbs algorithm.

When variances are known and G 0 is a conjugate prior (normal distribution N (m,Σ 1 )), there is no use of MH algorithm. The posterior is a normal distribution N (ϕ

post c = Σ(#(c)Σ -1 2 ȳc +Σ -1 1 m),Σ post c
)

where

Σ post c = (#(c)Σ -1 2 +Σ -1 1 ) -1 . The predictive posterior is a nor- mal distribution N (ϕ post c , Σ 2 + Σ post c
). In our context, Σ

EXPERIMENTS

The parallel experimental evaluation was conducted on a cluster of 32 machines, each operated by Linux, with 64 Gigabytes of main memory, Intel Xeon CPU with 8 cores and 250 Gigabytes hard disk. The centralized approach is an implementation DPM in Scala, and was executed on a single machine with the same characteristics.

The distributed algorithm we proposed is an approximation of a classic DPM, we will compare its properties to a centralized DPM implementation, on synthetic data and also in our use-case for digital agronomy. The first step of our process is a distributed K-means that sets the initial state (usually we set K to be one tenth of the dataset size).

Reproducibility : All our experiments are fully reproducible. We make our code and data available at https://github.com/khadidjaM/ DC-DPM.

In the rest of this section, we describe the datasets in Section 5.1 and our evaluation criteria in Section 5.2. Then, in Section 5.3, we measure the performances, in response time, of our approach compared to the centralized approach and also by reporting its scalability and speed-up. We evaluate the clusters obtained by DC-DPM in the case of real and synthetic dataset in Section 5.4 and Section 5.5 discusses the results and interest of our work in a real use-case of agronomy.

Datasets

We carried out our experiments on a real world and a synthetic dataset.

Our synthetic data are generated using a two-steps principle. In the first step we generate cluster centers according to a multivariate normal distribution with the same variance σ 2 1 for all dimensions. In the second step, we generate the data corresponding to each center, by using a multivariate normal distribution parameterized on the center with the same variance σ 2 2 for all dimensions. We generated a first batch of 5 datasets having size 20K, 40, 60, 80K and 100K with σ 2 1 = 1000 and σ 2 2 = 1. They represent 10 clusters. We generated a second batch of 5 datasets having size 2M, 4M, 6M, 8M and 10M with σ 2 1 = 100000 and σ 2 2 = 10. They represent 100 clusters. This type of generator is widely used in statistics, where methods are evaluated first on synthetic data before being applied on real data.

Our real data correspond to the use-case described in Section 1. The image used to test our algorithm was in RGB format. After preprocessing it contains 1,081,200 data points, described by a vector of 3 values (red, green and blue) belonging to [0,1].

Clustering Evaluation Criteria

There are two cases for evaluating the results of a clustering algorithm. Either there is a ground truth available, or there is not. In the case of an available ground truth, there are measures allowing to compare the clustering results to the reference, such as ARI, described below, for instance. This is usually exploited for experiments when one wants to check performances in a controlled environment, on synthetic data or labelled real data. In the case where there is no ground-truth (which is the usual case, because we don't know what should be discovered in real world applications of a clustering algorithm) the results may be evaluated by means of relative measures, like RSS, described below, for instance. In our experiments, we chose the following three criteria. 1. The Adjusted Rand Index (ARI): it is the corrected-for-chance version of the Rand Index [START_REF] Nguyen | Information Theoretic Measures for Clusterings Comparison: Variants, Properties, Normalization and Correction for Chance[END_REF], which is a function that measures the similarity between two data clustering results, for example between the ground truth class assignments (if known) and the clustering algorithm assignments. ARI values are in the range [-1,1] with a best value of 1. 2. The residual sum of squares (RSS): it is a measure of how well the centroids (means) represent the members of their clusters. It is the squared distance of each data from its centroid summed over all vectors. In the univariate case, the RSS value divided by the number of observations gives the value of the Mean Squared Error (MSE), an estimator of the residual variance. In multivariate dataset with independent variables, the RSS value divided by the number of observations gives an estimator of the sum of the variable variances. This sum represents its lower bound and also the best value to be observed in the clustering of synthetic data. To simplify, we give in the following the result of the RSS value divided by the number of data N and the variance. Therefore the lower bound is known and should be equal to the number of variables (for example 2 for our synthetic data). 3. K, the number of discovered clusters.

Response Time

In this section we measure the clustering time in DC-DPM and compare it to the centralized approach. Figure 3 reports the response times of DC-DPM and the centralized approach on our synthetic data, limited to 100K data points. Actually, the centralized approach does not scale and would take several days for larger datasets.The distributed approach is run on a cluster of 8 nodes. The results reported by Figure 3 are in logarithmic scale. The clustering time increases with the number of data points for all approaches. This time is much lower in the case of DC-DPM, than the centralized approach. On 8 machines (64 cores) and for a dataset of 100K data points, DC-DPM performs the clustering in 24 seconds, while the centralized approach needs more than 7 hours on a single machine.

Figure 4 reports an extended view on the clustering time, only for DC-DPM, and with a dataset having up to 10 million data points from the synthetic dataset. DC-DPM is run on a cluster of 16 machines. The running time increases with the number of data points. Let us note that the centralized approach does not scale and cannot execute on such dataset size. DC-DPM enjoys linear scalability with the dataset size.

Figures 5 and6 illustrate the parallel speed-up of our approach on 2M data points from the synthetic dataset and on more than 1 million data points obtained after preprocessing the image of our use-case. The results show optimal or near optimal gain. In Figure 6 we observe that the response time for 2 nodes is more than twice the response time for 4 nodes. That is unexpected when measuring a speed-up. However, the response times of our approach are very fast (a few minutes) and do not consider the time it takes for Spark to load-up, before running DC-DPM. The slight difference between an optimal speedup and the results reported in Figure 6 are due to that loading time.

Clustering Evaluation

In the following experiments, we evaluate the clustering performance of DC-DPM and compare it to the centralized approach.

Table 1 reports 1) the ARI value computed between the clustering obtained and the ground truth, 2) the RSS value divided by the number of data N and variance (σ 2 2 ), and 3) the number of clusters, obtained with DC-DPM and the centralized approach on our synthetic data while increasing the dataset size. The DC-DPM is run on a cluster of 8 nodes. DC-DPM performs as well as the centralized approach, there is a small gap in RSS values which is negligible compared to the gained time.

Table 2 reports an extended view on the ARI value, and the RSS value divided by the number of data N and by the variance (σ 2 2 ), and number of clusters obtained with DC-DPM, while increasing dataset size (up to 10 million data points). DC-DPM is run on a cluster of 16 machines. The performance keeps showing the maximum possible accuracy, even with a large number of data points. Figure 7 gives a visual representation of our 4M data points synthetic dataset. Each cluster is assigned to a color. Our goal is to retrieve these clusters. Figure 8 represents the results obtained by our approach on the data of Figure 7, with 16 nodes. Each cluster is assigned to a different color. shows the performance of our approach with almost perfect results where the discovered clusters are the same as the actual ones from the data. This is confirmed by Table 2, line 2.

Use-case

Phenotyping and precision agriculture use more and more information from sensors and drones, like aerial images, leading to the emerging domain of digital agriculture (see for example http://ec.europa. eu/research/participants/portal/desktop/en/opportunities/h2020/topics/ dt-rur-12-2018.html). An important challenge, in this context, is to be able to distinguish clusters of plants: status (normal, hydric stress, disease,...) or species for example. Clustering, applied to images, is a key in this domain.

We want to discover clusters in the image presented in Section 1 and transformed as described in Section 5.1. We set σ 2 1 = 1 because our data is in the range of [0,1]. For both parameter values of σ 2 2 = 0.01 and σ 2 2 = 0.0025, the clusters are extracted in approximately 3 minutes with DC-DPM running in parallel on 16 computing nodes. This is confirmed by Figure 6. The centralized approach does not scale on this data and we could not obtain results.

The number of clusters depends on the value of the variance error. A value of σ 2 2 = 0.01 gave a rough clustering with only K = 3 clusters. Those clusters identified the brightness in the RGB image (see figure 9). A lower value of σ 2 2 = 0.0025 gave a clustering with k = 12 clusters, which is enough to reconstruct the image (see figure 10). Depending on the aim of the clustering, different types of wavelength or data must be used for identification. The accuracy of the clustering (number of clusters) relies on the variance value (σ 2 2 ). Clustering is used to detect structures in the data (genetic, population, status) before processing. This group detection allows reducing data dimension and bias in further prediction analysis.

DC-DPM was compared to the centralized DPM on a part of the the Durum image in the experimental field. RGB Image with DC-DPM (top, 12 clusters) and centralized DPM (bottom, 17 clusters), σ 2 2 = 0.0025. The results were quite similar as shown in figure 11. The impact of σ 2 on the number of clusters varies for centralized and distributed approaches and may be adjusted by the end-user.

CONCLUSION

We proposed DC-DPM, a novel and efficient parallel solution to perform clustering via DPM on millions of data points. We evaluated the performance of our solution over real world and synthetic datasets. The experimental results illustrate the excellent performance of DC-DPM (e.g., a clustering time of less than 30 seconds for 100K data points, while the centralized algorithm needs several hours). The results also illustrate the high performance of our approach with results that are comparable to the ones of the centralized version. Overall, the experimental results show that by using our parallel techniques, the clustering of very large volumes of data can now be done in small execution times, which are impossible to achieve using the centralized DPM approach. A nice perspective of our approach is to open fundamental research tracks such as DPM clustering on complex data like, e.g. time series.

Figure 1 :

 1 Figure 1: Durum in an experimental field. RGB image.

  workerj k }), we replace the Dirichlet Process Mixture Models made Scalable and Effective by means of Massive DistributionSAC '19, April 8-12, 2019, Limassol, Cyprus prior predictive distribution ( ∫ F (y n,j ,ϕ)dG 0 (ϕ)) by a posterior predictive distribution. Eventually, we use the cluster mean value ( ȳj,k ) as an individual predictive value:

Figure 2 :

 2 Figure 2: Diagram/workflow of the DC -DPM

  the mean value ϕ post c was replaced by an individual drawn from the posterior.

Figure 3 :

 3 Figure 3: Logarithmic scale. Response time of DC-DPM and centralized DPM as a function of the dataset size.

Figure 4 :

 4 Figure 4: Response time (minutes) of DC-DPM as a function of the dataset size.

Figure 5 :

 5 Figure 5: Clustering time as a function of the number of computing nodes on the synthetic data.

Figure 6 :

 6 Figure 6: Clustering time as a function of the number of computing nodes on the image of our use-case.

Figure 7 :

 7 Figure 7: Visual representation of the synthetic dataset.

Figure 8 :

 8 Figure 8: Visual representation of the results obtained by DC-DPM

Figure 9 :

 9 Figure 9: Clustering of the Durum image by DC-DPM with σ 2 = 0.01.

Figure 10 :

 10 Figure 10: Clustering of the Durum image by DC-DPM with σ 2 2 = 0.0025.

Figure 11 :

 11 Figure 11: Clustering of a part of the Durum image by DC-DPM (top) and Centralized DPM (bottom).

  Mixture Models made Scalable and Effective by means of Massive DistributionSAC '19, April 8-12, 2019, Limassol, Cyprus

Table 1 :

 1 Clustering evaluation criteria obtained with the centralized DPM and with DC-DPM.

	Centralized DPM		DC-DPM
	ARI	RSS N ×σ 2 2	Clusters ARI	RSS N ×σ 2 2	Clusters
	20K 1.00	2.01	10	1.00	2.04	10
	40K 1.00	2.00	10	1.00	2.03	10
	60K 1.00	2.00	10	1.00	2.02	10
	80K 1.00	2.00	10	1.00	2.01	10
	100K 1.00	2.00	10	1.00	2.02	10

Table 2 :

 2 Clustering evaluation criteria obtained by DC-DPM .

	ARI RSS/(N*σ 2 2 ) Clusters
	2M 1.00	2.00	102
	4M 1.00	2.00	100
	6M 1.00	2.00	100
	8M 1.00	2.02	99
	10M 1.00	2.10	101

ACKNOWLEDGEMENTS

The research leading to these results has received funds from the European Union's Horizon 2020 Framework Programme for Research and Innovation, under grant agreement No. 732051.