Francois James
email: francois.james@univ-orleans.fr

Nicolas Vauchelet
email: vauchelet@ann.jussieu.fr

Numerical methods for one-dimensional aggregation equations

Keywords: aggregation equation, duality solutions, finite volume schemes, asymptotic preserving schemes, weak measure solutions, hydrodynamical limit. 2010 AMS classifications. 35B40, 35D30, 35L60, 35Q92, 65M08

We focus in this work in the numerical discretization of the one dimensional aggregation equation

, in the attractive case. Finite time blow up of smooth initial data occurs for potential W having a Lipschitz singularity at the origin. A numerical discretization is proposed for which the convergence towards duality solutions of the aggregation equation is proved. It relies on a careful choice of the discretized macroscopic velocity v in order to give a sense to the product vρ. Moreover, using the same idea, we propose an asymptotic preserving scheme for a kinetic system in hyperbolic scaling converging towards the aggregation equation in hydrodynamical limit. Finally numerical simulations are provided to illustrate the results.

Introduction

This paper is devoted to the numerical approximation of the so-called aggregation equation which writes in one space dimension

∂ t ρ + ∂ x a(W ′ * ρ)ρ = 0, (1.1)
complemented with some initial data ρ(0, x) = ρ ini (x). This nonlocal and nonlinear conservation equation is involved in many applications in physics and biology, where it describes the behaviour of a population of particles (in physical applications) or cells (in biological applications) interacting under a continuous interaction potential W . The quantity ρ denotes the density of these particles or cells. The function a is often linear (a(u) = ±u), see e.g. [START_REF] Benedetto | A kinetic equation for granular media[END_REF][START_REF] Burger | On an aggregation model with long and short range interactions[END_REF][START_REF] Li | Long time asymptotics of kinetic models of granular flows[END_REF][START_REF] Morale | An interacting particle system modelling aggregation behavior: from individuals to populations[END_REF][START_REF] Okubo | Diffusion and Ecological Problems: Modern Perspectives[END_REF][START_REF] Leverentz | Asymptotic dynamics of attractiverepulsive swarms[END_REF], but in several applications, such as pedestrian motion [START_REF] Colombo | A class of nonlocal models for pedestrian traffic[END_REF][START_REF] Crippa | Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow[END_REF] or chemotaxis (see [START_REF] James | On the hydrodynamical limit for a one dimensional kinetic model of cell aggregation by chemotaxis[END_REF] and Section 4.2 below) a specific nonlinearity has to be considered. Depending on the choice of the potential W and the function a, one can be in the repulsive or in the attractive case, the latter leading to aggregation phenomena.

In this work we focus on the case involving attractive forces. Individuals attract one another under the action of the potential W , assumed to be smooth away from 0 and bounded from below. More precisely, W satisfies the following properties:

1 Assumption 1.1 We assume that W ∈ C 1 (R \ {0}
) is a pointy attractive potential, i.e.

W ′′ = -δ 0 + w, w ∈ C 0 (R), with w L 1 (R) = w 0 < ∞, (1.2)
in the distributional sense, where δ 0 is the Dirac measure at 0.

Under Assumption 1.1, attractivity in the nonlinear case is ensured provided the function a satisfies a ∈ C 1 (R), 0 ≤ a ′ (x) ≤ α, ∀ x ∈ R. (1.3) This case has been extensively studied in the linear case a = id [START_REF] Bertozzi | Finite-time blow-up of L ∞ -weak solutions of an aggregation equation[END_REF][START_REF] Bertozzi | Blow-up in multidimensional aggregation equation with mildly singular interaction kernels[END_REF][START_REF] Bertozzi | Finite-time blow-up of solutions of an aggregation equation in R n[END_REF] and it is known that if the potential W has a Lipschitz singularity then weak solutions blow up in finite time (see e.g. [START_REF] Bertozzi | Finite-time blow-up of L ∞ -weak solutions of an aggregation equation[END_REF][START_REF] James | On the hydrodynamical limit for a one dimensional kinetic model of cell aggregation by chemotaxis[END_REF]), so that measure valued solutions arise. At the theoretical level, global in time existence has been obtained in the linear case a = id and in any space dimension by Carrillo et al. [START_REF] Carrillo | Global-intime weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF], in the space P 2 of probability measures with finite second moment, using the geometrical approach of gradient flows. In the nonlinear case, but in one space dimension, global existence of measure solutions has been obtained by completely different means in [START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF], namely thanks to the notion of duality solutions. It has also been proved in [START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF] that in the linear case a = id, both notions coincide. In the same respect, we refer to [START_REF] Bonaschi | Equivalence of gradient flows and entropy solutions for singular nonlocal interaction equations in 1D[END_REF], where gradient flow solutions are proved to be equivalent to entropy solutions of the Burgers equation, with a particular focus on the repulsive case. The key point leading to both uniqueness of solutions and equivalence between the two notions is the definition of the macroscopic velocity. In the framework of gradient flows, it is defined as the unique element with minimal norm in the subdifferential of the interaction energy associated to W , where P 2 is endowed with the Wasserstein distance (see [2,[START_REF] Carrillo | Global-intime weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF] for more details). In [START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF], the macroscopic velocity is defined using the chain rule for BV functions, and this is the viewpoint adopted for numerical analysis.

In some applications, the aggregation equation can be obtained as the hydrodynamic limit of some kinetic system. For instance, we consider here the following two velocities kinetic model in hyperbolic scaling which is used to model the so-called run-and-tumble process in bacterial chemotaxis (see [START_REF] Dolak | Kinetic models for chemotaxis: Hydrodynamic limits and spatio-temporal mechanisms[END_REF][START_REF] Filbet | Derivation of hyperbolic models for chemosensitive movement[END_REF][START_REF] James | On the hydrodynamical limit for a one dimensional kinetic model of cell aggregation by chemotaxis[END_REF]). The dynamics of the distribution function of cells f ± ε at time t, position x and with velocity ±1 is governed by

∂ t f ± ε ± ∂ x f ± ε = ± 1 ε φ(-W ′ * ρ ε )f - ε -φ(W ′ * ρ ε )f + ε . (1.4)
The left hand side corresponds to the free transport (run phase) whereas the right hand side models the velocity reorientation, called tumble process. The function φ is given by φ(x) = 1 2 (1 -a(x)) where the function a is assumed to be odd, in addition to (1.3). Existence of global in time L ∞ weak solutions for such a kinetic equation with fixed ε > 0 is well-known (see e.g. [START_REF] Bournaveas | Global existence for a kinetic model of chemotaxis via dispersion and Strichartz estimates[END_REF][START_REF] Vauchelet | Numerical simulation of a kinetic model for chemotaxis[END_REF]). However the bound is not uniform in ε and it actually turns out that limit when ε goes to zero of the total density ρ = lim ε→0 (f + ε + f - ε ) solves the aggregation equation (1.1), see Section 3 below for details.

Therefore an interesting issue consists in providing a numerical scheme for the kinetic system (1.4) which allows to recover the asymptotic limit when ε → 0. Such schemes are usually called asymptotic preserving (AP) [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF]. They are of great interest for kinetic equations since letting ε → 0 with the mesh size and time step fixed, the scheme becomes a scheme for the macroscopic limit (see e.g. [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]). In other words, AP schemes allow a numerical discretization whose time step is not constrained by some constant depending on ε. We refer to [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF] for a review on AP schemes.

The aim of this work is precisely to design numerical methods for (1.1) and (1.4) that are able to capture the measure solutions after blow-up. The main difficulty is that after blow up the velocity a(W ′ * ρ) is discontinuous, so that the definition of the flux has to be considered with great care. Following the principle that holds at the continuous level, the numerical velocity is obtained thanks a careful discretization of the Vol'pert calculus for BV functions. We emphasize that the numerical solution may depend upon the way of discretizing the velocity. For equation (1.1) we work directly on the definition of a(W ′ * ρ), for the kinetic model, the discretization is defined through the right-hand side of equation (1.4). The final scheme is obtained by a splitting technique, as in for instance [START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF][START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF], which is in this particular case very easy to implement. As we shall see, this is not completely satisfactory in the sense that the fluxes at equilibrium are not correctly computed. A more sophisticated technique consists in using well-balanced schemes [START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF][START_REF] Gosse | Computing qualitatively correct approximations of balance laws[END_REF] as it has been successfully used for chemotaxis models in [START_REF] Gosse | A well-balanced scheme for kinetic models of chemotaxis derived from onedimensional local forward-backward problems[END_REF][START_REF] Gosse | Asymptotic-preserving and well-balanced schemes for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF][START_REF] Gosse | Maxwellian decay for well-balanced approximations of a super-characteristic chemotaxis model[END_REF] for smooth solutions. However, it is not clear that such schemes allow to recover the solutions after blow up.

The outline of the paper is as follows. Section 2 is devoted to the aggregation equation (1.1). After recalling existence and uniqueness result for this system, we provide a numerical scheme and prove its convergence. In Section 3, we consider the kinetic equation (1.4). We first establish the rigorous derivation of the aggregation equation thanks to a hyperbolic limit ε → 0. Then we propose an asymptotic preserving scheme and prove its convergence. Finally, Section 4 is devoted to some numerical simulations.

Part of these results were announced in [START_REF] James | Numerical simulation of a hyperbolic model for chemotaxis after blow-up[END_REF].

2 Aggregation equation

Existence of duality solutions

For Y and Z two metric spaces we denote C(Y, Z) the set of continuous functions from Y to Z. We denote C 0 (R) the set of continuous functions that vanish at infinity and C c (R) the set of those with compact support. Let M b (R) be the set of bounded Radon measures and by P 1 (R) the set of positive measure in M b (R) such that R |x|dµ(x) < ∞. From now on, the space M b (R N ) is always endowed with the weak topology σ(M b , C 0 ). We denote

S M := C([0, T ]; M b (R N ) -σ(M b , C 0 )).
Duality solutions have been introduced in [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF] to solve scalar conservation laws with discontinuous coefficients. More precisely, it gives sense to measure valued solutions of the scalar conservation law

