Yann Brenier

A DOUBLE LARGE DEVIATION PRINCIPLE FOR MONGE-AMP ÈRE GRAVITATION

Monge-Ampère gravitation is a nonlinear modification of classical Newtonian gravitation, when the Monge-Ampère equation substitutes for the Poisson equation. We establish, through two applications of the large deviation principle, that the MA gravitation for a finite number of particles can be reduced, through a double application of the large deviation principle, to the simplest possible stochastic model: a collection of independent Brownian motions with vanishing noise.

Introduction

The purpose of this paper is twofold. We first want to make a short presentation of the Monge-Ampère gravitational (MAG) model, which has been introduced in [START_REF] Brenier | A modified least action principle allowing mass concentrations for the early universe reconstruction problem[END_REF], in close connection with earlier works such as [START_REF] Brenier | Derivation of the Euler equations from a caricature of Coulomb interaction[END_REF][START_REF] Brenier | A geometric approximation to the Euler equations: The Vlasov-Monge-Ampère equation[END_REF][START_REF] Loeper | Quasi-neutral limit of the Euler-Poisson and Euler-Monge-Ampère systems[END_REF][START_REF] Cullen | The semigeostrophic equations discretized in reference and dual variables[END_REF][START_REF] Ambrosio | Hamiltonian ODE in the Wasserstein spaces of probability measures[END_REF] as well as with optimal transport theory (see [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Ambrosio | Savaré Gradient flows in metric spaces and the Wasserstein spaces of probability measures[END_REF]). This model can be seen as a nonlinear modification of the classical model of Newtonian gravitation, for which we use the fully nonlinear Monge-Ampère equation as a substitute for the linear Poisson equation to derive the gravitational potential from the density field. We also briefly compare the Monge-Ampère and the Newton gravitational models, emphasizing that in the MAG model there is an absolute control of the acceleration of the gravitating particles, no matter how concentrated they can be, in sharp contrast with Newtonian gravitation, for which some particles may runaway to infinity in finite time. However, the main purpose of this paper is about deriving the MAG model from a completely elementary microscopic model in which a finite number of particles just move as independent Brownian trajectories without any interaction. In order to get the MAG model, we need two applications of the large deviation principle (LDP) [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF][START_REF] Dembo | Large deviations techniques and applications[END_REF]. (Let us mention the interesting connection between large deviation principles, gradient flows and optimal transport theory recently made in [START_REF] Ch | From the Schrödinger problem to the Monge-Kantorovich problem[END_REF][START_REF] Mielke | On the relation between gradient flows and the largedeviation principle, with applications to Markov chains and diffusion[END_REF], which was certainly influential for us.) Through the first application of the LDP, we get a first order (in time) dynamical system based on the indistinguishability of the gravitating particles. (In a rather paradoxical way, the indistinguishability principle leads to a model where particles do interact.) The second application of the LDP enables us to lift this first order dynamical system to a second order one which is nothing but the discrete version of the MAG model. Let us emphasize that this derivation is purely formal and further investigations are clearly needed to get a complete and rigorous theory. The first section of this paper will be devoted to a presentation of both the classical Newtonian and the Monge-Ampère models, with a suitable formulation of the MAG model relying on "optimal transport" tools [START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF][START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF][START_REF] Villani | Topics in optimal transportation[END_REF]. The second section is devoted to the double application of the LDP in order to recover the discrete version of the MAG model.

1. Monge-Ampère and Newton gravitations 1.1. Classical Newtonian gravitation. To describe classical Newton gravitation, let us attach to each gravitating particle a label a, which we suppose to belong to an abstract Borel measure space (A, λ), and its position at time t, X t (a) ∈ R d (classically d = 3). For each time t, we denote by ρ t the image measure of λ by X t , defined on R d by (1)

R d γ(x)ρ t (dx) = A γ(X t (a))λ(da) ∀γ ∈ C 0 c (R d ).
For the sake of simplicity, we further assume ρ t (dx) to be Z d periodic in x and of unit mean on the unit cube [0, 1] d , so that we can see ρ t (dx) as a probability measure on the flat torus

T d = R d /Z d . (
This kind of assumption is common in computational Cosmology [START_REF] Frisch | A reconstruction of the initial conditions of the Universe by optimal mass transportation[END_REF].) We are now ready to write the Newtonian model:

(2)

d 2 X t (a) dt 2 = -∇φ t (X t (a)), △φ t = ρ t -1,
where the average density 1 has been subtracted out from the field ρ t so that the "gravitation potential" φ t is well defined from the Poisson equation as a zero mean Z d -periodic function. (This is a perfectly meaningful assumption: there should be no resulting force ∇φ t out of a uniform ρ t .) If, at time t, ρ t is just a plain probability measure on T d , then ∇φ t , obtained through the Poisson equation △φ t = ρ t -1, is merely an L p function for p < (1d -1 ) -1 on T d . This makes very dubious the meaning of ODE

d 2 X t (a) dt 2 = -∇φ t (X t (a))
in the general case. This is why Newtonian's gravitation usually splits up into two, rather unrelated theories, according to the choice of its initial conditions: 1) As A is a discrete set of N particles, then Newtonian gravitation boils down to an "N -body problem", set on the flat torus T d , which is of form

(3) d 2 X t (a) dt 2 = R(X t (a)) + b =a G(X t (a), X t (b)),
where R is a fixed function depending on the torus, and the singular kernel G(x, y), associated with the operator -∇△ -1 , behaves as (yx)/|y -x| d , as x approaches y. This "bad" singularity makes the study of the "N -body problem" very hard, in particular because of the possibility that particles may collide or runaway at infinity, in finite time. (See [START_REF] Chenciner | A remarkable periodic solution of the three-body problem in the case of equal masses[END_REF] as a recent reference.) 2) As A corresponds to a "continuum" of particles, it is fruitful to introduce, at each time t, the nonnegative measure f t defined on the "phase space" T d × R d by (4)

T d ×R d γ(x, v)f t (dx, dv) = A γ(X t (a), dX t (a) dt )λ(da), ∀γ ∈ C 0 c (T d × R d ).
A straightforward calculation shows that f must satisfy, at least formally, the "Vlasov-Poisson" system [START_REF] Rein | Collisionless kinetic equations from astrophysics-the Vlasov-Poisson system[END_REF] (5)

∂ t f t + ∇ x • (vf t ) -∇ v • (∇φ t f t ) = 0, △φ t = ρ t -1 = R d f t (dv) -1.
The meaning of this equation is unclear in the general case, because of the dubious product of ∇φ t by f t . However, it is not hard to see that 0 ≤ f t (dx, dv) ≤ Cdxdv for some positive constant C is a consistent property for equation [START_REF] Boilliat | Recent mathematical methods in nonlinear wave propagation[END_REF]. It turns out that this property is sufficient to justify the multiplication of f t by ∇φ t at least when the initial condition f 0 vanishes for large values of v, say for |v| ≥ C ′ for some other constant C ′ . Then, the Newtonian model gets well defined and, at least for d ≤ 3, global weak solutions can be shown to exist globally in time, as f 0 is given, whatever are constants C and C ′ . In addition, such solutions are unique and smooth as long as f 0 has a smooth density with respect to the Lebesgue measure.

(See [START_REF] Pfaffelmoser | Global classical solutions of the Vlasov-Poisson system in three dimensions for general initial data[END_REF][START_REF] Lions | Propagation of moments and regularity for the 3-dimensional Vlasov-Poisson System[END_REF][START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF] for such results, stated on the whole space R d rather than on T d . See also [START_REF] Rein | Collisionless kinetic equations from astrophysics-the Vlasov-Poisson system[END_REF], as well as [START_REF] Deng | Solutions of Euler-Poisson equations for gaseous stars[END_REF] for the closely related Euler-Poisson system.)

1.2. Monge-Ampère gravitation. We now introduce the Monge-Ampère gravitation (MAG) model which differs from (2) just by the substitution of the Monge-Ampère equation for the Poisson equation:

(6) d 2 X t (a) dt 2 = -∇φ t (X t (a)), det(I + D 2 φ t ) = ρ t ,
where I denotes the d × d identity matrix and φ t is solution to the Monge-Ampère equation (in a suitable sense), which is Z d -periodic with zero mean and satisfies the (weak) ellipticity condition I+ D 2 φ t (x) ≥ 0, in the sense of symmetric matrices, for every x. [Notice that φ t is unique and Lipschitz continuous (resp. smooth) as soon as ρ t (dx) has an integrable (resp. smooth and positive) density with respect to the Lebesgue measure dx, see [START_REF] Cordero-Erausquin | Sur le transport de mesures périodiques[END_REF] for example.] We see that Newtonian gravitation can be formally recovered from the MAG model just by expanding the determinant about I and retaining only the linear part:

det(I + D 2 φ t ) ∼ 1 + TraceD 2 φ t = 1 + △φ t .
Notice that, as d = 1 (which is a case of limited interest, describing gravitating "parallel pancakes"), the MAG model coincide with Newtonian gravitation. (As a consequence, our derivation of the MAG model from a double application of the large deviation principle, obtained in the second section of the present paper, is also valid for the Newtonian gravitational model in one space dimension.) However, let us emphasize that, to the best of our knowledge, MAG has never been considered by any physicist in dimension larger than one. It has been so far only considered by mathematicians (see [START_REF] Brenier | A modified least action principle allowing mass concentrations for the early universe reconstruction problem[END_REF] and [START_REF] Brenier | A geometric approximation to the Euler equations: The Vlasov-Monge-Ampère equation[END_REF][START_REF] Loeper | Quasi-neutral limit of the Euler-Poisson and Euler-Monge-Ampère systems[END_REF][START_REF] Cullen | The semigeostrophic equations discretized in reference and dual variables[END_REF][START_REF] Ambrosio | Hamiltonian ODE in the Wasserstein spaces of probability measures[END_REF] for closely related topics), mostly because of its close connection with optimal transport theory (as discussed in the next subsection). It is easy to describe, as we did for the Newtonian model, the MAG model through a "kinetic equation", the so-called "Vlasov-Monge-Ampère" (VMA) system [START_REF] Brenier | A geometric approximation to the Euler equations: The Vlasov-Monge-Ampère equation[END_REF][START_REF] Loeper | Quasi-neutral limit of the Euler-Poisson and Euler-Monge-Ampère systems[END_REF]:

(7) ∂ t f t + ∇ x • (vf t ) -∇ v • (∇φ t f t ) = 0, det(I + D 2 φ t ) = ρ t = R d f t (dv).
As for the Vlasov-Poisson system (5), the existence of global weak solutions can be shown as soon as

0 ≤ f 0 (dx, dv) ≤ Cdxdv, f 0 (T d × {v ∈ R d , |v| ≤ C ′ }) = 0,
for some positive constants C, C ′ . Existence of a unique smooth solution, but only for a short time interval, has been proven by Loeper [START_REF] Loeper | Quasi-neutral limit of the Euler-Poisson and Euler-Monge-Ampère systems[END_REF] provided f 0 (dx, dv) (resp. ρ 0 (dx)) has a smooth density with respect to dxdv (resp. dx) and is uniformly compactly supported in v ∈ R d (resp. is strictly positive). These are still limited results, due to the full non-linearity of the Monge-Ampère equation that leads to analytic difficulties, in particular when seeking for smooth solutions. However, the Monge-Ampère equation enjoys remarkable properties, closely related to the theory of optimal transportation [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF][START_REF] Villani | Topics in optimal transportation[END_REF]. This is why we are going to introduce a related formulation of the MAG model, with interesting geometric features.

1.3. The MAG model written in optimal transportation terms. For the description of the MAG model in optimal transport terms, it is convenient to discuss the model in a (slightly) different and more abstract framework. We consider a metric measure space (A, λ), made of a compact subset A of R d equipped with a Borel probability measure λ. Two typical examples are, on one side, the unit cube with the Lebesgue measure, and, on the other side, any set of N points equipped with the (normalized) counting measure. In the first case, we will speak of the "continuous" case, while, in the second case, we will speak of the "discrete" case.

(Of course many others situations could be also considered, in particular the flat torus T d as we did in the previous subsections, but we will focus on these two cases only.) We introduce the separable Hilbert space H of all λ-square-integrable (Observe, in the discrete case, when A is made of

maps from A to R d , H = L 2 (A, λ; R d ),
N distinct points A(a) ∈ R d , for a = 1, • • • , N , S just reads S = {(A(σ(1)), • • • , A(σ(N )) ∈ H, σ ∈ S N } with N !
elements, where S N denotes the group of all permutations of the first N integers.)

According to Edelstein's theorem [START_REF] Edelstein | On nearest points of sets in uniformly convex Banach spaces[END_REF][START_REF] Aubin | Mathematical methods of game and economic theory[END_REF], given a separable Hilbert space H and a closed bounded subset S, there is, in the sense of Baire, a generic set (i.e. containing a countable intersection of dense open subsets of H) of points X for which there exists a unique closest point π(X) on S. In addition π(X) is nothing but the gradient, at point X, of the Lipschitz convex function Π defined on the Hilbert space H by ( 9)

Π(X) = sup s∈S ((X, s)) - ||s|| 2 2 , π(X) = ∇Π(X) = Arg inf s∈S ||X -s|| 2 2 .
In our particular case, we can say much more in the continuous case, when A is the unit cube with λ as the Lebesgue measure, thanks to "optimal transport theory" Theorem 1.1. ( [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]) Let X ∈ H be a non degenerate map, in the sense that the image measure ρ of λ by X is absolutely continuous with respect to the Lebesgue measure on R d . Then, X has a unique closest point π(X) on S. We also have

(10) π(X) = T • X,
where T is the unique map (in the ρ a.e. sense) y ∈ R d → T (y) ∈ A such that: i) there is a convex Lipschitz function ψ : R d → R, ρ-a.e. differentiable, with:

(11) T (y) = ∇ψ(y), ρ -a.e. y ∈ R d
ii) λ is the image of ρ by T :

(12) R ⌈ γ(T (y))ρ(dy) = A γ(a)λ(da), ∀γ ∈ C 0 (R d ).
As a matter of fact, equations [START_REF] Brenier | Derivation of the Euler equations from a caricature of Coulomb interaction[END_REF][START_REF] Brenier | Hydrodynamic sturcture of the augmented Born-Infeld equations[END_REF] form a generalized formulation of the Monge-Ampère problem on R d [START_REF] Brenier | A modified least action principle allowing mass concentrations for the early universe reconstruction problem[END_REF] ρ = det(D 2 x ψ), (∇ψ)(support(ρ)) = A, D 2 ψ ≥ 0, with a unique solution ∇ψ (in the ρ-a.e. sense). (See [START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF][START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF][START_REF] Villani | Topics in optimal transportation[END_REF] for more details.) This generalized formulation of the Monge-Ampère equation allows us to write the MAG model ( 6) in a much more geometric way. Indeed, (6) just reads

d 2 X t (a) dt 2 = X t (a) -∇ψ t (X t (a)), det(D 2 ψ t ) = ρ t ,
after setting ψ t (x) = |x| 2 /2 + φ t (x), that we complete with the weak ellipticity condition D 2 ψ t ≥ 0 and the range condition ∇ψ t (support(ρ t )) = A (which substitutes for the Z d -periodicity condition we used in writing ( 6)). Assume that, at time t, X t is non degenerate (or, in other words, ρ t is absolutely continuous with respect to the Lebesgue measure on R d ). Then, using Theorem 1.1, we may write

∇ψ t • X t = π(X t ) = ∇Π(X t )
and finally obtain ( 14)

d 2 X t dt 2 = X t -π • X t = X t -∇Π • X t ,
where Π is the Lipschitz convex function defined by [START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF], π its gradient and S is the set of measure preserving maps defined by [START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF]. In the rest of this paper, we will retain [START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF][START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF][START_REF] Brenier | Reconstruction of the early universe as a convex optimization problem[END_REF] as our definition of the MAG model. 1.4. The discrete Monge-Ampère gravitational model. The geometric formulation [START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF][START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF][START_REF] Brenier | Reconstruction of the early universe as a convex optimization problem[END_REF] of the MAG model is very convenient to get its discrete version, as A is just a subset of N points A(a) in R d , for a = 1, • • • , N , equipped with the counting measure, in which case, we will speak of the "discrete MAG model with N particles". Indeed, at the discrete level, a time-dependent map X t can be seen as a sequence of N "particles", with positions X t (a), moving in R d , for a = 1, • • • , N . (Notice that X t gets degenerate at time t in case of collisions.) In this discrete setting, the MAG model [START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF][START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF][START_REF] Brenier | Reconstruction of the early universe as a convex optimization problem[END_REF] reads as the finite dimensional dynamical system:

(15) d 2 X t (a) dt 2 = X t (a) -A(σ t (a)), σ t = Arg inf σ∈SN N a=1 |X t (a) -A(σ(a))| 2 ,
where σ ∈ S N denotes the set of all permutations of the first N integers. So we see that, at the discrete level, the MAG looks both very simple and very different from the classical Newtonian N body problem (3)! (Similar systems have been previously studied in [START_REF] Brenier | Derivation of the Euler equations from a caricature of Coulomb interaction[END_REF][START_REF] Cullen | The semigeostrophic equations discretized in reference and dual variables[END_REF].)

The MAG model, as defined by [START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF][START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF][START_REF] Brenier | Reconstruction of the early universe as a convex optimization problem[END_REF], can be seen as a dynamical system in a Hilbert space H with a force term F (x) = x -∇Π(x) with a trivial linear part and the gradient of a Lipschitz convex function Π. To the best of our knowledge, there is no theory for such an ODE in infinite dimension. Indeed, we are very far from the standard Cauchy-Lipschitz setting. This is a challenging open problem.

(See the related theory developed by Ambrosio and Gangbo for some infinite dimensional hamiltonian systems [START_REF] Ambrosio | Hamiltonian ODE in the Wasserstein spaces of probability measures[END_REF]. As a matter of fact, their main example is very similar to the MAG model, written in a different way. See also [START_REF] Brenier | A geometric approximation to the Euler equations: The Vlasov-Monge-Ampère equation[END_REF][START_REF] Gangbo | Hamilton-Jacobi equations in the Wasserstein space[END_REF].) However, for the discrete MAG model with N particles, H ∼ R N d is of finite dimension. Thus, equation ( 14), which also reads [START_REF] Brenier | A geometric approximation to the Euler equations: The Vlasov-Monge-Ampère equation[END_REF], can be neatly solved globally in time in the sense of Bouchut and Ambrosio [START_REF] Bouchut | Renormalized solutions to the Vlasov equation with coefficients of bounded variation[END_REF][START_REF] Ambrosio | Transport equation and Cauchy problem for BV vector fields[END_REF], for every initial condition (X 0 , dX0 dt ), except on a negligible subset of the "phase space" H × H, of zero 2N d dimensional Lebesgue measure. (A more accurate statement can be found in [START_REF] Ambrosio | Transport equation and Cauchy problem for BV vector fields[END_REF].) Indeed, the force term F (x) = x -∇Π(x) in the right-hand side of ( 14) is a smooth perturbation of a "bounded variation" function, since Π is a Lipschitz convex function (which implies that D 2 Π can be seen as a bounded nonnegative measure valued in the convex cone of nonnegative symmetric matrices, and, therefore, that ∇Π is of bounded variation). Notice that the exceptional set of bad initial conditions in the phase space is not empty, as it follows clearly from formulation [START_REF] Brenier | A geometric approximation to the Euler equations: The Vlasov-Monge-Ampère equation[END_REF], where we see that the evolution of the particles becomes ambiguous as different particles depart from the same position with exactly the same velocity.

1.5. Monge-Ampère versus Newtonian gravitational models. It is now interesting to compare the Monge-Ampère and the Newton gravitational models. We first observe that, according to the MAG model, particles may never runaway to infinity in finite time. Indeed, from the optimal transportation formulation (8,9,14), we immediately get

| d 2 X t dt 2 -X t | ≤ R = sup
|X t | + | dX t dt | ≤ (|X 0 | + | dX 0 dt | + 1)C cosh t, ∀t ∈ R,
where C depends only on R. Of course, there is nothing similar with Newtonian gravitation. Indeed, the Poisson equation in (2) does not behave well when particles concentrate. For instance, if some particles concentrate as a delta measure at some point y at time t, then ∇φ t (x) has a singularity as bad as (xy)/|x -y| d . In this way, the MAG model has a lot of similarity with the Born-Infeld (BI) theory of the electromagnetic field [START_REF] Born | Foundations of the new field theory[END_REF][START_REF] Boilliat | Recent mathematical methods in nonlinear wave propagation[END_REF][START_REF] Brenier | Hydrodynamic sturcture of the augmented Born-Infeld equations[END_REF][START_REF] Brenier | Derivation of particle, string, and membrane motions from the Born-Infeld electromagnetism[END_REF], in which any electrostatic field is bounded by a universal constant. (Notice that the BI model, which goes back to 1934, is no longer for use in Electrodynamics, but has enjoyed a remarkable revival in String Theory since the 1990s [START_REF] Polchinski | String theory[END_REF].)

In addition, the MAG model enjoys good properties of approximations by finite sums of Dirac measures (see [START_REF] Cullen | The semigeostrophic equations discretized in reference and dual variables[END_REF] for closely related results). In sharp contrast, such discrete approximations have never been justified, to the best of our knowledge, in the case of 3D Newtonian gravitation, because of the bad singularities of the Green function for the Poisson equation [START_REF] Hauray | N-particles approximation of the Vlasov equations with singular potential[END_REF]. As a matter of fact, the treatment of point particles in classical Electrodynamics remains an outstanding open problem in both Theoretical and Mathematical Physics (see [START_REF] Feynman | The Feynman lectures in Physics[END_REF][START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF] for instance).