∂ t ρ(t, x) + ∂ x (b(t, x)ρ(t, x)) = 0, where b ∈ L ∞ ((0, T ) × R) satisfies the so-called one-sided Lipschitz condition ∂ x b(t, .) ≤ β(t) for β ∈ L 1 (0, T ), in the distributional sense. (2.1)
This key point suggests that the velocity field should be compressive. We refer to [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF] for the precise definition and general properties of these solutions.

Let us first define a notion of duality solution for the aggreagation equation (1.1) in the spirit of [START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF][START_REF] James | Chemotaxis: from kinetic equations to aggregate dynamics[END_REF]:

Definition 2.1 We say that ρ ∈ C([0, T ]; M b (R)) is a duality solution to (1.1) if there exists a ρ ∈ L ∞ ((0, T ) × R) and α ∈ L 1 loc (0, T ) satisfying ∂ x a ρ ≤ α in D ′ ((0, T ) × R), such that for all 0 < t 1 < t 2 < T , ∂ t ρ + ∂ x ( a ρ ρ) = 0
in the sense of duality on (t 1 , t 2 ), and a ρ = a(W ′ * ρ) a.e.

From now on, we denote by A the antiderivative of a such that A(0) = 0. Using the chain rule, a natural definition of the flux is

J := -∂ x A(W ′ * ρ) + a(W ′ * ρ)(w * ρ).
(2.2)

In fact, a formal computation shows that

-∂ x A(W ′ * ρ) = -a(W ′ * ρ)(W ′′ * ρ) = a(W ′ * ρ)(ρ -w * ρ),
where we use (1.2) for the last identity.

Then we are in position to state the existence and uniqueness result of [START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF].

Theorem 2.2 ([27], Theorem 3.9) Let us assume that ρ ini is given in P 1 (R). Under Assumption 1.1 on the potential W and (1.3) for the nonlinear function a, for all T > 0 there exists a unique duality solution ρ of (1.1) in the sense of Definition 2.1 with ρ ≥ 0, ρ(t) ∈ P 1 (R) for t ∈ (0, T ) and which satisfies in the distributional sense:

∂ t ρ + ∂ x J = 0, (2.3)
where J is defined in (2.2). Moreover, there exists a, called universal representative, such that a = a(W ′ * ρ) a.e. Then ρ = X # ρ ini , where X is the Filippov flow associated to the velocity a.

Numerical discretization

Let us consider a uniform space discretization with step δx and denote by δt the time step; then t n = nδt and x i = x 0 + iδx, i ∈ Z. We assume that the initial datum ρ ini is given in P 1 (R).

For n ∈ N, we assume to have computed an approximation (ρ n i ) i∈Z of (ρ(t n , x i )) i∈Z , we denote by (S n i ) i∈Z an approximation of (W * ρ(t n , x i )) i∈Z and by (ν n i ) i∈Z an approximation of (w * ρ(t n , x i )) i∈Z . Let us denote λ = δt/δx and M the total mass of the system, M = |ρ ini |(R). We obtain an approximation of ρ(t n+1 , x i ) denoted ρ n+1 i by using the following Lax-Friedrichs discretization of equation (2.3)-(2.2):

ρ n+1 i = ρ n i -λ(J n i+1/2 -J n i-1/2 ) + λ 2 c(ρ n i+1 -2ρ n i + ρ n i-1 ) (2.4) J n i+1/2 = - A(∂ x S n i+1 ) -A(∂ x S n i ) δx + a n i+1/2 ν n i+1 + ν n i 2 , (2.5)
where we have defined

c := max x∈[-M (1+w 0 ),M (1+w 0 )] |a(x)|. (2.6)
In this scheme, we use the discretization

∂ x S n i+1 = S n i+2 -S n i 2δx , (2.7)
and the approximation

a n i+1/2 =      0 if ∂ x S n i+1 = ∂ x S n i , A(∂ x S n i+1 ) -A(∂ x S n i ) ∂ x S n i+1 -∂ x S n i otherwise.
(2.8)

We need now a scheme for S n i . From Assumption 1.1, we deduce by taking the convolution of (1.2) with ρ that -W ′′ * ρ + w * ρ = ρ. This equation is discretized by using a standard finite difference scheme:

- S n i+1 -2S n i + S n i-1 δx 2 + ν n i = ρ n i .
(2.9)

For the computation of (ν n i ) i , there are multiple ways; here we propose to use a piecewise constant approximation for ρ on each interval [x i , x i+1 ), so that

ν n i = k∈Z x k+1 x k ρ n i w(x i -y) dy,
which can be rewritten

ν n i = k∈Z ρ n k w ki , w ki = x i -x k x i -x k+1 w(z) dz = (i-k)δx (i-1-k)δx w(z) dz.
(2.10)

We turn back now to the discrete macroscopic flux (2.5), and notice that using (2.8) we have

J n i+1/2 = a n i+1/2 - ∂ x S n i+1 -∂ x S n i δx + ν n i+1 + ν n i 2 .
Using then (2.7) and (2.9) we obtain the following result.

Lemma 2.3 With the choice (2.8) for the discretization of the macroscopic velocity a n i+1/2 , the macroscopic flux (2.5) rewrites