2. Formal derivation of the discrete Monge-Ampère gravitation model from a double application of the large deviation principle

In this second section, we provide a f ormal derivation of the discrete MAG model ( 15) from a double application of the large deviation principle [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF][START_REF] Dembo | Large deviations techniques and applications[END_REF].

2.1.

A basic model of N independent Brownian particles. We fix a positive integer N and a finite set

A of N distinct points A(a) ∈ R d , for a = 1, • • • , N .
Using the notations of subsection 1.3, we set H = R N d and ( 16)

S = {(A(σ(1)), • • • , A(σ(N )) ∈ H, σ ∈ S N },
where S N denotes the set of all possible permutations of the first N integers. Given σ 0 ∈ S N , we consider the motion of N "particles" that move in the Euclidean space R d , according to

(17) X ε t (a) = s 0 (a) + εB t (a), s 0 (a) = A(σ 0 (a)), ∀a = 1, • • • , N.
where

(t ∈ R + → B t (a) ∈ R d ) a=1,••• ,N denote N independent realizations of the "standard" (i.e. normalized) Brownian motion in R d . Fixing t * > 0 and Y * (a) ∈ R d , ∀a = 1, • • • , N,
it is easy to compute the probability that, at time t * , each particle occupies the position given by Y * :

Prob(X ε t * (a)≈Y * (a), ∀a = 1, • • • , N ) ≈ N a=1 [exp( -|Y * (a) -s 0 (a)| 2 2εt * )(2πεt * ) -d/2 ],
or, in other words

Prob(X ε t * ≈Y * ) ≈ exp( -||Y * -s 0 || 2 2εt * )(2πεt * ) -N d/2 .
Here we have denoted by | • | and || • || the Euclidean distance on respectively R d and H = (R d ) N . We have also used symbol ≈ just to make notations lighter. [What we precisely mean is: for any Borel B subset of (R d ) N , the probability that

X ε t * belongs to Y * + B is given by Y ∈Y * +B exp( -||Y -s 0 || 2 2εt * )(2πεt * ) -N d/2 dY,
but we hope that our simplified notation is acceptable. Indeed, we will use it again, without further notice.] 2.2. First application of the large deviation principle. We now want to compute the probability that X ε t * and Y * coincide, up to a permutation, property that we denote by (X ε t * ) ≈ perm (Y * ). We find

(18) Prob[(X ε t * ) ≈ perm (Y * )] ≈ 1 N ! σ∈SN exp( -||Y * • σ -s 0 || 2 2εt * )(2πεt * ) -N d/2 .
When the level of noise ε goes to zero, we immediately get

-lim ε→0 ε log Prob[(X ε t * ) ≈ perm (Y * )]≈ inf σ∈SN ||Y * • σ -s 0 || 2 2t * ,
which is a rather trivial illustration of the Laplace method and the large deviation principe (LDP) [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF]. There is a more sophisticated aspect of the LDP: as the level of noise goes to zero, the Brownian trajectories of the particles conditioned by (X ǫ t * ) ≈ perm (Y * ) behave more and more as constant speed minimizing geodesic curves. More precisely, for all t ∈ [0, t * ],

(

) X ε t ∼ ε→0 X t = s 0 + t t * (Y * • σ * -s 0 ), σ * = Arg inf σ∈SN ||Y * • σ -s 0 || 2 2t * . which implies (20) dX t dt = Y * • σ * -s 0 t * = X t * -s 0 t * = X t -s 0 t , ∀t ∈]0, t * ]. 19
Now, we are going to translate this large deviation principle into a self-consistent dynamical system for the particles. Let us first denote by π(Y ) (as in subsection 1.3) the unique closest point on S of a "generic" point

Y ∈ H = (R d ) N (21) π(Y ) = Arg inf s∈S ||Y -s|| 2 2 = ∇Π(Y ), Π(Y ) = sup s∈S ((Y, s)) - ||s|| 2 2
(where ((•, •)) denotes the inner product attached to || • || on H). Then we state:

Proposition 2.1. Equation ( 20), that we have derived from (17) (through a large deviation principle), can be written as a self-consistent ordinary differential equation for X t at least for t > 0,

t dX t dt = X t -π(X t ). (22)
Proof. To get this result, our simple but crucial observation is that, along the geodesic curve defined by [START_REF] Cordero-Erausquin | Sur le transport de mesures périodiques[END_REF], namely

X t = s 0 + t t * (Y * • σ * -s 0 ), σ * = Arg inf σ∈SN ||Y * • σ -s 0 ||, s 0 is the closest point in S not only of the end-point X t * = Y * • σ * but also of all points X t , ∀t ∈ [0, t * ].
[Then, we can write s 0 = π(X t ) in equation ( 20) which immediately leads to [START_REF] Deng | Solutions of Euler-Poisson equations for gaseous stars[END_REF] and completes the proof of our Proposition.] Although this property is geometrically quite obvious, let us provide a comprehensive proof for the sake of completeness. By definition [START_REF] Dembo | Large deviations techniques and applications[END_REF]

X * = X t * = Y * • σ * , σ * = Arg inf σ∈SN ||Y * • σ -s 0 ||.
Expanding the squares, we obtain

κ = ||s 0 -s 0 • σ|| 2 + 2t t * ((X * -s 0 , s 0 -s 0 • σ))
. By definition of X * = Y * • σ * and σ * , we get from ( 23)

||X * -s 0 || = ||Y * • σ * -s 0 || ≤ ||Y * • σ * • σ -1 -s 0 || = ||X * -s 0 • σ||. Thus ((X * , s 0 -s 0 • σ)) ≥ 0. So, we deduce κ ≥ ||s 0 -s 0 • σ|| 2 - 2t t * ((s 0 , s 0 -s 0 • σ)) = 2(1 - t t * )[||s 0 || 2 -((s 0 , s 0 • σ))
] ≥ 0 (using ||s 0 || = ||s 0 • σ|| and the Cauchy-Schwarz inequality), which is the desired inequality and completes the proof of Proposition 2.1.

Equation [START_REF] Deng | Solutions of Euler-Poisson equations for gaseous stars[END_REF] is clearly singular at time t = 0. We can lift this singularity with an exponential rescaling of time t = exp(θ), θ ∈ R, and finally obtain:

(24) dX θ dθ = X θ -π(X θ ) = X θ -∇Π(X θ ).
So we have obtained a first order (in time) dynamical system, of "gradient type".

Miscellaneous remarks. i) Since π is the gradient of a Lipschitz convex function (namely Π defined by ( 9)), equation ( 24) is uniquely solvable in the framework of "maximal monotone operator theory" [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. As already discussed in [START_REF] Brenier | A modified least action principle allowing mass concentrations for the early universe reconstruction problem[END_REF], this is a way to introduce a dissipative mechanism in the motion of particles, such as sticky collisions when particles stick to each other while conserving their momentum (but not their kinetic energy which decreases). At the level of the present paper, we do not want to enter such considerations and leave [START_REF] Feynman | The Feynman lectures in Physics[END_REF] just as a f ormal equation. ii) Notice that ( 24) is a gradient flow in the variable X, valued in the Hilbert space H, which is the counterpart of the gradient flow in the variable Law(X), in the so-called "Wasserstein space" of half the negative squared "Wasserstein distance", as discussed in [START_REF] Ambrosio | Savaré Gradient flows in metric spaces and the Wasserstein spaces of probability measures[END_REF]. For us, it is important to keep a formulation in terms of X ∈ H and not in terms of Law(X), because, in the sequel of our discussion, it is crucial to restore the "individuality" of the gravitating particles during their motion, in order to get a second order dynamical system for them. iii) Quite remarkably, as explained in [START_REF] Brenier | A modified least action principle allowing mass concentrations for the early universe reconstruction problem[END_REF], equation ( 22) is nothing but the Zeldovich model used in Cosmology [START_REF] Zeldovich | Gravitational instability: An approximate theory for large density perturbations[END_REF][START_REF] Shandarin | The large-scale structure of the universe: turbulence, intermittency, structures in a self-gravitating medium[END_REF][START_REF] Frisch | A reconstruction of the initial conditions of the Universe by optimal mass transportation[END_REF][START_REF] Brenier | Reconstruction of the early universe as a convex optimization problem[END_REF] as an approximation of semi-Newtonian gravitation in an Einstein-de Sitter space! iv) Let us finally provide a possible interpretation of equation ( 24), vaguely related to the so-called "anthropic principle": What is observed at out "present" time t * is just a random output of the independent Brownian trajectories of a large number N of indistinguishable particles initially located on the set A. As the noise vanishes, the motion of these particles, conditioned by what we can observe now, just looks driven by the deterministic law [START_REF] Feynman | The Feynman lectures in Physics[END_REF]. Of course, this is a highly speculative and questionable interpretation coming from a mathematician and not from an authorized physicist.

2.3. Second application of the large deviation principle. Equation ( 24) is a f irst order dynamical system, a so-called "gradient flow" since [START_REF] Feynman | The Feynman lectures in Physics[END_REF] can also be written dX θ = ∇Φ(X θ )dθ, where Φ is half of the squared distance function to S:

(25) Φ(Y ) = ||Y || 2 2 -Π(Y ) = inf s∈S ||Y -s|| 2 2 ,
where Π is defined by [START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF]. Describing a gravitational theory by a gradient flow does not sound reasonable. We would rather like to get a second order, conservative, dynamical system. Here again, the large deviation principle turns out to be useful to get such a second order system out of [START_REF] Feynman | The Feynman lectures in Physics[END_REF]. Fixing η > 0, we introduce the "noisy" version of (25) defined by: (26) dX η θ = ∇Φ(X η θ )dθ + η dB θ , where B θ (again) denotes a Brownian process in (R d ) N . Given two points Y 0 and Y 1 in (R d ) N , we expect from the LDP (or, more precisely, the Freidlin-Wentzell theorem [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF][START_REF] Dembo | Large deviations techniques and applications[END_REF]) that, as the level of noise η goes to zero, Strictly speaking, this is correct when ∇Φ is Lipschitz continuous, which is not true in our case (where ∇Φ is not even continuous), without further assumptions on the data. However, from the pure modeling viewpoint, it is very tempting to find the second order dynamical system linked to the least action principle [START_REF] Gangbo | Hamilton-Jacobi equations in the Wasserstein space[END_REF]. Since Φ, as defined by [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF], is half of a squared distance function, we have

-η log Prob[X η θ0 ≈Y 0 and X η θ1 ≈Y 1 ] ∼ η→0 A(θ 0 , Y 0 , θ 1 , Y 1 ) (27)
(28) 1 2 ||∇Φ(Y )|| 2 = Φ(Y ),
for every Y ∈ H \ N , where N is the set on which Φ is not differentiable, which is a negligible subset of H, both in the Lebesgue almost everywhere sense and in the Baire category sense. Thus, at least for each curve X that stays away from N for Lebesgue almost every θ ∈ [θ 0 , θ 1 ], we have From this equivalent LAP, we find as optimality equation (30) d 2 X θ dθ 2 = ∇Φ(X θ ) = X θπ(X θ ), which is just the second order version of the gradient flow equation [START_REF] Feynman | The Feynman lectures in Physics[END_REF]. More explicitly, we have obtained which is nothing but the discrete version [START_REF] Brenier | A geometric approximation to the Euler equations: The Vlasov-Monge-Ampère equation[END_REF] of the Monge-Ampère model of gravitation discussed in the first section. So we have achieved, at least at a formal level, the main goal of our paper which was, through a double application of the large deviation principle, the derivation of the (discrete) Monge-Ampère gravitational model from one of the simplest thinkable model of particles: N independent Brownian trajectories, with an intriguing interplay between their indistinguishability (for the first application of the LDP) and their individuality (for the second application of the LDP).

  with norm and inner product respectively denoted by || • || and ((•, •)). (Notice that, when A is a finite set of N points in R d , then H ∼ R N d is of finite dimension.) We crucially consider the subset S of all measure-preserving maps s of A: (8) S = {s ∈ H, A γ(s(a))λ(da) = A γ(a)λ(da), ∀γ ∈ C 0 (R d )}.

  A(θ 0 , Y 0 , θ 1 , Y 1 ) = inf{ X θ )|| 2 dθ, X θ0 = Y 0 , X θ1 = Y 1 }.In addition, as η goes to zero,X η θ ∼ X θ , ∀θ ∈ [θ 0 , θ 1 X θ )|| 2 dθ, X θ0 = Y 0 , X θ1 = Y 1 }.

|| 2 +|| 2 || 2

 222 ||∇Φ(X θ )|| 2 )dθ -+ Φ(X θ ))dθ -Φ(X θ1 ) + Φ(X θ0 ).Therefore, the action principle (27) is equivalent to(29) Ã(θ 0 , Y 0 , θ 1 , Y 1 ) = inf{ + Φ(X θ ))dθ, X θ0 = Y 0 , X θ1 = Y 1 }.

( 31 ) d 2

 312 X θ (a) dθ 2 = X θ (a) -A(σ θ (a)), whereσ θ = Arg inf σ∈SN N a=1 |X θ (a) -A(σ(a))| 2 ,

  of the closest point operator π, it is enough to show thatκ = ||X ts 0 • σ|| 2 -||X ts 0 || 2 ≥ 0for each permutation σ ∈ S N and each t ∈ [0, t * ]. By definition[START_REF] Cordero-Erausquin | Sur le transport de mesures périodiques[END_REF],κ = ||s 0 + t t * (X *s 0 )s 0 • σ|| 2 -||

	t t * (X * -s 0 )|| 2
	where we set, still according to (19),
	(23)

Acknowledgments. This work has been partly supported by the grant ISOTACE ANR-12-MONU-0013 (2012-2016). The author wants to thank the hospitality of the Radon Institute (RICAM, Linz) and the Hausdorff Institute (HIM, Bonn) during the special semester on New Trends in Calculus of Variations (Oct-Dec 2014) and the Junior Hausdorff trimester program on Optimal Transportation (Jan-Apr 2015). http://www.ricam.oeaw.ac.at/specsem/specsem2014 http://wt.iam.uni-bonn.de/conference-new-trends-in-optimal-transport/ He also thanks Christian Léonard for very stimulating discussions and comments bringing together large deviation and optimal transport theories.


Guillaume Perrin
email: guillaume.perrin@enpc.fr

Denis Duhamel
email: denis.duhamel@enpc.fr

Christian Soize
email: christian.soize@univ-paris-est.fr

Christine Funfschilling
email: christine.funfschilling@sncf.fr

Quantification of the influence of the track geometry variability on the train dynamics

Keywords: Railway track geometry, random fields, train dynamics, risk assessment

published or not. The documents may come

Quantification of the influence of the track geometry variability on the train dynamics

Introduction

To face the always more demanding challenges of the railway field, the expected benefits of simulation versus experiment are multiple, as it would allow cheaper, shorter, and more

T (s) B(s) N (s) x 1 (s) T (s) B(s) N (s) x 2 (s) T (s) B(s) N (s) 2x 3 (s) T (s) B(s) N (s) E+ 2x 4 (s)
Figure 1: Parametrization of the track irregularities, where E the theoretical track gauge, s is the curvilinear abscissa of the track and (O NT (s), T(s), N(s), B(s)) is the Frenet frame (for each rail, the mean position is represented in black, whereas the real position is in grey). practical procedures. However, if simulation has to be used in security, maintenance and comfort prospects, it has to be representative of the physical behavior of the track-vehicle system. The numerical models of the train and of the wheel-rail contact have indeed to be fully validated, and a procedure to accurately investigate the diversity of running conditions that the train can be confronted to during its life cycle has to be defined. Such running conditions refer, in particular, to a double scale description of the track geometry. Indeed, whereas the mean line position of the perfect track, which is chosen at the building of a new track, is characterized by three curvilinear quantities, which are the vertical curvature c V , the horizontal curvature c H , and the track superelevation c L , the description of four kinds of track irregularities has to be added to define the actual positions of the two rails, which are the lateral and vertical alignment irregularities x 1 and x 2 , the cant deficiencies x 3 and the gauge irregularities x 4 (see Figure 1). These four track irregularities, which are the main source of excitation of the train dynamics, are moreover in constant evolution, which is due to the train dynamics, to the modifications of the track substructure and to maintenance operations. Numerical methods to quantify the influence of this track irregularity variability on the train dynamics are therefore needed.