J n i+1/2 = a n i+1/2 ρ n i+1 + ρ n i 2 . ( 2
ρ n+1 i = ρ n i 1 -λc + λ 4 (a n i-1/2 -a n i+1/2 ) + λ 2 c + a n i-1/2 2 ρ n i-1 + λ 2 c - a n i+1/2 2 ρ n i+1 . (2.12)
Remark 2. [START_REF] Bertozzi | Finite-time blow-up of L ∞ -weak solutions of an aggregation equation[END_REF] The choice of the discretization (2.8) for the macroscopic velocity can be seen as a consequence of the chain rule (or Vol'pert calculus) for BV functions [START_REF]pert, The spaces BV and quasilinear equations[END_REF] (see also remark 3.98 of [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]): for a BV function u, the fonction a V defining the chain rule

∂ x A(u) = a V ∂ x u is constructed by a V (x) = 1 0 a(tu 1 (x) + (1 -t)u 2 (x)) dt,
where

(u 1 , u 2 ) =        (u, u) if x ∈ R \ S u , (u + , u -) if x ∈ J u ,
arbitrary elsewhere.

We have denoted by S u the set of x ∈ R where u does not admit an approximate limit and by J u ⊂ S u the set of jump points. Applying that to u = ∂ x S, we obtain,

a V (x) =          a(∂ x S(x)) if x ∈ R \ S u , A(∂ x S(x + )) -A(∂ x S(x -)) ∂ x S(x + ) -∂ x S(x -) if x ∈ J u ,
arbitrary elsewhere.

Numerical analysis

In this subsection, we prove the convergence of the numerical scheme defined in (2.4)-(2.8) towards the unique duality solution of Theorem 2.2. We first state a Lemma which proves a CFL-like condition for the scheme:

Lemma 2.5 Let us assume that (1.3) holds and that the condition

λ := δt δx ≤ 1 2c , (2.13)
is satisfied with c defined in (2.6). Let us assume that ρ ini ∈ P 1 (R) is given, compactly supported and nonnegative, and we define

ρ 0 i = 1 δx x i+1 x i ρ ini (dx) ≥ 0.
Then for all i and n ∈ N, the sequences computed thanks to the scheme defined in (2.4)-(2.9) satisfy

ρ n i ≥ 0, |a n i+1/2 | ≤ c. Proof.
• Preliminaries. Let us define M n i = δx j≤i ρ n j . Since the scheme (2.4) is conservative, the total mass of the system is conserved, then

M n ∞ = M 0 ∞ = M . Clearly, ρ n i = (M n i -M n i-1 )/δx and from (2.11) we have J n i+1/2 = a n i+1/2 (M n i+1 -M n i-1
)/(2δx). Then we deduce from (2.4) that

M n+1 i = (1 -λc)M n i + λ 2 c -a n i+1/2 M n i-1 + λ 2 c + a n i+1/2 M n i+1 . (2.14)
Thus M n+1 i can be written as a convex combination of M n i-1 , M n i and M n i+1 provided the CFL condition is satisfied and |a n i+1/2 | ≤ c, which we prove below by induction. We first establish some estimates. By definition of ν n i in (2.10), we deduce from (1.2) that for all i ∈ Z,

δx j≤i ν n j ≤ δx k∈Z ρ n k w 0 = M w 0 . (2.15)
Moreover, from the definition of ∂ x S n i in (2.7) and using equation (2.9), we deduce

- ∂ x S n i+1 -∂ x S n i δx + ν n i+1 + ν n i 2 = ρ n i+1 + ρ n i 2 .
Summing this latter equation over i, we obtain

∂ x S n i+1 = - 1 2 M n i+1 + M n i -δx 2 j≤i ν n j + ν n i+1 -ν n 0 .
(2.16)

• Induction on n. We are now in position to prove the lemma by an induction on n. For n = 0, by construction of the initial data, we have ρ 0 i ≥ 0. Then for all i ∈ Z, we have 0 ≤ M 0 i ≤ M and with (2.16) and (2.15) we deduce that

|∂ x S 0 i+1 | ≤ M (1 + w 0 ), for all i ∈ Z.
Futhermore, since we have

A(∂ x S 0 i+1 ) -A(∂ x S 0 i ) ∂ x S 0 i+1 -∂ x S 0 i = a(θ 0 i ), θ 0 i ∈ (∂ x S 0 i , ∂ x S 0 i+1 ) ⊂ (-M (1 + w 0 ), M (1 + w 0 )),
we deduce with (2.8) that |a 0 i+1/2 | ≤ c, which proves the result for n = 0. Let us assume that ρ n i ≥ 0 and |a n i+1/2 | ≤ c, for some n ∈ N. From condition (2.13) and the induction assumption |a n i+1/2 | ≤ c, we deduce that in the scheme (2.12), all the coefficients in front of ρ n i-1 , ρ n i and ρ n i+1 are nonnegative. Thus ρ n+1 i ≥ 0 for all i. Moreover, we have clearly by definition that 0 ≤ M n i ≤ M . Then, from the condition (2.13) and induction assump-

tion |a n i+1/2 | ≤ c, we deduce with (2.14) that M n+1 i is a convex combination of M n i+1 , M n i and M n i-1 . Then 0 ≤ M n+1 i ≤ M
. Thus, as above, using (2.16) with n + 1 instead of n, we have

|∂ x S n+1 i+1 | ≤ M (1 + w 0 ), which implies |a n+1 i+1/2 | ≤ c.
For any given sequence (u n i ) i , we define the corresponding piecewise constant reconstruction

u δ (t, x) = n∈N i∈Z u n i 1 [nδt,(n+1)δt)×[x i ,x i+1 ) (t, x).
(2.17)

Then we have the following convergence result:

Theorem 2.6 Let us assume that ρ ini ∈ P 1 (R) is given, compactly supported and nonnegative and define

ρ 0 i = 1 δx x i+1 x i ρ ini (dx) ≥ 0. Under assumption (1.3), if (2.13
) is satisfied, then the discretization ρ δ converges in S M towards the unique duality solution ρ of Theorem 2.2 as δt and δx go to 0.

Proof. The proof is divided into several steps.

(i) Estimates. Recall that from the proof of Lemma 2.5, we have 0

≤ M n i ≤ M . Therefore, since 0 ≤ ρ n i = (M n i -M n i-1 )/δx, equation (2.12) implies a BV (R) estimate on (M n i ) i , provided (2.13) is satisfied. More precisely the scheme is TVD for the sequence (M n i ) i . (ii) Convergence. Standard techniques imply a L ∞ ∩BV ((0, T )×R) estimate on the function M δ defined from (M n i ) by (2.17). It implies the convergence, up to a subsequence, of M δ in L 1 loc (R + × R) towards a function M ∈ L ∞ ∩ BV ((0, T ) × R)
when δt and δx go to 0 and satisfy (2.13).

Let us define

ρ = ∂ x M ∈ L ∞ ((0, T ); M b (R)). Obviously, noting that ρ n i = (M n i -M n i-1
)/δx, we deduce that ρ is the limit in S M of ρ δ . By definition (2.10), we have that ν δ = w * ρ δ . Therefore, the sequence (ν δ ) δ converges, up to a subsequence, towards ν := w * ρ for a.e. t > 0 and x ∈ R.

From (2.16), we deduce that we have the same bound on the sequence (∂ x S n i ) i,n as on (M n i ) i,n . We conclude that the sequence (∂

x S n i ) i,n is bounded in L ∞ ∩ BV ((0, T ) × R).
As above, we get the convergence, up to a subsequence, in

L 1 loc (R + × R) of ∂ x S δ towards a function ∂ x S belonging to L ∞ ∩ BV ((0, T ) × R)
as δt and δx go to 0 and satisfy (2.13). By definition of ∂ x S δ , we have the strong convergence up to a subsequence in

L 1 loc (R + , W 1,1 loc (R)) of S δ towards S.
(iii) Passing to the limit. Passing to the limit in the equation (2.9) we deduce that S and w satisfy in the weak sense the equation

-∂ xx S + w = ρ.
Moreover, from Lemma 2.5, we deduce that the sequence (a δ ) δ is bounded in L ∞ , thus we can extract a subsequence converging in L ∞ -weak * towards a. From the L 1 loc convergence of (∂ x S δ ) δ , we deduce that a = a(∂ x S) a.e. Then, from (2.5), we have the convergence in the sense of distribution of J δ towards J = -∂ x (A(∂ x S)) + aw a.e. Finally, taking the limit in the distributional sense of equation (2.4) we deduce that ρ is a solution in the sense of distribution of (2.3)-(2.2). By uniqueness of this solution, we deduce that ρ is the unique duality solution of Theorem 2.2.

Finally, we notice that, as in the continuous case (see [START_REF] James | Chemotaxis: from kinetic equations to aggregate dynamics[END_REF]), the nonnegativity of the density ρ allows to ensure an one-sided Lipschitz condition on the discretized macroscopic velocity.

Proposition 2.7 With the notations and assumptions of Theorem 2.6, the discrete macroscopic velocity in (2.8) satisfies the discrete one-sided Lipschitz condition:

1 δx a n i+1/2 -a n i-1/2 ≤ 2M α w ∞ ,
where M is the total mass of the system and α is defined in (1.3).

Proof. From definition (2.8) we have, applying the mean value Theorem:

1 δx a n i+1/2 -a n i-1/2 = 1 δx a(θ n i+1/2 ) -a(θ n i-1/2 ) = a ′ (γ n i ) δx θ n i+1/2 -θ n i-1/2 , where θ n i+1/2 ∈ (∂ x S n i , ∂ x S n i+1
) and γ n i ∈ (θ n i-1/2 , θ n i+1/2 ) (where the notation (µ, ν) stands for the interval (ν, µ) when ν < µ). Then, using assumption (1.3),

1 δx a n i+1/2 -a n i-1/2 ≤ a ′ (γ n i ) δx max 0, ∂ x S n i+1 -∂ x S n i , ∂ x S n i -∂ x S n i-1 , ∂ x S n i+1 -∂ x S n i-1 .
From the definition (2.7) and with (2.9), we have

∂ x S n i+1 -∂ x S n i δx = ν n i+2 + ν n i 2 - ρ n i+1 + ρ n i 2 ≤ ν n i+2 + ν n i 2 ,
and

∂ x S n i+1 -∂ x S n i-1 δx = ν n i+2 + 2ν n i + ν n i-1 2 - ρ n i+1 + 2ρ n i + ρ n i-1 2 ≤ ν n i+2 + 2ν n i + ν n i-1 2 ,
where we use the nonnegativity on the sequence (ρ n i ) i,n . By definition (2.10), the sequence

(ν n i ) i,n is uniformly bounded in L ∞ by w ∞ M . We deduce that (∂ x S n i+1 -∂ x S n i )/δx and (∂ x S n i+1 -∂ x S n i-1
)/δx are uniformly bounded from above by w ∞ M . Using moreover assumption (1.3) allows to conclude the proof.

Remark 2.8 In some applications, we have w = W . In this case, we prefer to set ν n i = S n i instead of (2.10). The sequence (S n i ) i,n is then entirely determined by solving (2.9). The convergence result in Theorem 2.6 still holds true, with a straightforward adaptation of the proof.

Asymptotic preserving scheme

In this section, we consider an asymptotic preserving scheme allowing to recover the numerical discretization (2.4)-(2.5) from a kinetic model (1.4).

Asymptotic preserving (AP) schemes have been widely developed since the 90s for a wide range of time-dependent kinetic and hyperbolic equations. The basic idea is to develop a numerical discretization that preserves the asymptotic limits from the microscopic to the macroscopic models [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF]. Moreover, in the definition of [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF], an AP scheme should be implemented explicitely (or at least more efficiently than using a Newton type solvers for nonlinear algebraic systems).

Hydrodynamic limit

As already mentioned, aggregation equation (1.1) can be derived as the hydrodynamic limit of the kinetic system (1.4). We recall that the function φ satisfies

φ(x) = 1 2 (1 -a(x))
with a an odd function satisfying (1.3) and a ∞ ≤ 1.

(3.1)

We introduce the potential S ε = W * ρ ε , which, due to (1.2) is a weak solution to

-∂ xx S ε + w * ρ ε = ρ ε . (3.2)
Using the fact that φ(x) + φ(-x) = 1, we can rewrite these equations as

∂ t f ± ε ± ∂ x f ± ε = 1 ε φ(∓∂ x S ε )ρ ε -f ± ε , (3.3)
Remark 3.1 For the sake of clarity, and since our examples are limited to this case, we consider only the two-velocities model (3.3). However the results can be adapted to more general models with continuous and bounded velocities.

As above, we denote by A an antiderivative of a and by Φ an antiderivative of φ, so that, from (3.2),

Π ± ε := φ(±∂ x S ε )ρ ε = ∓∂ x (Φ(±∂ x S ε )) + (w * ρ ε )φ(±∂ x S ε ). (3.4)
Now we recall briefly the derivation of the aggregation equation (1.1) from the two velocities kinetic system (3.3). We first notice that the momentum equations obtained by adding and substracting the two equations in (3.3) are given by

∂ t ρ ε + ∂ x J ε = 0, ( 3.5)
∂ t J ε + ∂ x ρ ε = 1 ε Π - ε -Π + ε -J ε , (3.6)
where

ρ ε = f + ε + f - ε and J ε = f + ε -f + ε . From (3.1
) and (3.4), we deduce

Π - ε -Π + ε = a(∂ x S ε )ρ ε = -∂ x A(∂ x S ε ) + (w * ρ ε )a(∂ x S ε ). (3.7)
Letting formally ε → 0 in (3.6), we get

J ε → J 0 := Π - 0 -Π + 0 = -∂ x A(∂ x S 0 ) + (w * ρ 0 )a(∂ x S 0 ).
Inserting in (3.5), we recover the aggregation equation

∂ t ρ 0 + ∂ x J 0 = 0, J 0 = -∂ x A(∂ x S 0 ) + (w * ρ 0 )a(∂ x S 0 ), (3.8)
where S 0 = W * ρ 0 . The preceding formal computations can be made rigorous. Following the lines of proof of Theorem 3.10 of [START_REF] James | Chemotaxis: from kinetic equations to aggregate dynamics[END_REF], we can indeed state the following theorem. Theorem 3.2 Assume φ is given by (3.1) and W as in (1.2). Let f ini ≥ 0 be given such that

ρ ini := f ini,+ + f ini,-belongs to P 1 (R). Let f ± ε be solutions to (3.3)-(3.
2). Then, as ε → 0, f ± ε converge in the following sense:

ρ ε := f + ε + f - ε ⇀ ρ in S M := C([0, T ]; M b (R) -σ(M b , C 0 )),
where ρ is the unique duality solution of Theorem 2.2.

An AP numerical scheme

As above, we consider a time discretization of step δt and a uniform space discretization of step δx. We denote f n,± h,ε an approximation of f ± ε at time t n = nδt, for n ∈ N. Before giving the numerical scheme, we explain formally the main idea which is based on the following time splitting argument:

• Assuming that approximations f n,± h,ε and S n h,ε of f ± ε and S ε are known at time t n , we set

ρ n h,ε = f n,+ h,ε + f n,- h,ε .
We have now everything at hand to compute Π n,± h,ε from (3.4):

Π n,± h,ε := ∓∂ x Φ(±∂ x S n h,ε ) + (w * ρ n h,ε )φ(±∂ x S n h,ε ). (3.9)
We solve in this first step, during a time step δt, the relaxation equation

∂ t f ± h,ε = 1 ε (Π ∓ h,ε -f ± h,ε ), (3.10
= S n h,ε . Therefore Π n,± h,ε is constant during this time step: Π n+1/2,± h,ε = Π n,± h,ε
, so that we can solve exactly equation (3.10) by

f n+1/2,± h,ε = e -δt/ε (f n,± h,ε -Π n,∓ h,ε ) + Π n,∓ h,ε . (3.11)
• In a second step, we discretize during a time step δt the free transport equation:

∂ t f n,± h,ε ± ∂ x f n,± h,ε = 0.
Denoting by D x some discrete derivative with respect to x, we obtain

f n+1,± h,ε = f n+1/2,± h,ε ∓ δtD x f n+1/2,± h,ε . (3.12) Then we compute ρ n+1 h,ε = f n+1,+ h,ε + f n+1,- h,ε
and solve (3.2) to obtain S n+1 h,ε . We can start then a new iteration.

The main advantage of this method is that the small parameter ε is taken into account only in the first step. Letting ε → 0, we deduce easily from (3.11) that at the limit ε → 0, we have

f n+1/2,± h,ε → f n+1/2,± h,0 = Π n+1/2,± 0 = Π n,± 0 , since as explained above ρ n+1/2 h,ε = ρ n h,ε and S n+1/2 h,ε = S n h,ε
. Moreover, with (3.9) we obtain

Π n,± h,0 := -∂ x Φ(±∂ x S n h,0 ) + (w * ρ n h,0 )φ(±∂ x S n h,0 ).
Then by applying the first step (3.12), we have

f n+1,± h,0 = Π n+1/2,± h,0 ∓ δtD x Π n+1/2,± 0
.

Summing, we deduce, using the notations

ρ h,0 = f + h,0 + f - h,0 and J h,0 = f + h,0 -f - h,0 that ρ n+1 h,0 = ρ n+1/2 h,0 -δtD x J n+1/2 h,0
.

Moreover, we have ρ

n+1/2 0 = ρ n 0 and J n+1/2 0 = J n 0 . Then, ρ n+1 h,0 = ρ n h,0 -δtD x J n h,0 ,
which is an explicit in time discretization of the conservation equation (3.8). Numerical discretization. We recall the time and space discretization t n = nδt, n ∈ N and x i = iδx, i ∈ Z. As above, (ρ n i ) i,n , (S n i ) i,n and (ν n i ) i,n are approximations of resp. (ρ(t n , x i )) i,n , (S(t n , x i )) i,n and (W * ρ(t n , x i )) i,n . Moreover, we denote by (f n,± i ) i,n an approximation of (f ± ε (t n , x i )) i,j,n and by (Π n,± i ) i,n an approximation of (Π ± ε (t n , x i )) i,n defined in (3.4). Assuming (f n,± i ) i,n are known for some n ∈ N, we compute the approximated density

ρ n i = f n,+ i + f n,- i , J n i = f n,+ i -f n,- i .
Then we compute (ν n i ) i as in previous Section thanks to (2.10) and the macroscopic potential (S n i ) i by solving

- S n i+1 -2S n i + S n i-1 δx 2 + ν n i = ρ n i . (3.13)
We have seen in the previous Section that the flux and therefore the corresponding macroscopic velocity should be defined with care. Then we use the following discretization of the macroscopic velocity

a n i =        0, if ∂ x S n i+1/2 = ∂ x S n i-1/2 , A(∂ x S n i+1/2 ) -A(∂ x S n i-1/2 ) ∂ x S n i+1/2 -∂ x S n i-1/2 , otherwise. (3.14)
In this expression, we use the notation ∂ x S n i+1/2 = (S n i+1 -S n i )/δx. We define then the approximations φ n,± i of φ(±∂ x S ε (t n , x i )) by

φ n,± i = 1 2 1 ± a n i . (3.15)
Then the quantity Π defined in (3.4) is approximated by

Π n,± i = φ n,± i ρ n i . (3.16)
Using (3.14) and (3.13), we also have

Π n,± i = 1 2 (ρ n i ± a n i ν n i ) ∓ 1 2 A(∂ x S n i+1/2 ) -A(∂ x S n i-1/2 ) δx . (3.17)
Then we compute

f n+1/2,± i = e -δt/ε f n,± i + (1 -e -δt/ε )Π n,± i . ( 3
+ φ n,- i = 1
. Then, we obtain (f n+1 ij ) i,j by, for instance, applying an upwind discretization for the step 2:

   f n+1,+ i = f n+1/2,+ i -λ(f n+1/2,+ i -f n+1/2,+ i-1
),

f n+1,- i = f n+1/2,- i + λ(f n+1/2,- i+1 -f n+1/2,- i ), (3.19)
where we recall that λ = δt/δx.

The following theorem states the AP property of the scheme defined above.

Theorem 3.3 Let us assume that the initial distributions f ini,± are nonnegative and such that ρ ini := f ini,+ + f ini,-belongs to P 1 (R). Let us assume that φ satisfies (3.1). Let us consider the sequence (f n,± i ) i,n computed thanks to (3.13)- (3.19). Then, under the CFL condition,

λ := δt δx ≤ 1, (3.20)
as ε → 0, the sequence (ρ n i ) i,n := (f n,+ i + f n,- i ) i,n converges weakly, up to a subsequence, towards the sequence ( ρ n i ) i,n , computed by a Lax-Friedrichs discretization as in (2.4)-(2.5) of the equation

∂ t ρ + ∂ x J = 0, J = -∂ x A(∂ x S) + a(∂ x S) S.
Proof. The proof is divided into several steps. (i) Nonnegativity. Assume f n,± i ≥ 0 for all i ∈ Z. By the assumption a ∞ ≤ 1, we have that the function φ is nonnegative and therefore φ n,± i defined in (3.15)

f n+1,± i is a convex combination of f n+1/2,± i , f n+1/2,± i-1 and f n+1/2,± i+1 . Then f n+1,± i ≥ 0.
(ii) Mass conservation. We recall that we have from (3.18)

that ρ n+1/2 i = ρ n i . Summing (3.19) over i ∈ Z, we deduce that δx i∈Z ρ n+1 i = δx i∈Z ρ n+1/2 i .
Then the scheme is conservative:

δx i∈Z ρ n+1 i = δx i∈Z ρ n i = M. (iii) Estimates. Since ρ n+1/2 i = ρ n i ,
ρ n+1 i = ρ n+1/2 i + λ 2 (ρ n+1/2 i+1 -2ρ n+1/2 i + ρ n+1/2 i-1 ) - λ 2 (J n+1/2 i+1 -J n+1/2 i-1
), (3.22) where

J n+1/2 i = f n+1/2,+ i -f n+1/2,- i
. We introduce the quantity

γ n+1/2 i+1/2 =        J n+1/2 i+1 + J n+1/2 i ρ n+1/2 i+1 + ρ n+1/2 i , if ρ n+1/2 i+1 + ρ n+1/2 i = 0 0 otherwise
We clearly have that |γ

ρ n+1 i = ρ n i + λ 2 (ρ n i+1 -2ρ n i + ρ n i-1 ) - λ 2 γ n+1/2 i+1/2 (ρ n i+1 + ρ n i ) -γ n+1/2 i-1/2 (ρ n i + ρ n i-1 ) . (3.23)
As in the proof of Theorem 2.6, we introduce the quantity

M n i = δx k≤i ρ n k . By definition of M n i , we have δx(ρ n i+1 + ρ n i ) = M n i+1 -M n i-1
. Therefore, summing (3.23), we deduce

M n+1 i = M n i (1 -λ) + λ 2 M n i+1 (1 -γ n+1/2 i+1/2 ) + λ 2 M n i-1 (1 + γ n+1/2 i+1/2 ). (3.24)
Thus M n+1 i is a convex combination of M n i-1 , M n i and M n i+1 , provided the CFL condition (3.20) is satisfied. It is then obvious by an induction on n to deduce that for all i and n, 0 ≤ M n i ≤ M and that we have a BV -estimate on the sequence (M n i ). Using (2.10), we deduce, by using ρ n k = (M n k -M n k-1 )/δx, that

ν n i = k∈Z M n k (i-k)δx (i-1-k)δx w(z + δx) -w(z) δx dz.
Since the function w is bounded, we deduce that the integral in the right hand side is bounded. Therefore, from the bound on (M n i ) i,n , we deduce a BV -bound on the sequence (ν n i ) i,n . From (3.13), we have

- ∂ x S n i+1/2 -∂ x S n i-1/2 δx + ν n i = ρ n i . (3.25)
Summing this latter equation on i, we deduce,

∂ x S n i+1/2 = k≤i δxν n k -M n i . (3.26)
Using moreover (3.21), we have the bound for all i, n

|∂ x S n i+1/2 | ≤ M (1 + w 0 ).
(iv) Passing to the limit ε → 0. We will denote with a tilde all the limits of considered quantities when they exist.

From the L ∞ ∩ BV bound independant of ε on the sequences (M n i ) i,n and (ν n i ) i,n , we deduce that we can extract a subsequence that converges strongly in L 1 loc as ε → 0 to ( M n i ) i,n and ( ν n i ) i,n respectively. We deduce using moreover (3.26) that, up to a subsequence, (∂ x S n i+1/2 ) i,n converges strongly in L 1 loc as ε → 0. Moreover, since the sequence (γ

n+1/2 i+1/2 ) i,n is bounded in L ∞ independantly of ε, we can extract a subsequence converging in L ∞ -weak * as ε → 0 to ( γ n+1/2 i+1/2 ) i,n .
Taking the limit ε → 0 in (3.24), we deduce that the limit sequences ( M n i ) i,n and ( γ n+1/2 i+1/2 ) i,n satisfy the same relation (3.24). Defining

ρ n i = ( M n i -M n i-1
)/δx, this sequence satisfies equation (3.23). Moreover, we have the weak convergence in S M of (ρ n i ) i,n towards ( ρ n i ) as ε → 0. Letting ε → 0 the relation (3.25) is still available for the limit quantities.

From the strong convergence, we can pass to the limit in (3.14) to deduce the L ∞ -weak * convergence of (a n i ) i,n . Passing to the limit in (3.17), we deduce that the sequence (Π n,± i ) i,n converges to a limit still satisfying (3.17) with a tilde on all quantities. This equation can also be rewritten

Π n,± i = φ n,± i ρ n i = 1 2 (1 ± a n i ) ρ n i . Finally, from (3.18), we have J n+1/2 i = e -δt/ε J n i + (1 -e -δt/ε )(Π n,+ i -Π n,- i ). Since |J n i | ≤ M ,
we can take the weak limit as ε → 0 and deduce that

J n+1/2 i ⇀ Π n,+ i -Π n,- i = a n i ρ n i .
We conclude the proof by passing to the limit ε → 0 in (3.22).

Numerical simulations

We present in this Section some numerical examples to illustrate our results. In particular, we present two examples with applications in biology or plasma physics, where w = 0 or w = W .

The analytical results presented above are setted on the whole real line R. For the numerical simulations, we have to restrict our domain of computation to a finite domain of R. Then boundary conditions should be fixed. This point is really challenging due to the convolution products which requires to know the solution on the whole real line. Since we are focusing on aggregation equation, it is known (see e.g. [START_REF] James | Chemotaxis: from kinetic equations to aggregate dynamics[END_REF][START_REF] Carrillo | Global-intime weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]) that initial data with compact support are compactly supported and collapse in finite time into a single Dirac delta. Therefore analytical solutions never reach the boundaries. Then, since we expect the density to be 0 outside of the domain, we impose homogeneous Dirichlet boundary conditions on the density (and on the distribution function for the kinetic equation) and homogeneous Neumann boundary conditions for the potential S.

Simulation of an aggregation equation

In this subsection we consider the case W = -1 2 |x| and a = id. Then the equation writes

∂ t ρ + ∂ x ((W ′ * ρ)ρ) = 0.
This equation appears in several applications in biology or physics, see for instance [START_REF] Nieto | High field limit for Vlasov-Poisson-Fokker-Planck equations[END_REF] where this system is the high field limit of a Vlasov-Poisson-Fokker-Planck system, the quantity S = W * ρ being the solution of the Poisson equation. In biology, it can be seen as a Patlack-Keller-Segel model without diffusion, the quantity S being the chemoattractant concentration.

Numerical scheme (2.4)-(2.5) is implemented. We notice that in the case a(x) = x, (2.8) rewrites a n i+1/2 = 1 2 ∂ x S n i+1 + ∂ x S n i . Numerical results are display in Figure 1 for two different initial data. In Figure 1 left, we take ρ ini (x) = e -10x 2 . We observe that the initial bump stiffens and the solutions blows up in finite time to form one stationary single Dirac. In Figure 1 right we take ρ ini (x) = e -10(x-1.25) 2 + 0.8e -20x 2 + e -10(x+1) 2 . As for the previous initial data, the initial bumps blow up and collapse in one single Dirac mass in finite time. In this case W ′′ = -δ 0 . Then, setting S = W * ρ, we have that -∂ xx S = ρ and -∂ x S is an antiderivative of ρ. Then, integrating the aggregation equation, we can rewrite it as

∂ t ∂ x S + 1 2 ∂ x ∂ x S 2 = 0.
We recognize the Burgers equation for ∂ x S. Moreover, in this particular case where ν n i = 0, we can deduce from (2.4)-(2.9) a scheme on (∂ x S n i ) i,n . First, (2.5) rewrites

J n i+1/2 = - 1 δx A(∂ x S n i+1 ) -A(∂ x S n i ) .
Moreover, denoting ρ n i+1/2 = 1 2 (ρ n i + ρ n i+1 ), we deduce from (2.9) and (2.7)

ρ n i+1/2 = - 1 δx ∂ x S n i+1 -∂ x S n i .
We deduce k≤i ρ n k+1/2 = -

1 δx ∂ x S n i+1 .
Equation (2.4) implies straightforwardly

ρ n+1 i+1/2 = ρ n i+1/2 (1 -λc) + λ 2 c(ρ n i-1/2 + ρ n i+3/2 ) + λ 2 (J n i-1/2 -J n i+3/2 ).
Summing this latter equation, we deduce

∂ x S n+1 i = ∂ x S n i (1 -λc) + λ 2 c(∂ x S n i-1 + ∂ x S n i+1 ) - λ 2 A(∂ x S n i+1 ) -A(∂ x S n i-1
) .

In the case at hand where a = id, we have A(x) = x 2 /2, we recognize the well-known Lax-Friedrichs discretization for the Burgers equation. Here we have c = M where M is the total mass of the system. Then the numerical results of Figure 1 were expected; we recover the convergence in finite time in a single Dirac mass as established for instance in [START_REF] Carrillo | Global-intime weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]Section 4].

A kinetic model for chemotaxis

Let us consider the so-called Othmer-Dunbar-Alt model, describing the motion of cells by chemotaxis, in one dimension. This model has been used since the 80's when it has been observed that the motion of bacteria is due to the alternance of straigth swim in a given direction, called run phase, with cells reorientation to choose a new direction, called tumble phase. This system governs the dynamics of the distribution function f ε . In the hyperbolic scaling, it reads:

∂ t f ε + v∂ x f ε = 1 ε V T [v ′ → v]f ε (v ′ ) -T [v → v ′ ]f ε (v) dv ′ .
In this equation T [v ′ → v] is the turning rate, corresponding to the probability of cells to change their velocities from v ′ to v during a tumble phase. In [START_REF] Dolak | Kinetic models for chemotaxis: Hydrodynamic limits and spatio-temporal mechanisms[END_REF] the following expression has been taken for the turning rate

T [v ′ → v] = φ(v ′ ∂ x S ε ).
In this equation, the quantity S ε corresponds to the chemoattractant concentration which solves the elliptic equation

-∂ xx S ε + S ε = ρ ε , (4.1)
where ρ ε = V f ε (v) dv is the density of cells. This latter equation can be rewritten S ε = W * ρ ε for W = 1 2 e -|x| ; therefore we have W = w in (1.2). The velocity v is assumed to have a constant modulus which is fixed to 1 by normalization. Then, in one dimension, the kinetic equation reduces to:

∂ t f ± ε ± ∂ x f ± ε = 1 ε (φ(∓∂ x S ε )f ± ε -φ(±∂ x S ε )f ± ε ). (4.2)
The density of cells is defined by

ρ ε := f + ε + f - ε .
We assume that the turning rate satisfies (3.1). Then, we can rewrite (4.2) as (3.3). Applying Theorem 3.2, the limiting model when ε → 0 is given by

∂ t ρ + ∂ x J = 0, J = -∂ x A(∂ x S) + a(∂ x S)S. (4.3)
This equation has been studied in [START_REF] James | Chemotaxis: from kinetic equations to aggregate dynamics[END_REF].

Remark 4.1 As in the first example of this Section, we can recover an equation for the potential S = W * ρ, which turns out here to be nonlocal. Indeed, taking the convolution with W = 1 2 e -|x| of (4.3), we obtain

∂ t S + A(∂ x S) -W * A(∂ x S) + ∂ x W * (a(∂ x S)S) = 0.
Then by recombining (4.1) and (4.3), this latter equation can rewrite

∂ t S -∂ txx S + ∂ x -∂ x A(∂ x S) + a(∂ x S)S = 0.
It bears some resemblance with the well-known Camassa-Holm equation [START_REF] Camassa | An integrable shallow water equation with peaked solitons[END_REF], and exhibits the same peakon-like solutions. However the underlying dynamics is completely different and in the present case peakons collapse. Notice also that there are no anti-peakons because of the positivity of ρ.

The computational domain is assumed to be [-2.5, 2.5] and the velocity v is normalized to 1. We consider the function a(x) = 2/π Arctan(10x), which clearly satisfies (1.3), and recall that W = 1 2 e -|x| .

Macroscopic model

First we consider (4.3) and discretize the system thanks to (2.4)-(2.9), with ν n i replaced by S n i in (2.9). Figure 2 displays the numerical results for the following two initial data:

ρ ini (x) = e -10(x-0.7) 2 + e -10(x+0.7) 2 , (4.4)

ρ ini (x) = e -10(x-1.25) 2 + 0.8e -20x 2 + e -10(x+1) 2 . (4.5)
As expected, we have a fast blow up of the regular solution and a finite-time collapse in a single Dirac mass, at a rate much higher than in the case W = -1 2 |x|, compared with Figure 1-right. The behaviour of such Dirac solutions can be recovered by studying solutions in the form ρ(t, x) = n i=1 m i δ(x -x i (t)). Then we have S(t, x) = W * ρ(t, x) = 1 2 n i=1 m i e -|x-x i (t)| . After straightforward computations, we deduce from the expression in (4.3) that

J = - n i=1 [A(∂ x S)] x i δ(x -x i (t)),
where the notation [f ] x i denotes the jump of the function f at the point x i . In particular, we have that ρ satisfies system (4.3) provided,

m i x ′ i (t) = -[A(∂ x S)] x i .
Moreover, the function a being increasing and odd, the function A is strictly convex and can be chosen even. Then, equilibrium states satisfy

-[A(∂ x S)] x i = A 1 2 m i + j =i m j e -|x j -x i | -A 1 2 -m i + j =i m j e -|x j -x i | = 0.
This equality is true only if j =i m j e -|x j -x i | = 0, which implies n = 1. Therefore stationary states are given by a single stationary Dirac mass. Convergence towards this equilibrium is proved in [START_REF] James | Chemotaxis: from kinetic equations to aggregate dynamics[END_REF][START_REF] Carrillo | Global-intime weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF].