The general scheme for probabilistic analysis is usually divided in three steps [START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF][START_REF] Dutfoy | Reference Guide, Tutorial, Open TURNS[END_REF][START_REF] Funfschilling | Propagation of variability in railway dynamic simulations: application to virtual homologation[END_REF]. First, the mechanical model and the associated input parameters and output criteria have to be defined precisely. Three kinds of inputs are needed for a railway model to be defined: the vehicle model, the track model, and the wheel-rail contact model. Given these three inputs, the train response can be computed as the solution of a system of coupled equations that are strongly nonlinear. Once these equations have been solved, the spatial accelerations of each mass body, as well as the internal and external loads are available. These railway outputs can then be post-processed to define safety, comfort and maintenance criteria. Hence, thanks to the coupling of a series of breakthroughs in the modeling of complex mechanical systems [START_REF] Ding | Second order adjoint sensitivity analysis of multibody systems described by differential-algebraic equations[END_REF][START_REF] Hickey | Higher-order spectra for identification of nonlinear modal coupling[END_REF][START_REF] Iwnicki | Handbook of railway vehicle dynamics[END_REF][START_REF] Kerschen | Past, presentand future of nonlinear system identification in structural dynamics[END_REF][START_REF] Li | Comparison between karhunen-loève expansion and translation-based simulation of non-gaussian processes[END_REF][START_REF] Petsounis | Parametric time-domain methods for the identification of vibrating structures -a critical comparison and assessment[END_REF][START_REF] Shabana | On the computer formulations of the wheel/rail contact problem[END_REF][START_REF] Serban | Analytical derivatives for multibody system analysis[END_REF][START_REF] Sextro | Dynamical contact problems with friction[END_REF][START_REF] Wringen | Track to Car body vibration transfer for a bogie rail vehicle[END_REF] with an increase of the computational resources, it is now possible to compare the simulated dynamics of a real train on a measured track geometry and the measured dynamics of the same train on the same track [START_REF] Kraft | Parameter identification for a TGV model[END_REF]. Then, much attention has to be paid to the modeling of the input variability, as any error on the input will be propagated to the output. Only the uncertainty in the track irregularities is addressed in this work. Finally, the variability has to be propagated through the mechanical model. Noticing that the railway mechanical system is based on a very high number of variable input parameters, that the train response is very sensitive, very non linear, and very fuzzy with respect to these input parameters, it appears that the best method to propagate the track variability to the train dynamics is the Monte Carlo (MC) method [START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF]. If such a MC method is used, a method to generate independent realistic and representative track conditions is thus needed. In this prospect, methods based on the one-sided power spectral density (PSD) functions of the track geometry [START_REF] Fryba | Dynamics of railway bridges[END_REF][START_REF] Au | Impact study of cable-stayed railway bridges with random rail irregularities[END_REF][START_REF] Majka | Dynamic response of bridges to moving trains: A study on effects of random track irregularities and bridge skewness[END_REF] could be introduced to statistically characterize the random irregularities. From a time-frequency transformation based on a spectral representation [START_REF] Fries | A state-space approach to the synthesis of random vertical and crosslevel rail irregularities[END_REF][START_REF] Claus | Modeling and simulation of railway bogie structural vibrations[END_REF], track irregularities can indeed be generated from these PSD functions, which can be either estimated from measured track irregularities or from track safety standards such as the ones given by the Federal Railroad Administration [START_REF] Garg | Dynamics of railway vehicle systems[END_REF]. In that case, the railway track irregularities are however seen as four independent Gaussian random fields that are assumed to be ergodic and stationary in space. These track irregularities are nevertheless strongly dependent in practice, and due to the strong coupling of the train dynamics with the track degradation, they are actually neither stationary nor Gaussian. More recently, taking advantage of the recent developments in the modeling of non-Gaussian and non-stationary vector-valued random fields [START_REF] Das | Polynomial chaos representation of spatio-temporal random field from experimental measurements[END_REF][START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Ghanem | On the construction and analysis of stochastic models: Characterization and propagation of the errors associated with limited data[END_REF][START_REF] Le Maître | Spectral Methods for Uncertainty Quantification[END_REF][START_REF] Nouy | Generalized spectral decomposition method for stochastic nonlinear problems[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF][START_REF] Perrin | Karhunen-loève expansion revisited for vector-valued random fields: scaling, errors and optimal basis[END_REF][START_REF] Perrin | A posteriori error and optimal reduced basis for stochastic processes defined by a finite set of realizations[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-gausian tensor-valued random fields using partial and limited experimental data[END_REF], an original construction of the track irregularities distribution has been proposed in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF]. Such a modeling allows then the generation of running conditions that are realistic from frequency and statistical points of view, and representative of the measured track quality.

Hence, the idea of this work is to show to what extent such a complex modeling of the track geometry can be coupled with a rigid-body modeling of a complete train to analyze the influence of the track irregularity variability on the train dynamics. To this end, Section 2 introduces first the railway stochastic modeling and its validation from experimental data. Section 3 underlines then how such a stochastic modeling of the railway system can be used to perform robust comparisons between trains with different mechanical properties and to quantify the consequences of an increase of the train speed on the train stability and aggressiveness.

For confidentiality reasons, only qualitative analysis will be presented in this work as very few numerical values will be given.

Railway stochastic modeling

This section aims at presenting the numerical frame, on which the quantification of the influence of the track variability on the train dynamics is based. First, a description of the railway deterministic problem is introduced. Then, the railway stochastic modeling is described. The relevance of both deterministic and stochastic problems is moreover evaluated from comparisons to experimental data.

Railway deterministic modeling 2.1.1. Definition of the deterministic modeling

As presented in Section 1, a deterministic railway simulation is based on the introduction of three kinds of inputs.

• The vehicle model V. Multibody simulations are usually employed to model the train dynamics. Carbodies, bogies and wheelsets are therefore modeled by rigid bodies linked with connections represented by rheologic models (dampers, springs, ...). For 1 ≤ i ≤ N DoF , and t in [0, T ], we denote by u i (t) the position at time t of the coordinate associated with each degree of freedom of the rigid body modeling of the train, and by ui (t) = du i dt (t) its time derivative. For instance, for a classical one-carriage TGV, which is made of 10 coaches, 13 bogies and 52 wheelsets that are linked by a series of suspensions and bumpstops, N DoF is about two hundreds (see Figure 2 for a simplified representation of the TGV).

• The track geometry, T . As described in Section 1, the track geometry refers to a double scale description. On the first hand, the track design gathers the horizontal curvature c H , the vertical curvature c V and the cross level c L , and corresponds to the description of a perfect track without irregularities. On the other hand, four track irregularities, x 1 , x 2 , x 3 and x 4 , have been added to characterize the deviation of the real railway track toward this perfect track. It is moreover recalled that whereas the track design is chosen constant, such irregularities are in constant evolution.

• The contact model, C, allows the computation of the contact forces between the rails and the wheels. In the railway community, these contact forces are almost always computed from the wheel profile and the rail profile thanks to the Hertz and Kalker theories [START_REF] Hertz | Ueber die beruhrung fester elasticher körper[END_REF][START_REF] Kalker | Three-Dimensional elastic bodies in rolling contact[END_REF].

Introducing the vector of the generalized coordinates, at time t,

U (t) = (u 1 (t), u 2 (t), . . . , u1 (t), u2 (t), . . .) , (1)
the train dynamics can therefore be determined by solving the Euler-Lagrange equation, which is written as:

d dt ∂E c ∂ ui - ∂E c ∂u i = L i (U , T , C), 1 ≤ i ≤ N DoF , (2)
with E c the total kinetic energy of the train, and L i (U , T , C) the general load that is applied to the degree of freedom i, which depends on the track geometry T , on the wheel-rail contact C and on the generalized coordinated U . Eq. ( 2) can be rewritten in a matricial form as:

[A(U )] U = F (U , T , C), (3)
with [A] and F two strongly nonlinear operators. This system is usually solved with an explicit time scheme. In the following, the commercial black-box software, Vampire (see [START_REF] Evans | Rail vehicle dynamic simulation using vampire[END_REF][START_REF] Evans | Vampire User guide[END_REF] for further details about this software), is used. The chosen time step of this explicit scheme was identified according to a convergence analysis and is generally taken equal to 10 -4 second. The generalized coordinates vector U is then post-treated to define the final comfort and safety criteria associated with the railway system. In this work, five representative outputs are chosen to characterize the train dynamics, which can be classified in two categories.

1. First, the maximal values of the vertical and lateral accelerations in the train coaches, zmax and ÿmax , are controlled to guarantee the comfort of the passengers. 2. Secondly, the safety and maintenance criteria of the track-vehicle system are based on the analysis of the wheel-rail contact forces. In this prospect, the following three criteria are generally introduced to characterize the vehicle dynamics on a given track geometry of total length S tot :

• a shifting criterion:

(Y ℓ + Y r ) max = max wheelset w max 0≤s≤S tot {Y w ℓ (s) + Y w r (s)} , (4)
• a derailment criterion:

(Y /Q) max = max wheel q max 0≤s≤S tot {Y q (s)/Q q (s)} , (5)
• a wear criterion:

(T γ) = wheel q S tot 0 T q (s)γ q (s)ds , (6)
where:

• 0 ≤ s ≤ S tot is the curvilinear abscissa of the track;

• Y w ℓ and Y w r are the left and right lateral forces of the same wheelset w, such that the higher (Y ℓ + Y r ) max , the more chance for a shifting of the track;

• Y q and Q q are the lateral and vertical components of the wheel-rail contact force at wheel q, such that the higher (Y /Q) max is, the more on the flange a wheel of the train can be;

• T q and γ q are respectively the creep force and the slip at wheel q, such that the higher (T γ), the higher the contact wear for one run of the complete train.

Finally, given a model of the wheel-rail contact C, the deterministic railway problem corresponding to the dynamics of a vehicle V on a track geometry T can be expressed as:

(V, T , C) → c = g (V, T , C) , c = (z max , ÿmax , (Y ℓ + Y r ) max , (Y /Q) max , (T γ)) , (7)
where it is reminded that g is a complex and nonlinear operator. These nonlinearities are mostly due to the train suspensions (especially the airsprings between the bogies and the coaches), to a series of bumpstops in the train description and to the wheel-rail contact forces. Moreover, due to the train dynamics, to the track irregularities and to the specific wheel and rail profiles, the contact positions between each wheel of the train and the rails keep changing. The wornest the track geometry is, the more discontinuous these changes are likely to be. For instance, the diversity of these contact positions and contact forces can be seen in Figure 4, which is based on the run of a train on a measured track geometry around a curve, whose design is shown in Figure 3.

Domain of validity for the deterministic problem

As a first comment on the validity of the deterministic railway modeling, it is important to point out that all European railway reference standards and reference maintenance guides only consider the low-frequency content, f ≤ f c , of the train dynamic quantities of interest (either simulated or measured). As presented in the former Section, the software Vampire is used to solve the railway deterministic problem. The train being constituted of rigid bodies, the simulated high-frequency response of the train cannot be physical. As an illustration, Figure 5 compares the measured and simulated frequency properties of a bogie of a TGV. As shown in Figure 6, although the transverse and vertical accelerations of the bogie are low-pass filtered at the reference cut-frequency f = f c , it can be seen that the low-frequency response is well reproduced both in the time and frequency domains by the deterministic model. As a consequence, in agreement with the work achieved in [START_REF] Kraft | Parameter identification for a TGV model[END_REF], it is assumed that the proposed railway deterministic model is valid on the frequency band 0 ≤ f ≤ f c . In the following, each output of the train dynamics (whether measured or simulated) will thus be low-pass filtered at frequency f c before being analyzed. ), the vertical curvature c V (×km -1 ) and of the cross level c L (×m -1 ) with respect to the track curvilinear abscissa s.

Railway stochastic problem 2.2.1. Definition of the stochastic problem

As presented in Section 1, only the track geometry variability is analyzed in this work. The wheel and rail profiles of high speed trains and lines being checked and maintained very regularly, only perfect wheel and rail profiles will be considered in the following, such that the contact properties, C, are chosen to be constant. It is moreover supposed that the track irregularities can be separated from the track design. Hence, in the following, the track design is also supposed to be constant, while the four track irregularities can vary. As a consequence, vector c, which is defined by Eq. ( 7), becomes a random vector that is denoted by C = (C 1 , C 2 , C 3 , C 4 , C 5 ). It is reminded that by definition of vector c, C 1 and C 2 refer to the vertical and lateral maximal accelerations in the train coaches, C 3 is the maximal value of the sum of the transverse loads of the wheelsets, C 4 is the maximal value of the Y /Q ratio, and C 5 is the cumulated wear along the track. At last, given a fixed description of the track design, (c H , c V , c L ), and a normalized model of train, V, for which mechanical parameters are also fixed and have been accurately identified, the railway stochastic problem can be witten:

X → C = G (X | c H , c V , c L , V, C) , (8)
where X = (X 1 , X 2 , X 3 , X 4 ) is a non-Gaussian and non-stationary four-dimensional random field, which gathers the evolutions with respect to the curvilinear abscissa of the considered track design of the four track irregularities that were introduced in Section 1. In this work, the distribution of X is moreover supposed to be identified according to the developments achieved in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF], from the measurement of the track geometry of a whole railway network. Hence, coupled with any track design, we admit that such a construction for X allows the generation of independent running conditions that are realistic and representative of the available measured track irregularities.

Validation of the stochastic problem

Two validations for the track geometry stochastic modeling, based on the train dynamics, are proposed in this section. In a first step, it is shown that the track generator coupled with the Vampire software allows the simulation of train accelerations that are similar to accelerations that have been recorded on a real high speed train on a real track. In a second step, another type of validation is proposed to evaluate the relevance of the track stochastic modeling for the analysis of the wheel-rail forces, as these forces can hardly be measured.

Relevance of the track stochastic modeling for the analysis of the train accelerations.

Since 2007, the measurement train IRIS-320 keeps monitoring the track geometry of the French high speed lines. The rigid body modeling of this train has therefore been achieved, simulations have then been performed at constant speed S on ν = 500 variable running conditions of total length S tot . For each simulation, the track design is therefore chosen to be constant and to correspond once again to the one given by Figure 3, whereas the track irregularities, which correspond to a particular realization, X(Θ n ), of X, are supposed to be different from one simulation to another. Hence, for all s in [0, S tot ], we respectively define C sim z (Θ n , s) and C sim y (Θ n , s) as the vertical and lateral maximal values, at position s, of the accelerations in all the coaches of the train that is excited by the track irregularity X(Θ n ). Given these two sets of train responses,

C sim z (Θ n ), 1 ≤ n ≤ ν and C sim y (Θ n ), 1 ≤ n ≤ ν , let {D z i (s), s ∈ [0, S tot ], 1 ≤ i ≤ 10}
and {D y i (s), s ∈ [0, S tot ], 1 ≤ i ≤ 10} be the decile functions, such that at each position s, i/10 of the values of C sim z (Θ n , s) and C sim y (Θ n , s) are in D z i (s) and D y i (s) respectively. These decile functions, whose representations are shown in Figure 7, allow us to evaluate the influence of the track irregularity variability on such maxima.

The IRIS-320 train is moreover equipped with accelerometers that record the vertical and transverse accelerations at three coaches,

ÿ(1) C , ÿ(2) C , ÿ(3) C , z(1) C , z(2) C , z (3)

C

.

In order to evaluate the relevance of the former results for the maximal accelerations in the train coaches, we define C exp z and C exp y , such that for any value of the curvilinear abscissa of the track, s, we get: , are then extracted from the experimental database collected by measurement train IRIS 320. These measurements were chosen as their dynamic characteristics were the most comparable to the simulated one, in terms of cross level, horizontal and vertical curvatures, speed of the train and length of the curve. If the chosen simulated dynamic characteristics were not similar to the extracted dynamic characteristics on the complete domain [0, S tot ], non-valid domains were added to these figures. The evolutions of these measured accelerations are compared to the simulated ones in Figure 8. In the light of these results, the track generator coupled with the Vampire software seems to be able to simulate realistic and representative runs of the IRIS-320 train to analyze the link between the two first quantities of interest of the stochastic modeling, C 1 and C 2 , and the track geometry variability.

C exp z (s) = max i∈{1,2,3} z(i) C (s) , (9)
C exp y (s) = max i∈{1,2,3} ÿ(i) C (s) . ( 10
) Vertical acceleration 0 0 S tot D z 10 (s) D z 9 (s) D z 8 (s) D z 7 (s) D z 6 (s) D z 5 (s) D z 4 (s) D z 3 (s) D z 2 (s) D z 1 (s) C sim z (Θ 1 ) Abscissa s Transverse acceleration 0 0 S tot D y 10 (s) D y 9 (s) D y 8 (s) D y 7 (s) D y 6 (s) D y 5 (s) D y 4 (s) D y 3 (s) D y 2 (s) D y 1 (s) C sim y (Θ 1 )

Dynamic validation of the track generator for the analysis of the wheel-rail contact forces.

No on-track measurements of the contact forces between the train and the track at high speed being available, an other approach is proposed to evaluate the relevance of the track generator to simulate realistic and representative values for C 3 , C 4 and C 5 . To this end, the particular curve of total length S tot shown in Figure 3, is once again considered. From the available measurements of the track geometry, ν exp = 400 different track conditions of total length S tot , {X exp (θ 1 ), . . . , X exp (θ ν exp )}, are gathered. The same normalized high-speed train V, for which mechanical parameters are supposed to be accurately identified, is thus made run first on the ν exp measured track conditions, and then on the ν = 500 formerly generated track conditions, {X(Θ 1 ), . . . , X(Θ ν )}, at the same speed S. Eight quantities of interest that are representative of the train dynamics are then compared:

• the left and right transverse contact forces at the first wheelset of the first bogie of the motor car,

Q 1 = Y ℓ M C and Q 2 = Y r M C ;
• the left and right transverse contact forces at the second wheelset of the second bogie of the second passenger car,

Q 3 = Y ℓ P C and Q 4 = Y r P C ;
• the left and right Y /Q ratio at the first wheelset of the first bogie of the motor car,

Q 5 = (Y /Q) ℓ M C and Q 6 = (Y /Q) r M C ;
• the left and right Y /Q ratio at the second wheelset of the second bogie of the second passenger car,

Q 7 = (Y /Q) ℓ P C and Q 8 = (Y /Q) r P C .
For 1 ≤ i ≤ 8, we are now interested in the mean power spectral densities of Q i and the mean numbers of upcrossings of the level u by Q i over the length S tot , which are respectively denoted by P SD mes (Q i ) and N mes up (Q i , u, S tot ) when these quantities are computed from the measured track geometries and P SD gen (Q i ) and N gen up (Q i , u, S tot ) when these quantities are computed from the generated track geometries. The comparisons between these quantities are represented in Figure 9, where a very good fit can be seen (for confidentiality reasons, replacemen Abscissa it is reminded that no numerical values for the abscissa of these graphs are given). The stochastic modeling of the track geometry, based on the developments carried out in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF], is thus relevant from the train response point of view. In the following, it is therefore supposed that the track stochastic modeling, coupled with the software Vampire is also relevant to investigate the relation between the track variability and the three quantities of interest C 3 , C 4 and C 5 .

s Vertical 0 0 S tot D z 10 (s) D z 1 (s) C exp,(1) z NVD Abscissa s Transverse 0 0 S tot D y 10 (s) D y 1 (s) C exp,(1) y NVD Abscissa s Vertical 0 0 S tot D z 10 (s) D z 1 (s) C exp,(2) z Abscissa s Transverse 0 0 S tot D y 10 (s) D y 1 (s) C exp,(2) y Abscissa s Vertical 0 0 S tot D z 10 (s) D z 1 (s) C exp,(3) z NVD Abscissa s Transverse 0 0 S tot D y 10 (s) D y 1 (s) C exp,(3) y NVD Abscissa s Vertical 0 0 S tot D z 10 (s) D z 1 (s) C exp,(4) z Abscissa s Transverse 0 0 S tot D y 10 (s) D y 1 (s) C exp,(4) y Abscissa s Vertical 0 0 S tot D z 10 (s) D z 1 (s) C exp,(5) z Abscissa s Transverse 0 0 S tot D y 10 (s) D y 1 (s) C exp,(5) y

Influence of the track irregularities variability on the train dynamics

As explained in Section 1, a better understanding of the specific link between the track irregularities and the train response is needed to optimize the maintenance, and to better anticipate the consequences of modifications of the running conditions. To this end, we denote by P C (dx) = p C (x)dx the multidimensional distribution of random vector C, where p C is the associated density. This distribution is strongly related to the distribution of the track irregularity random field, P X (see Eq. ( 8)). Assuming that the latter distribution has been accurately identified, the track variability has now to be propagated through the railway model to characterize P C . As the statistical dimension of X is very high and as the relation between P C and P X is very complex and strongly nonlinear, the Monte Carlo method is chosen to evaluate P C . From ν independent realizations of X, {X(Θ 1 ), . . . , X(Θ ν )}, which correspond to ν independent realistic and representative running conditions that a train can be confronted to during its lifecycle, ν independent realizations of C, {C(Θ 1 ), . . . , C(Θ ν )}, can be deduced as:

C(Θ n ) = G (X(Θ n ) | c H , c V , c L , V, C) , 1 ≤ n ≤ ν. (11)
The statistical properties of C are finally deduced from the analysis of this ν-dimensional set of independent realizations of C. Three applications of this stochastic modeling are now presented. These are based on the track design of total length S tot shown in Figure 3, and on ν = 4, 000 independent running conditions. First, the compared influence of the track design and of the track irregularities is illustrated. Then, it is shown to what extent such a method can be used to quantify the influence of an increase of the train speed on C, and to compare the safety and the aggressiveness of three different high speed trains.

Influence of the track design

The idea of this section is to quantify the importance of the track irregularities and of the track design on vector C. In this prospect, the response of a normalized high train V to the former ν track conditions of total length S tot is analyzed. Four categories are considered: the alignment (A), the curve entrance (CE), the established curve (C) and the curve exit (SC). The response of the train is therefore sorted with respect to these four curve categories, such that, for 1 ≤ i ≤ 5, four values of the railway quantities of interest

C A i (Θ n ), C EC i (Θ n ), C C i (Θ n ) and C SC i (Θ n ) can be computed.
Based on these sets of ν independent realizations, the PDFs of the components of C are estimated from a kernel smoothing method [START_REF] Dutfoy | Reference Guide, Tutorial, Open TURNS[END_REF], and are represented in Figure 10. From these graphs, it can be seen that the influence of the track design on the wear criterion, C 5 is very high. The other dynamic quantities, C 1 , C 2 , C 3 and C 4 seem however to be much more dependent on the track irregularities than on the track design.

Influence of an increase of the speed on the quantities of interest

The second application deals with the influence of the train speed on the PDFs of the five considered criteria. Only the established curve configuration is shown. Railway simulations are therefore performed on the same ν realistic and representative track geometries, at the four speeds S1 = S, S2 = 1.1S, S3 = 1.2S and S4 = 1.3S. Two other sets of simulations have then been carried out for a different value of the track superelevation, c * L , at speeds S3 = 1.2S and S4 = 1.3S in order to quantify the importance of this track design parameter with respect to the five studied criteria. In other words, whereas c L is chosen to compensate the train inertial acceleration in curve at speed S1, c * L allows the compensation of the train inertial acceleration in curve at speed S3. For each speed, the PDFs of each component of C are once again estimated using a kernel smoothing method based on the ν = 4, 000 independent railway simulations. These PDFs are represented in Figure 11. In this figure, the nonlinearity of the system can be noticed, as the consequences of an increase of the speed of 10% to 30% are much higher than 30% for each criterion. In particular, an increase of 30% of the speed of the train can yield an increase of more than 500% of the contact wear if the track superelevation is not adjusted. In addition, these figures emphasize the importance of the adjustment of the track superelevation to the speed, in terms of minization of the wear, of the shifting and of the risk of derailment.

Comparison of three high speed trains

In this section, it is supposed that three different models of three competitive high speed trains, V 1 , V 2 and V 2 are available. The mechanical parameters of these trains are very different and were carefully identified from experimental measurements. These three trains are thus made run on the same ν track geometries at the same speed S. The PDFs of each criterion C i associated with each train are then shown in Figure 12. Hence, the stochastic modeling allows us to compare the dynamic response of these three trains when excited by a representative set of realistic track conditions. In particular, criteria C 3 and C 5 could be interesting indicators to compare the aggressiveness of each train.

Conclusions and prospects

A method to propagate the track geometry variability through railway mechanical simulations is nowadays of great interest to face always more challenging railway issues. In this prospect, this paper has presented a method based on the coupling of a track geometry stochastic modeling with a rigid-body modeling of a train. Based on experimental data, this stochastic modeling of the track geometry allows the generation of running conditions that are realistic and representative of the quality of the measured railway network. In this work, it has therefore been shown to what extent such a complex modeling, which takes into account the dependencies between the different track irregularities, can be used to analyze the complex link between the track geometry variability and the train dynamics and stability. In particular, three applications have been described in this paper. First, the compared influence of the track design and of the track irregularities has been analyzed. Then, the impact of an increase of the speed on the train stability has been presented. At last, it has been shown to what extent such an approach could be used to compare competitive high speed trains with respect to their response on a set of representative track conditions.

Figure 2 :

 2 Figure 2: Simplified description of a multibody model of a TGV.

Figure 3 :

 3 Figure 3: Evolution of the horizontal curvature c H (×km -1), the vertical curvature c V (×km -1 ) and of the cross level c L (×m -1 ) with respect to the track curvilinear abscissa s.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Evolution of the contact conditions with respect to the train dynamics and to the horizontal curvature (CP=contact point, CF=contact force).

Figure 7 :

 7 Figure 7: Influence of the track variability on the maximal values of the accelerations in the train coaches.

Figure 8 :

 8 Figure 8: Comparison between simulated and measured maximal values of the vertical and transverse accelerations in the train coaches (NVD = non-valid domain). 13

Figure 10 : 1 PDFsFigure 11 :

 10111 Figure 10: Influence of the track design on the marginal PDFs of vector C.

Level u

Upcrossings


Keywords:

Timing of rare-elements (Li-Be-Ta-Sn-Nb) magmatism in the European Variscan belt

Jérémie Melleton, Eric Gloaguen

Timing of rare-elements (Li-Be-Ta-Sn-Nb) magmatism in the European Variscan belt J. MELLETON, E. GLOAGUEN BRGM,Direction Géoressources,ISTO,UMR 7327,B.P. 36009, High-phosphorus peraluminous rare-elements granites and rare-elements LCT (Lithium, Caesium, Tantalum) pegmatites are the most important sources of raw materials for some critical metals like tantalum (1,2)represent important economic storehouses for industrial minerals like feldspar, quartz, mica or kaolin. They principally emplace in orogenic settings (3). Afast overview of three mainEuropean Variscan districts, i.e. the Moldanubian domain of the Bohemian massif, the French Massif Central (FMC) and the NW Iberia provides a basis for questioning the origin of rare-elements magmatism and the actual classification of rare-elements pegmatites, in particular the LCT pegmatites.