Kinetic framework

We turn now to the kinetic framework and implement the scheme described in Section 3. In Figure 3 and 4 we display the dynamics of the cell density ρ for the two regular bumps initial data (4.4) and an initial distribution function given by f ini,± (x) = 1 2 ρ ini (x). We plot in the left part of the figures the numerical results corresponding to the macroscopic model (4.3), whereas the right part corresponds to numerical solution of (4.2) with ε = 0.1. The macroscopic behaviour is the same as in Fig. 2-left, while in the kinetic case (right), the solution does not blow up, as it is expected. However, the behaviour is similar to the one for the macroscopic model: we observe the formation of two interacting aggregates that attract one another to collapse in a single aggregate in finite time. This behaviour can be interpreted formally, since in the kinetic case, stationary states for (4.2) are given by

±∂ x f ± ε = 1 ε (φ(∓∂ x S ε )ρ ε -f ± ε ).
Summing these two equations, we easily deduce that ∂

x (f + ε -f - ε ) = ∂ x J ε = 0. Since we expect J ε = 0 at infinity, we can take f + ε = f - ε = f ε .
Then, using the expression of φ in (3.1), the previous equation rewrites ∂

x f ε = 1 ε a(∂ x S ε )f ε .
Formally, if we consider that the function a is an approximation of the sign function, then the latter equation is a linear ODE which can be solved easily and implies that f is given by a sum of exponential functions with tail ± 1 ε . This behaviour corresponds to what is observed in Figure 4 right.