Granitic pegmatites are widespread in most of the Bohemian Massifbut LCT pegmatites are most common in the Moldanubian domain. In this area, their emplacements seem mainly controlled by migmatitic domes and shear zones and correspond to two events(4). The older at ~ 333 ± 3 Ma just follow HT-MP event of the end of the Moravo-Moldanubian phase and the younger at ~ 325 ± 4 Ma is contemporaneous with beginning of the Bavarian phase (U-Pb ages on colombite and tantalite).

In the FMC, most of the actually known rare-elements magmatic bodies form a province in the North Limousin area, which represents the northwestern part of the FMC.U-Pb dating of columbite-group minerals from Beauvoir, Montebras and Chèdeville rare-elements magmatic bodies leads to emplacement ages at 317 ± 6 Ma, 314 ± 4 Ma and 309 ± 5 Ma respectively. The contemporaneous Marche fault system (5), which crosscuts in a general E-W trend all the northern part of the Limousin, seems to be a key-structure for the rare-elements magmatism of the area.

Although rare-elements pegmatites are known in the different Variscan massifs of Iberia, the northwest part of the Iberian Variscan belt contains numerous fields that represent the first economic targets in Europe, particularly in the Central Iberian Zone (CIZ) and in Galicia-Trás-os-Montes Zone (GTOMZ)Three events of rare-elements magmatism have been recognized in Northwest Iberia (U-Pb on colombite and tantalite; 6): (i)emplacement of the Argemela rare-element granite, in the CIZ, with an age of 326 ± 3 Ma; (ii) intrusion of rare-element pegmatites from the GTOMZ at an average age of 310 ± 5 Ma; (iii)emplacements of rare-element pegmatites in the CIZ and in the southern GTOMZ at about 301 ± 3 Ma. Moreover, the observed southward propagation of ages of emplacement matches the propagation of deformation, metamorphism and magmatism in the two different geotectonic zones.

The spatial and temporal distributions of rare-element pegmatites in the Variscan beltsuggest that rare-elements magmatism could be related to local specific conditions like particular sources, tectonic and thermal regimes.

1) Černý&Ercit(2005). Can.Mineral., 4, 2005Mineral., 4, -2026. . 2) Linnen&Cuney(2005). In: Linnen RL, Samson IM (Eds.) Rare-Element Geochemistry and Mineral Deposits. Geological Association Canada Short Course Notes, 17, 45-68. 3) Martin & De Vito (2005). Can. Mineral., 2027-2048.4) Melleton et al. (2012). Can. Mineral., 50, 1751-1773. 5) Gébelin et al. (2007). Tectonics, TC2008, doi:10.1029/2005TC001822. 6) Melleton et al. Submitted to Mineral.Pet.


Abdelhamid Makiou
email: makiou@telecom-paristech.fryouceibegriche@ieee.orgahmed@telecom-paristech.fr

Youcef Begriche

Ahmed Serhrouchni

Improving Web Application Firewalls to Detect Advanced SQL Injection Attacks

Keywords:

Injections flaws which include SQL injection are the most prevalent security threats affecting Web applications[1].

To mitigate these attacks, Web Application Firewalls (WAFs) apply security rules in order to both inspect HTTP data streams and detect malicious HTTP transactions. Nevertheless, attackers can bypass WAF's rules by using sophisticated SQL injection techniques. In this paper, we introduce a novel approach to dissect the HTTP traffic and inspect complex SQL injection attacks. Our model is a hybrid Injection Prevention System (HIPS) which uses both a machine learning classifier and a pattern matching inspection engine based on reduced sets of security rules. Our Web Application Firewall architecture aims to optimize detection performances by using a prediction module that excludes legitimate requests from the inspection process.

I. INTRODUCTION

Structured Query Language (SQL) injection is one of the most devastating vulnerabilities that impacts DataBase Man agement Systems (DBMS), as it can lead to the exposure of all the sensitive information stored in an application's database [START_REF] Ansari | SQL Injection in Oracle: An exploration of vulnerabilities International Journal on Computer Science and Engineering[END_REF]. In order to confront SQL injection attacks, various methodologies and techniques have been used. On one hand, Web application developers adopted safe coding and applied input validation functions. They developed filters to protect the application's entries from SQL code injections. These filters block inputs that contain SQL keywords or special characters commonly used in malicious SQL code injection. On the other hand, Web Application Firewalls protect application's database by inspecting HTTP traffic and applying a set of security rules. The language used for expressing security rules can explicitly describe a signature of an SQL injection attack, or implicitly describe the way of detecting these attacks. It can also express an anomaly score value which increases every time a malicious pattern appears in an HTTP request. If the anomaly value reaches a predefined threshold, the request will be rejected. In spite of the robustness of the above methods, attackers can bypass them by substituting malicious pattern characters or varying its format. However, a basic solution is to write a specific rule for each type of evasion technique, but it requires a high mastery of both HTTP protocol and regular expressions programming. Furthermore, pattern matching algorithms, used by security rules in order to inspect complex patterns, decrease the overall performances of the detection engine. Since pattern matching algorithms require a lot of resources, some WAFs [START_REF] Ristic | ModSecurity Handbook: The Complete Guide to the Popular Open Source Web Application Firewall[END_REF] [9] are configured to only inspect POST request. In this paper, we propose a hybrid approach to detect SQL injection attacks and their evasion techniques. Our proposal enhances both the inspection process of HTTP streams and security rules management.

Our contributions can be summarized as follows: 1) A brief survey on filters evasion techniques including code obfuscations and headers-level SQL injections.

2) A novel method to dissect and parse the HTTP protocol, that enhances security rules management, which makes the formalism of writing security rules less complex.

3) A Hybrid Injection Prevention System (HIPS) architecture, by introducing a supervised machine learning module. The classifier (machine learning module) predicts SQL injections and forwards suspicious headers to the rules-based detection engine in order to be deeply analyzed . 4) An evaluation of the effectiveness of the proposed method based on false negatives impact instead of false positives which are commonly used in Anti-spams cost evaluation.

II. RELATED WORK

In [START_REF] Tajpour | Comparison of SQL Injection Detection and Prevention Techniques, 2nd International Conforence on Education Technology and Computer (ICETC)[END_REF][4], Kruegel and Vigna propose an anomaly-based in trusion detection system for web applications. It characterizes HTTP requests using a number of statistical characteristics de rived from parameter's length, character distribution, structure, presence and order. This method focuses only on the incoming query parameters whereas it ignores the corresponding HTTP response. These results are either causing unnecessary false positives or missing certain attacks. AMNESIA [START_REF] Halfond | AMNESIA: Analysis and Monitoring for NEutralizing SQL-Injection Attacks[END_REF] is an SQL injection detection and prevention system which combines static analysis and run-time monitoring. It uses a model-based approach to detect illegal queries. Nonetheless, it requires web application's source code reviewing. In SQLrand [3][6] instead of normal SQL keywords developers create queries using randomized instructions. In this approach a proxy filter intercepts queries to the database and de-randomizes the keywords. By using the randomized instruction set, attacker's injected code could not have been constructed. As it uses a secret key to modify instructions, security of the approach is dependent on attacker ability to seize the key. It requires the integration of a proxy for the database in the system as the same as developer training. ModSecurity WAF, proposed by Ivan Ristic [START_REF] Ristic | ModSecurity Handbook: The Complete Guide to the Popular Open Source Web Application Firewall[END_REF], is an open source solution based on signature attack detection. ModSecurity is widely used and has medium performances. Though, this system is strongly related to some types of web servers and it only analyses POST queries in order to avoid performance deterioration. In addition, the rules formalism is very complex which requires a high expertise in HTTP protocol and in regular expressions coding. IronBee [START_REF]The IronBee Project[END_REF], a new project similar to ModSecurity, aims to improve detection performance and facilitate the expression of security rules by introducing the LUA scripting language. An other recent open source project NAXSI [START_REF]Naxsi project (Nginx Anti Xss Sqllnjection)[END_REF] uses a heuristic ap proach for the detection of XSS and SQL injection attacks. Its performances are acceptable but requires a learning process to define white-lists. SQLi and XSS rules are static and use a simple cumulative scoring system based on the appearance of some special strings.

III. SQLI AT TACKS AND EVASION TECHNIQUES

Halfond, Viegas, and Orso researches [START_REF] Halfond | AMNESIA: Analysis and Monitoring for NEutralizing SQL-Injection Attacks[END_REF] proposed a tax onomy of SQL injection attacks. Depending on the goal, attackers can append a syntactically correct SQL code to the original query, or forge their own malicious SQL commands and introduce them to the DBMS via a vulnerable web appli cation inputs. They will use one of these classes: Tautology, Incorrect Queries, UNION Queries, Piggy-backed Queries, Stored Procedure, Inference (Blind Injection, Timing Attacks).

In the following, we will present two techniques attackers can use in order to bypass the WAF rules.

A. SQL injections by code obfuscations

Despite the deployment of SQLi attacks detection systems, attackers can manage to overcome these systems by trying to disguise the appearance of their requests. This is due to both the flexibility of the SQL language and the complexity of writing security rules wich express all scenarios of a masked attack. The combined use of ASCII encoding, Hexadecimal, and trans-coding functions, which are available by the CGIs, makes the pattern matching inspection process more complex. Indeed, the pattern matching algorithms are able to detect variations in the forms but they can't detect trans-coding of string characters. Thus, attacks detection systems start with the transformation (sanitization) of the query before submitting it to the pattern matching engine. Despite these efforts, it is always possible to bypass filters posed by these systems, because both the flexibility of DBMS languages and the semantic content of the query are not known by the detection systems. More details are shown in the following example : -Uppercase/lowercase blend: Some filters often ignore the situation of mixed case, such as the input AND or and which can be replaced by AnD -Tautology: (where name = '='), (where false = " "), (OR true like true) -SQL keywords: SELECT becomes sel/**/ect , AND becomes && or AandND -The use of trans-coding functions : (char(44) = ') , (exec(char(Ox73687574646j776e» is a system shutdown.

B. Headers-level SQL injection

The injection of SQL codes is not only restricted to GET and POST methods through URLs or request's body. Unfortunately, other more injection techniques exploit the open paradigm of the HTTP protocol in order to avoid web applications' security mechanisms. Attackers exploit HTTP protocol headers because most defense systems analyze only few or probably do not analyze HTTP headers due to the huge resources requirement.

1) User-agent Header: This header is used by web applications to store information about the client application (web browser, audit application, etc). Some e-commerce applications store this information in databases for clients profiling. An attacker can exploit this header to insert SQL code that overcome security checks. He can insert a tautology in order to make the check condition always valid. Example : user-agent : Firefox' OR 1=1.

2) Referer Header: Referer is another HTTP header which can be vulnerable to SQL injection once the application is storing it in database without sanitizing it. An attacker may inject arbitrary SQL conunands to the database by exploiting the Referer header. With the same way as in the User-Agent, the attacker can insert a tautology to make the control over the completely useless Referer or insert SQL malicious code after separator: Referer: 1','0'); SQL 3) X-Forwarded-For Header : This header is used to get the original IP address of a client who is connected through a proxy or other load-balancing device. An attacker can also make the control obsolete by inserting a tautology. Example : X-Forwarded-For :ip-address or 1=1 4) Cookie Header : Web applications use the cookie header in order to ease the authentication process of their clients. It often happens that the application does not validate the cookie passed by the client and inserts it directly in an SQL statement that selects the user session. An attacker can then inject SQL code in this header and completely avert the application authentication mechanism. He can use tools such as cookie manager to forge a valid cookie. Cookie variables sometimes are not properly sanitized before being used in SQL query. The cookie contains base64 encoded form identifier, a field that is unknown and a password. If we use as a cookie 12345 UNION SELECT mypass :: mypass base64 encoded, the SQL query becomes: SELECT user-password FROM nk-users WHERE user id=12345 UNION SELECT mypass IV. PROPOSED SOLUTION

A. HTTP Protocol Dissection

In most security solutions, traffic dissection process is the first operation before applying any security control. In HAKA project [10], HTTP stream is divided in two tables, HTTP request and HTTP-response, each table is associated with one hook and each hook contains all security rules declared for either requests or responses. In this section we introduce an other manner to dissect HTTP streams.

a) Typical HTTP request: HTTP protocol is expressed in a human-readable ASCII text. Headers use text to describe a request from a client (browser) or a response from the server. An HTTP request begins usually with a GET or POST method, followed by the URL and the protocol version. The following headers provide various information about the client, connection, content, etc. These headers are separated by \r\ n to distinguish each header. which are separated by \r\ n characters. However, before making the dissection, it has to get information about security rules. Indeed, users are obliged to declare security rules for the body and for each header. With the knowledge of headers involved in the inspection process, the dissector will only extract and parse these headers,