We illustrate now the asymptotic preserving property of the scheme: in Figure 5 are displayed numerical results for different values of ε and for the three regular bumps initial data (4.5). As above, we notice the formation of 3 aggregates that merge into one single aggregate. When ε → 0, we observe that the numerical solutions converges to the one computed in the macroscopic case, which is an illustration of the result of Theorem 3.3, compare with Figure 2-right.

Specific issues

We conclude this section with some important remarks concerning the choice of the schemes above. First, we emphasize the importance of the choice of the discretized macroscopic velocity. For instance, in the aggregation equation (1.1), if instead of defining the discretization (2.8) we take a n i+1/2 = a((S n i+1 -S n i )/δx), we obtain Figure 6, to be compared with Figure 3. Concerning the kinetic model, we display in Figure 7 the results obtained when the discretization of a n i in (3.14) is replaced by a n i = a(∂ x S n i ), with ∂ x S n i = (S n i+1 -S n i-1 )/2δx. We notice in Figure 7 (left) that for ε = 10 -3 the behaviour of the density remains comparable with the macroscopic model (see Figure 2). When ε goes to zero, namely here ε = 10 -5 , we observe the same kind of result as for the macroscopic case, compare Figures 7 (right) and 6. This emphasizes that in some sense the asymptotic preserving property has to be complemented with a careful definition of the nonconservative product when Dirac masses appear in the solution. Notice also that for weak relaxation this definition is not so crucial.

Next we recall that the time-splitting strategy we have chosen here for its simplicity, which allows us to prove rigorous convergence results, is known to have some drawbacks concerning the macroscopic fluxes. Indeed these are expected to be flat for stationary solutions: this is the so-called well-balanced property. This is not the case here, as evidenced in Figure 8, where are displayed the numerical fluxes (J ε = f + ε -f - ε ) at the final time corresponding to the numerical simulations presented in Figure 5. This illustrates the fact that the flux converges to zero merely in a weak sense. To obtain better results, well-balanced schemes have to be considered, as in e.g. [START_REF] Gosse | Computing qualitatively correct approximations of balance laws[END_REF][START_REF] Gosse | A well-balanced scheme for kinetic models of chemotaxis derived from onedimensional local forward-backward problems[END_REF][START_REF] Gosse | Asymptotic-preserving and well-balanced schemes for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF]. However these results are obtained mostly for smooth solutions, hence likely can be applied for weak relaxation (that is "large" values of ǫ), but their behaviour when ǫ goes to zero is far from clear. The problem is actually twofold: first developing an asymptotic preserving scheme with the well-balanced property for smooth solutions, next making such a scheme compatible with the nonconservative product when measure-valued solutions appear. [2] L. Ambrosio, N. Gigli, G. Savaré, Gradient flows in metric space of probability measures, Lectures in Mathematics, Birkäuser, 2005

  | ≤ 1. Then, equation (3.22) rewrites

Figure 1 :

 1 Figure 1: Dynamics of the density ρ with two different intial data in the case W = -1 2 |x|

Figure 2 :

 2 Figure 2: Macroscopic model (4.3): cell density ρ for initial data (4.4) (left) and (4.5) (right)

Figure 3 :

 3 Figure 3: Dynamics of the cell density ρ for an initial data given by a two regular bumps: comparison between the macroscopic model (left) and the kinetic model for ε = 0.1 (right).

Figure 4 :

 4 Figure 4: Dynamics of the cell density ρ for an initial data given by a two regular bumps: comparison between the macroscopic model (left) and the kinetic model for ε = 0.1 (right).

Figure 5 :

 5 Figure 5: Dynamics of the cell density for an initial data given by a sum of 3 regular bumps: simulation of the kinetic model with ε = 0.1 (left) and ε = 10 -3 (right).

Figure 6 :

 6 Figure 6: Wrong velocity discretization for (1.1)

Figure 7 :

 7 Figure 7: Wrong velocity discretization for kinetic model: ε = 10 -3 (left) -ε = 10 -5 (right)

Figure 8 :

 8 Figure 8: Flux at the final time T = 6 corresponding to the simulations in the Figure 5: kinetic model with ε = 0.1 (left) and ε = 10 -3 (right).

  we have S Due to the mass conservation, we still have the bound in (2.15), i.e. for all i and n,

	n+1/2 i i . δx = S n i and ν n+1/2 i = ν n |ν n i | ≤ M w 0 .	(3.21)
	j≤i
	Adding the equations (3.19), we deduce,

Acknowledgments. The authors would like to acknowledge warmly the anonymous referee for the careful reading and for the constructive remarks that helped to improve this paper.

The second author is supported by the French "ANR programme blanc" project Kibord: ANR-13-BS01-0004.


Tightly Coupling GPS with Lane Markings for Autonomous Vehicle Navigation

Zui Tao 1,2 , Philippe Bonnifait 1,2 Abstract-Tightly coupling GPS pseudorange and Doppler measurements with other sensors is a way to increase accuracy and integrity of the positioning information particularly when it is computed autonomously. Highly accurate digital maps are also more and more key components for autonomous vehicle navigation and can enhance the localization system. In this paper, a video camera is used to get relative information with respect to lane markings and dead-reckoning sensors are also integrated to provide positioning information with high availability. A reduced order state space modeling of the observation problem is proposed in order to get an efficient realtime computing system. Thanks to an adequate error model of the GPS pseudoranges and a measured input state space formalization, a Kalman filter with correlated noises is developed. Experimental results show that this tightly coupled approach clearly improves the performance in terms of accuracy and integrity compared to a loosely coupled method that uses GPS fixes computed by an external receiver.

I. INTRODUCTION

Recently, impressive demonstrations of self-driving cars using production-based sensors and enhanced maps have been made on both rural and urban routes [START_REF] Ziegler | Making bertha drive an autonomous journey on a historic route[END_REF][2] [START_REF] Schreiber | Laneloc: Lane marking based localization using highly accurate maps[END_REF]. In most cases, an on-road assumption and the use of an informative digital map make it possible to realize low-cost autonomous driving. In the research community of Intelligent Transportation Systems, generation of detailed digital maps has done significant progress [START_REF] Betaille | Creating enhanced maps for lanelevel vehicle navigation[END_REF] [START_REF] Brummer | Evaluation of a mapping strategy based on smooth arc splines for different road types[END_REF] to the point that mapaided perception [START_REF] Kurdej | Controlling remanence in evidential grids using geodata for dynamic scene perception[END_REF] and localization are now intensively considered for autonomous navigation.

Global Navigation Satellite Systems are widely used in vehicle navigation applications. When integrating GPS with other sources of information, there are two main kinds of approaches: a loosely coupled approach which uses directly the GPS fix and a tightly coupled one which uses the GPS observables like pseudoranges, Dopplers and carrier phase measures. Loosely coupled methods have to estimate biases on the GPS fixes to improve the fusion process. The bias error is typically due to atmospheric effects and can be hardly observed by a stand-alone GPS loosely coupled with deadreckoning (DR) sensors. In [START_REF] Jo | Gps-bias correction for precise localization of autonomous vehicles[END_REF], the GPS bias is estimated using the perception of lanes and crosswalks stored in the map data. In [START_REF] Tao | Mapping and localization using gps, lane markings and proprioceptive sensors[END_REF], different GPS bias models are proposed and compared. However, loosely coupled methods depend on the GPS receiver and are behind its software processing which is particularly difficult to model. In this paper, we thus propose to estimate the bias on every pseudorange.

The benefits of using a map in tight coupling has been proved for a long time. In [START_REF] Taylor | Road reduction filtering for gps-gis navigation[END_REF], a Road Reduction filter is