B. The Hybrid Injection Prevention System (HIPS)

• 1 Request Dissector URL: The information that the URL must be inspected is known by the dissector, since the user has declared a security rule on the hook : HTTP URL. The dissector extracts URL string and passes it to the classifier.

• 2 If legitimate: The result of the classification is negative, which means that the URL string does not include any potential SQL injection code,

• 3 If SQLi: The classifier decides to forward the URL string to the security rules detection engine, because the URL contains probably an SQL injection code, In the result section, we will show the minimum threshold for which a legitimate content is considered suspicious.

• 4 If SQLi: The detection engine loads all the security rules hooked to the HTTP URL hook

• 5 If matched: The detection engine uses a pattern matching algorithm to inspect the URL content, if one security rule matches with a specific pattern, the detection engine rejects the HTTP request.

• 6 No Rule matched: After applying all security rules, no rule matches with the analyzed content. where Xl," . , Xn are the values taken by the random vari ables Xl, ... , Xn that are assumed conditionally independent relative to the category c (SQLi, Legitimate). Each random variable gives an information about a pattern type in a dis sected header. In this model, all random variables are binary:

X i = 1 if pattern of type i noted pa i is present, otherwise X i = O.
Consequently, each random variable X i = 1 follows a Bernoulli distribution with parameter P i = p(pa i )'

D. The Machine Learning Model

A naive Bayesian model, which we have adopted in our classifier, is a simple classification scheme that estimates the class-conditional probability by assuming that features are conditionally independent. In fact, we have a binary classification problem of identifying the normal http stream as legitimate requests class and malicious stream as SQLi requests class. According to the Bayes theorem [START_REF] Robert | Le choix Baysien. Principes et pratiques[END_REF] and the total probabilities theorem, for a vector :t = (Xl"", Xn), the probability to belong to the class c is defined as follows:

(1)

Using the theorem of the total probabilities, we deduce:

p(C = c) -p(x =:t jC=c) L p(C=c) . p(x =:t jC=c) (2)
cE{SQLi,Leg} 1) Cost Evaluation: Works of [12] on cost-sensitive evalu ation measures consist of evaluating the false positives effect on the total cost, in term of time wasting by users to delete spams. But in our case, the cost is the impact of false negatives on the trust of users granted to our classifier.

2) The False Negatives Effect: A false negative is mistak enly classifying an SQLi attack as a legitimate content, and a false positive is a legitimate content mistakenly classified as an SQLi attack. In our model, the cost of a false negative is much higher than the cost of a false positive. Indeed, wasting time in analyzing legitimate requests is more acceptable than passing a malicious code to the Web application. The two error types are defined as such follows:

• Classifying an SQLi attack as legitimate content:

(SQLi -7 Leg)

• Classifying a legitimate content as an SQLi attacks:

(Leg -7 SQLi)

In our classifier model, the first error is more serious than the second error. To illustrate this idea, we introduce the parameter A, as its objective is to give more importance to the first error by assuming that Leg -7 SQ Li is A times more costly than SQLi -7 Leg.

3) Classification criteria: According to the above two error types, the selection criteria is as follows:

The content of a header 55 is classified as legitimate if and The parameters defined above do not take into consideration the notion of weight for the two error types introduced in the previous paragraph. This leads us to introduce the weighted accuracy (Wacc) and weighted error (Werr = 1-Wacc). We assumed that SQLi -7 Leg is A times more devastating for our system than Leg -7 SQLi. To make accuracy and error rate sensitive to this cost, we should treat each SQL injection as if it was A inputs; when an SQL injection is misclassified, this counts as A errors; and when it is classified correctly, this counts as A successes

W

A nsqli--+sqli + nleg--+leg ac e = _----'c_----'c __ ----"_---"- A Nsqli + Nleg W _ A _ n ..c:. s ,,: ql -c i--+ :-: :--l e ", g ,-+ _ n -:-l -:-e ,,-g --+ ,--,-s ,,: ql ..:. i err =

A Nsqli + Nleg [START_REF] Ristic | ModSecurity Handbook: The Complete Guide to the Popular Open Source Web Application Firewall[END_REF] To have a precise idea of the filter's performance, we compare it to a non-filtered system in which all requests are considered as legitimate. We introduce the definition of the base-line weighted error and the base-line weighted accuracy (respectively noted Wacc b and Werr b ) which are defined as follows:

W b

A Nleg ac e = ��----"-:-::--

- A Nsqli + Nleg W b Nsqli err = ----'---- A Nsq1i + N1eg (8)
The TC R (Total Cost Ratio) value measures the performance of a machine learning classifier to the same environment without a classifier. In the case where the TC R value is negligible, the best approach is to not use a classifier and to send all requests to the rules based detection engine. An effective filter which could be used in real environments should have a TC R value higher than 1.

The TC R formula is defined as follows:

TC R = Werr b = Nsqli Werr A nleg--+sqli + nsqli--+leg V. RESULTS EVALUATION (9)
We have collected a training data set from the framework previously presented in this paper. The training set is a mixture of SQIi attacks and legitimate requests shown in the below table with different proportions (80-20)%,(66-34)%,(50-50)%, then we varied, for each scenario, the A parameter and finally calculated values of Wacc, Wacc b and TCR.

By increasing the a (threshold) value from 99% to 50%, we have an increase in number of false positives, which means that the rules engine will analyze legitimate requests, but at the same time, the evaluation has shown an increase in TCR value. However, in practice, false positives (legitimate requests classified as SQLi attacks) will be forwarded to the rules engine that will apply all security rules. This will decrease the overall performances of our system. To find a compromise between lower number of false positives and false negatives, the TCR values should be medium values. Acceptable values of TCR are related to a threshold values higher than 50% and lower than 90%.

Data (%)

.\ .66-34 5,00 0,00

"

Fig. 4: Total Cost Ratio VI. COMPARATIVE ANALYSIS Security Rules based WAFs do not have a module that predicts SQL injection attacks. In this case, our approach will improve the inspection performance because the classifier module will not forward legitimate traffic to the detection engine. Only suspected requests are deeply inspected by apply ing pattern matching algorithms. In [START_REF] Komiya | CJassification of malicious web code by machine learning[END_REF], Authors worked on Classification of Malicious Web Code by Machine Learning. They implemented and evaluated two classifiers both SQLIAs and XSS which can use TF-IDF method for weight calculation, and three machine learning approach among SVM, Naive Bayes, k-Nearest Neighbor Algorithm. They obtain good pre cision values for their classifier (99.16 %) by using SVM with Gaussian Kernel. Our classifiers obtains 97.6% by using Bayesian algorithm. However, they do not provide a solution to handle false negatives and false positives. On the other hand, the use of SVM algorithms with Gaussian Kernel may require significant CPU resources and decreases the classifier performances in multi-gigabits rates networks.

20i4 international Conference on information Assurance and Security (iAS)

VII. CONCLUSION

In this paper, we focused on the problem of detecting com plex SQL injections. We proposed a novel approach to dissect HTTP requests in order to cover most evasion techniques and improve security rules management process. We also provided an Injection Prevention System architecture which includes a machine learning classifier. Based on the TCR results, we have shown the effectiveness of the classifier by tuning its values in order to reduce false negatives. We were also able to show that false positives did not impact the overall performances of our system. A key element of future work is to apply the same approach in order to develop an anti XSS and SQL attacks solution.

Fig. 1 :

 1 Fig. 1: Raw HTTP Request b) HTTP Request Dissection: Our dissection module is able to recognize request's components (headers and the body)

Fig, 2 :Fig, 3 :

 23 Fig, 2: Detection Engine Architecture

4 )

 4 only if: p(C = LeglX = 55) > A. p(C = SQLilX = 55) (3) given that p(C = LeglX = 55) + p(C = SQLilX = 55) = 1, the selection criteria becomes as follows: Where : p( C = LeglX = 55) Method and parameters evaluation: In this section, we define the parameters that allow us to evaluate our filter. To this end, two evaluation parameters are used: accuracy (Acc) and error (Err = 1 -Acc) [12]. They are defined as follows: Acc Err where: nsqli--+sqli + nleg--+leg Nsqli + Nleg nsqli--+leg + nsqli--+sqli Nsqli + Nleg • Nsqli = nsqli--+leg + nsqli--+sqli • Nleg = nleg--+leg + nleg--+sqli(6)• n y --+ z denotes the number of patterns of class y that are mistakenly classified in class z.

  table shows the result of our features vector selection.

	Variables	Tokens (key words)
	xl	SELECT
	x2	UNION
	x3	UPDAT E
	..	Other SQL keywords
	xn
		Other Special Symbols
	xm	UNHEX
	..	Other Evasion Keywords
	x45	% (percent)

* (Asterisk)

..

TABLE I : Features Table 3 )

 IFeatures3 Request

representation: Each request's header is char acterized by a vector :t defined by :t = (Xl"", xn)

TABLE II :

 II Costs Table

	TCR
	30,00
	25,00
	20,00
	15,00
	• 50-SO
	10,00

International Conference on Information Assurance and Security (lAS)

20i4 international Conference on information Assurance and Security (iAS)