The authors are with 1 Université de Technologie de Compiègne (UTC), 2 CNRS Heudiasyc UMR 7253. developed to calculate pseudorange corrections derived from the digital road network. Since no perception information is required, this macro-scale approach introduces a map pseudo-measurement by projecting the GPS fix onto the road center-line. A navigation system incorporating this approach can achieve an horizontal position accuracy of 10m (1σ) typically [START_REF] Groves | Principles of GNSS, inertial, and multi-sensor integrated navigation systems[END_REF]. An alternative is to use the map as a heading sensor by considering that the map precision is often better than its accuracy [START_REF] Fouque | On the use of 2d navigable maps for enhancing ground vehicle localization[END_REF]. In order to benefit from a highly accurate map as described in [START_REF] Betaille | Creating enhanced maps for lanelevel vehicle navigation[END_REF] and to remove non-zero mean errors introduced by map matching on poly-lines, a vehicle needs to measure its relative position w.r.t. the lane in which it is traveling by integrating a perception information such as lane marking detections. With this objective in mind, this paper proposes a tightly coupled approach that integrates a camera observations and a lane marking map similar to the one used in [START_REF] Tao | Lane marking aided vehicle localization[END_REF] by using an extended Kalman filter with correlated noises (EKF-CN). Recent results on tightly coupling can be found in [START_REF] Rabaoui | Dirichlet process mixtures for density estimation in dynamic nonlinear modeling: Application to gps positioning in urban canyons[END_REF] [START_REF]Pseudoranges error correction in partial gps outages for a nonlinear tightly coupled integrated system[END_REF].

Loosely coupled methods as the one proposed in [START_REF] Tao | Lane marking aided vehicle localization[END_REF] can improve the lateral localization accuracy up to a decimeterlevel. The longitudinal error is also controlled thanks to the heading variation of the trajectory. In the tightly coupling case, the longitudinal uncertainty can be further corrected by the pseudoranges whose line-of-sight directions are along the driving direction as shown in [START_REF] Vu | Realtime computer vision/dgps-aided inertial navigation system for lanelevel vehicle navigation[END_REF].

Important performance metrics are not limited to accuracy. Integrity is also crucial for autonomous vehicles. The integrity risk is defined as probability that the position error exceeds predefined alert limits. In this paper, the integrity risk of the estimated position is set to 10 -2 (1%). By ignoring the faults of the sensors and by assuming that the noises follow Gaussian distributions, the corresponding bound of the estimated position is ±2.58σ. Therefore, an integrity failure happens when the estimated horizontal positioning error (HPE) is beyond this bound centered on the estimate. As discussed in [START_REF] Drevelle | Localization confidence domains via set inversion on short-term trajectory[END_REF], Bayesian state filtering is often overconfident, and so we are interested in studying here this issue and comparing the integrity performance of both loosely and tightly coupled methods.

The remainder is organized as follows. Section II introduces the system modeling and reminds the lane marking aided vehicle localization. Section III introduces the GPS measurement models and a reduced-order Kalman filtering is proposed for efficient embedded computation. Results of outdoor experiments are presented and analyzed in section IV. Conclusions on the proposed tightly coupled method compared to a loosely coupled one are given in section V.

II. LANE MARKING LOCALIZATION

If the vehicle approximately knows where it is located in the map, it is possible to implement a lane marking localization method that corrects DR estimates based for instance on wheel odometry. We explain in this section the key concepts to implement such an approach.

A. Frames

A local navigation frame R O tangent to the Earth is defined to have its x axis pointing to the East, y to the North and z oriented upwards w.r.t. the WGS84 ellipsoid (Fig. 1a). By choosing the origin O close to the navigation area, the vehicle trajectory is locally planar. Two more frames have to be defined. In Fig. 1b, R M denotes the mobile vehicle body frame with its longitudinal axis x M pointing forward and its lateral axis y M pointing left hand. The camera frame R C is located at the front of the vehicle even if the camera is located behind the windscreen with a position offset (Cx, Cy) in R M . P x is the translation from point M up to the bumper. In the following, the coordinates are expressed by default in R O .

B. Dead-reckoning

The linear velocity of each rear wheel is measured by the ABS speed sensors. The linear velocity is calculated by

v m = (v m rl + v m rr ) /2
, where v m rl and v m rr denote the measured linear velocity of the left and right rear wheels respectively.

The angular velocity of the vehicle (ω m ) is measured by the ESP yaw rate gyro. An estimate of the pose of the vehicle is provided by integrating these measurements from an initial known pose.

Please note that the measurements are affected by noise (v m = v + γ v ; ω m = ω + γ ω ) and so the dead-reckoned estimate drifts with respect to the traveled distance.

C. Lane marking map

The lane marking map used in this paper mainly consists of two-lane roadways with dashed lane markings in the center of the road and solid lane markings on the both sides of the road. Lane markings are expressed by polylines. The map is expressed in R O and reaches a centimeter-level accuracy.

D. Camera measurement

Let L denote the lane marking detection located at ordinate C 0 in R C (see Fig. 1b). Its coordinates in R O are given by (x, y being the position and θ the heading angle):

x L y L = P x • cosθ + C 0 • sinθ + x P x • sinθ -C 0 • cosθ + y (1)
In Fig. 1b, [AB] represents the detected lane marking segment. The coordinates of point A and B are (x A , y A ) and 2), we have:

(x B , y B ) in R O . Let us define a vector V = (x AB , y AB ) T with x AB = x B -x A and y AB = y B -y A . Point L on segment [AB] meets: x L = x A + λ • x AB y L = y A + λ • y AB (2) E X Y Z O x y z λ Ф (a) B A C C 0 O x y M C 1 L y O x O x C y C y M x M P x C y C x (b)
P x • cosθ + C 0 • sinθ + x = x A + λ • x AB P x • sinθ -C 0 • cosθ + y = y A + λ • y AB (3)
Because x AB •cosθ+y AB •sinθ = 0 (the detected segment is not perpendicular to the vehicle), the camera observation model is derived by Eq. ( 3) as:

C 0 = (Px • sinθ + y -y A ) x AB -(Px • cosθ + x -x A ) y AB x AB • cosθ + y AB • sinθ + βc (4)
where β c is the measurement noise.

E. Map-matching

Map-matching consists to determine which is the lane marking segment [AB] that the camera has detected. In this paper, the marking type (e.g. dashed or solid) given by the video camera is used to map-match the detected lane marking segment. In a first stage, a set S of candidate segments is selected based on the following conditions:

• The lane marking type is consistent with the one estimated by the camera; • The orientation of the segment is close to the heading of the vehicle; • The distance dist between point L and the candidate segment is less than the width of the road. In a second stage, the segment s which has the minimal dist is chosen as the map-matching result:

Map_matched_segment = arg min s∈{S} {dist} (5)

III. TIGHTLY COUPLING GPS

In this section, we present a strategy for tightly coupling the GPS raw measurements with the lane marking localization.

A. GPS measurements

The GPS raw measurements considered here are L1 pseudoranges and Dopplers on the visible satellites.

The position of the GPS antenna with respect to the body frame is taken into account.

   xa = x + cosθ • tx -sinθ • ty ya = y + sinθ • tx + cosθ • ty za = z + tz (6)
Where [x, y, z] T is the 3D position of the vehicle, [tx, ty, tz] T is the position of the receiver antenna in R M and Xa = [xa, ya, za] T is the position of the receiver antenna in R O .

For a given satellite i, its position vector X i s = x i , y i , z i T is reconstructed from the received navigation message (GPS data). The corresponding pseudorange is:

ρ i = R i + c • dt u -c • dt s + δ T GD + δ rel + δ iono + δ tropo (7)
with R i being the geometrical distance between the observed satellite and the receiver antenna R i = Xa -X i s , c • dt u the range equivalent of the receiver clock offset, c • dt s the range equivalent of the satellite clock offset, δ T GD the timing Group Delay, δ rel the range equivalent of the relativity effect on the satellite clock, δ iono the ionospheric delay and δ tropo the tropospheric delay.

For every satellite, the GPS receiver gets also, in its ephemeris data, some information to compute estimates of dt s , δ T GD and δ rel . It can also estimate δ iono and δ tropo by using classical atmosphere models implemented in open source software like [START_REF] Harris | The gpstk: an open source gps toolkit[END_REF]. By correcting the pseudorange with these estimated parameters, equation 7 becomes:

ρ i = R i + c • dt u + ε i pr + β i (8)
where ε pr represents the residual (non-white) errors of the pseudorange. β i is the measurement noise. c is the speed of light.

The Doppler shift is caused by the relative motion between the satellite and the receiver antenna. By defining the lineof-sight vector u i los of satellite i as,

u i los = X a -X i s /R i (9)
the Doppler shift is measured and linked to the GPS antenna trough the following equation [START_REF] Fouque | Matching raw gps measurements on a navigable map without computing a global position[END_REF]:

ρi = V r -V i s • u i los + c δ t u + β i d (10)
Where Vr = [ ẋ, ẏ, ż] T is the velocity of the receiver and

V i s = ẋi , ẏi , żi
T is the one of satellite i. • denotes the dot product. β i d is the measurement noise. For a numerically good computation, the filter estimates ranges equivalent values of the receiver clock parameters:

d = c • dt u ḋ = c •dt u

B. System modeling

Let us look for a planar model of the system in order to have an efficient state observer in terms of computations with a limited number of sensors. The localization problem can be reduced to 2D by considering z a as constant during the current navigation stage. The main sources of range error are (see Fig. 2):

• Pseudoranges residuals errors ε pr ;

• Inaccurate satellite position estimates due to the use of real-time navigation messages ε Sat ; • The fact that the evolution area is not exactly planar ε z . We propose to estimate a range-error parameter ε that results from the combination of ε Sat , ε pr and ε z . A first order auto-regressive shaping filter is used by (the subscript k being the time stamp):

ε i k = a • ε i k-1 + α ε k-1 (11)
Where α ε is the driving noise.

Let us consider a front-wheel drive vehicle with the assumption of a slip-free motion of the rear wheels. The speed vector is therefore collinear to x M .

With n satellites in view, the state vector X becomes:

X = x, y, θ, b, d, ḋ, ε 1 , • • • , ε n T ( 12
)
Where b is the bias of the gyro. With a 2D unicycle kinematic model, the evolution model is given by:

X k = f X k-1 , U m k , α k ⇔                              x k = x k-1 + T • v m k • cosθ k-1 y k = y k-1 + T • v m k • sinθ k-1 θ k = θ k-1 + T • (w m k -b k-1 ) b k = b k-1 + α b k d k = d k-1 + T • ḋk-1 + α d k ḋk = ḋk-1 + α ḋ k ε 1 k = a • ε 1 k-1 + α ε1 k . . . ε n k = a • ε n k-1 + α εn k (13)
Where U m = [v m , ω m ] T denotes the measured input vector and T the sampling period. A random constant model has been considered to handle the bias b.

The wheel speed sensor measurement noise υ v and the gyro measurement noise υ ω are supposed to be zeromean independent white noises. N denotes their covariance matrix. The covariance of the model noise 13) is denoted by Q.

α k = α b k α d k α ḋ k α ε1 k . . . α εn k T of equation (

C. GPS measurement validation

At time instant k, when the GPS measurements are available, a validation step is performed on the measurements of every satellite to avoid the use of badly tracked satellites and to reject multipath.

For a Doppler measurement, the validation process is:

• C/N 0 gating: Check that C/N 0 is high enough (e.g. 38 dBHz).

• Elevation mask: The elevation angle of satellite i is calculated using X i s and the current estimate X k|k-1 . The elevation mask angle is usually set to 15 degrees.

• Innovation gating: Check the Normalized Innovation Squared [START_REF] Fouque | Enhancement of global vehicle localization using navigable road maps and dead-reckoning[END_REF]. A Doppler is in general more reliable than a pseudorange. Therefore, a pseudorange can be used only when its corresponding Doppler measurement is valid. For every pseudorange, its innovation η i is computed and used to validate or reject the measurement.

By

denoting A k = ∂f ( X k-1|k-1 ,U m k ) ∂X , B k = ∂f ( Xk-1|k-1 ,U m k ) ∂U m
, the tightly coupling solver is described by Algorithm 1 in which the Dopplers are used at first. Additionally a Mahalanobis distance threshold is set to reject the mismatch before using C 0 to update the state vector. The process is time-triggered with the CAN bus data which has the highest rate. The latencies of the GPS and camera measurements are, in the current implementation, neglected. Since the camera provides essentially lateral corrections, neglecting its latency has almost no impact on the accuracy.

Algorithm 1 An iteration stage of the method

In out: X, P //estimated state and covariance matrix X, P =Update( X,P ,Good_Doppler) 13:

1: U m = [v m , ω m ] T =Get(DR measurements) 2: X = f X, U m 3: P = A k • P • A T k + B k • N • B T k + Q 4: if New GPS data is available then 5: ρ 1,...,n , ρ1,•••n =Get(GPS measurements
for j = 1, . . . , n do 14:

if ( ρj is valid )&(η i < Threshold) then 15:
Add(ρ j ) to the Good_Pr list [C 0,k , laneT ype] =Get( camera measurements ) 22:

[AB]=Map_match X, C 0,k , laneT ype, map end if 26: end if One can note that the filter uses a measured linear velocity v m in the Doppler observation model (Eq. 10). Therefore, the estimation process has correlated noises. Please refer to [START_REF] Bonnifait | A system for 3d localization of civil-engineering machines[END_REF] for more details about the EKF-CN equations.

Model noises variances

Measurement noises variances V ar α b = 5 × 10 -10 V ar (βc) = 0.16

V ar α d = 1 × 10 -3 V ar (β) = S • 10 -C/N i 0 10 V ar α ḋ = 1 × 10 -4 V ar (β d ) = 0.05 V ar (α ε ) = 1 × 10 -4 V ar (γ v ) = 1 × 10 -4
V ar (γ ω ) = 2.5 × 10 -3

IV. EXPERIMENTAL RESULTS

Outdoor experiments have been carried out near Paris France in May 2013. Three tests were performed on the same road, in urban conditions, with an automotive experimental vehicle. Test 1 and Test 2 were conducted in the morning and test 3 in the afternoon. During the test, the vehicle passed in a strong urban canyon of 300 meters with satellites outages and multipaths. The traveling distance of every test was about 2km with a typical speed of 30 km/h in the straight lines. The experimental vehicle was equipped with an IMU Oxford RT3000 which provided ground truth data at 100Hz rate. A CAN-bus gateway was used to access to the wheel speed sensors and to the yaw rate gyro. The measured input

[v m , w m ]
T from the CAN bus were available at 100Hz. A Mobileye camera was used to detect the lane markings at 10Hz. A low-cost U-blox 6T GPS receiver with a patch antenna providing pseudoranges and Dopplers measurements at 5Hz was used.

The measurement noise β i in Eq. ( 8) is supposed to be white and zero mean. As it is not stationary as indicated in [START_REF] Wieser | Improved positioning accuracy with high-sensitivity gnss receivers and snr aided integrity monitoring of pseudo-range observations[END_REF], its variance σ 2 i is modeled by a Wieser's model in which the measured carrier-to-noise density ratio C/N 0 is used as a quality indicator:

σ 2 i = S • 10 -C/N i 0 10 (14)
Where S = 60000m 2 Hz. Table I specifies every model and measurement noise. 10 satellites were in view during the different tests. The GPS satellite visibility was sometimes very constrained due to the urban canyon and due to buildings near the test area.

Localization performance can be studied in terms of longitudinal and lateral performances. Fig. 3 shows the lateral and longitudinal positioning errors of test 2 for loosely and tightly coupling. After t = 350s, the vehicle is back to the parking lot where lane markings are not stored in the map and so the localization system is not map-aided in this area.

The bias on every pseudorange is initialized to zero. For each subplot of Fig. 4, the abscissa expresses time on seconds, the ordinate gives every estimated bias ε i on meters. Satellites 4 and 13 are under the elevation mask (15 • ) and so their measurements are not used in the localization solver. Moreover, the measurements of satellite 26 became available only after 100s and were lost after 200s. The estimated biases are quite smooth, in the order of few meters and stay bounded during all the trial.

Table II gives performance metrics for both the loosely coupled solver (LC) proposed in [START_REF] Tao | Lane marking aided vehicle localization[END_REF] and the tightly coupled solver (TC) proposed in this paper. Fig. 5 shows the cumulative distribution of horizontal positioning errors of the three tests by four methods: standalone GPS, loosely coupling DR/GPS/Map, tightly coupling DR/GPS and tightly coupling DR/GPS/Map. It is noticeable that tightly coupling GPS with DR improves the accuracy compared with stand-alone GPS. The accuracy of the tightly coupled method is further improved by integrating the lane marking map. The 95th percentiles of the four methods are respectively 0.88m, 1.54m, 2.12m and 4.64m. 96.8% of the tightly coupling result reaches sub-meter accuracy.

The loosely coupled method relies on the GPS fixes computed by stand-alone GPS receiver. When the GPS receiver suffers from big errors during a long time, the loosely coupled solver is not able to estimate the biases efficiently. Although the 95th percentile is 1.54m, the biggest error reaches 5.29m.

Fig. 6a shows the global integrity performance of the loosely and tightly coupled methods. The gray area corresponds to the overconfident occasions of the filters. The global failure rate of the loosely coupled method is about 22% (Table II). The tightly coupled method reduces the integrity failure rate to 7.6% which is the same order as the predefined integrity risk of 1% even if it is a bit higher. The tightly coupled method is therefore more reliable. Now let us evaluate the precision of the results by considering the distribution of the estimated 2.58σ bound. Fig. 7 shows the tightly coupled method provides confidence domains that are smaller than 1.26m during 95% of the time which is a good improvement compared to the 2.06m of the loosely coupled one. As confidence is in practice compared to a threshold to indicate "use" or "don't use" to the client application, it is important, in terms of availability of the positioning information, to provide as small as possible confidence zones. To resume this integrity analysis, one can notice the confidence domains computed by the tightly coupled approach are smaller than the ones of the loosely coupled approach and they are more reliable in the sense that the ground truth is more often included in the confidence zone.

V. CONCLUSION

In this paper, we have proposed and studied a method to merge raw GPS measurements and lane marking measurements. In order to reduce the dimension of the state vector, the linear and angular velocity are used as measured input. The vehicle is also assumed to travel on a flat surface thanks to the use of a navigation frame close to the navigation area. An error model has been proposed to estimate the errors on GPS pseudoranges. A tightly coupled EKF-CN solver has been developed to reduce the calculation and improve the real-time performance. For testing the approach, DR sensors from the ABS and ESP systems, a L1-GPS receiver and the lane marking measurements from the lane departing warning system of the experimental vehicle have been used. Outdoor experiments have been carried out to compare the loosely coupled and tightly coupled methods. The tightly coupled method gives a better accuracy and a better reliability. This system provides a very cost-effective solution for the localization task of autonomous vehicles.

Fig. 1 :

 1 Fig. 1: Frames with λ ∈ [0, 1]. Plugging Eq. (1) into Eq. (2), we have:

Fig. 2 :

 2 Fig. 2: Range-error sources

  )

  for j = 1, . . . , n do 8:if ( ρj is valid) then // Please refer to Section III

  AB] is consistent with the vehicle state) then 24:X, P = Update X, P, C 0,k , [AB] 25:

Fig. 3 :Fig. 4 :

 34 Fig. 3: Longitudinal/lateral positioning errors. TC: tightly coupling; LC: loosely coupling

Fig. 5 :

 5 Fig. 5: Cumulative distribution function (CDF) of the HPE for the different methods

Fig. 6 :Fig. 7 :

 67 Fig. 6: Integrity plot

TABLE I :

 I Tuned EKF-CN parameters (International System Units)

TABLE II :

 II Horizontal positioning error statistics

	1 95% 96.8%
	87.3%
		0.8			Stand-alone GPS
					Loosely coupling DR/GPS/Map
					Tightly coupling DR/GPS
					Tightly coupling DR/GPS/Map
		0.6
	53.7%
	F(HPE)	0.4
	24.9%
		0.2
		0 0	1 0.88	2 1.54	3 Horizontal positioning error: HPE (m) 4 5 4.64 2.12	6	7	8

Acknowledgments

The authors would like to thank Javier-Ibañez Guzman, Nabil Chaari and Stéphane Bonnet for their support in the experiments.