G Costeseque
email: costeseg@cermics.enpc.fr

J P Lebacque

A variational formulation for higher order macroscopic traffic flow models: numerical investigation

Keywords: Traffic flow, Hamilton-Jacobi equation, Lax-Hopf algorithm, Lagrangian data

est destinée au dépôt

Introduction 1.General background

In order to get a realistic estimation of the real-time traffic states on networks, traffic operators and managers need macroscopic traffic flow models. These models must be simple, robust, allowing to get solutions at a low computational cost. The main macroscopic models are based on conservation laws or hyperbolic systems (see [START_REF] Hoogendoorn | State-of-the-art of vehicular traffic flow modeling[END_REF][START_REF] Kuhne | Traffic Flow Theory: A State of the Art Report -Revised Monograph on Traffic Flow Theory[END_REF] for traffic aspects and [START_REF] Garavello | Traffic flow on networks[END_REF] for mathematical aspects). The seminal LWR model (for Lighthill-Whitham and Richards) was proposed in [START_REF] Lighthill | On kinetic waves. II. Theory of Traffic Flows on Long Crowded Roads[END_REF][START_REF] Richards | Shock Waves on the Highway[END_REF] as a single conservation law with unknown the vehicle density. This model based on a first order Partial Differential Equation (PDE) is very simple and robust but it fails to recapture some empirical features of traffic. In particular, it does not allow to take into account non-equilibrium traffic states mainly in congested situation. More sophisticated models referred to as higher order models were developed to encompass kinematic constraints of real vehicles or also the wide variety of driver behaviors, even at the macroscopic level. In this paper we deal with models of the Generic Second Order Modeling (GSOM) family. Even if these models are more complicated to deal with, they permit to reproduce traffic instabilities (such as the so-called stop-and-go waves, the hysteresis phenomenon or capacity drop) which move at the traffic speed and differ from kinematic waves [START_REF] Zhang | A non-equilibrium traffic model devoid of gas like behavior[END_REF] (see also [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF] and references therein).

Before the wide propagation of internet handsets, traffic monitoring has mainly been built on dedicated infrastructure which imply quite important installation and maintenance costs. Traffic flow monitoring and management has been deeply modified with the development of new technologies in mobile sensing aiming to provide a quite important quantity of floating car data. Traffic flow models are needed to be well suited such that managers could use both Eulerian and Lagrangian data for improving traffic state estimation. The term Eulerian refers to "classical" fixed equipment giving records of occupancy or flow of vehicles on a freeway section. This kind of measurements come from e.g. fixed inductive loop detectors, Radio Frequency Identification (RFID) transponders, radars or video cameras. By opposite, the term Lagrangian is used to characterize data coming from sensors which move within the measured field of interest. Lagrangian data are provided by on board mobile sensors such as Global Positioning Systems (GPS) or GPS-enabled smartphones.

While the idea of monitoring traffic using mobile sensors appeared less than ten years ago with the popularization of the mobile internet devices, there exists a fast growing literature about how to integrate Lagrangian data into classical macroscopic traffic flow models. The process of incorporating Eulerian and Lagrangian data into a mathematical model to improve the modeling is called data estimation or equivalently inverse modeling. According to the major UC Berkeley field experiment named Mobile Century and then Mobile Millennium [55] investigating Lagrangian sensing, it has been shown that even a 2% to 5% penetration rate of probe vehicles into the driver population provides sufficient and accurate data for estimating traffic velocity or density on highways [START_REF] Herrera | Incorporation of Lagrangian measurements in freeway traffic state estimation[END_REF][START_REF] Herrera | Evaluation of traffic data obtained via GPS-enabled mobile phones: The Mobile Century field experiment[END_REF][START_REF] Work | Lagrangian sensing: traffic estimation with mobile devices[END_REF]]. Nevertheless, it has been demonstrated in [START_REF] Piccoli | Second Order Models and Traffic Data from Mobile Phones[END_REF] that the quality of estimation for higher-order traffic quantities including vehicle acceleration/deceleration, emission and fuel consumption rates is dramatically affected when the penetration rate of probe vehicles or the sampling frequency of the current mobile sensors decrease. However on board devices propose a real breakthrough in traffic monitoring by providing a very cheap and efficient way to collect traffic data.

Motivation

In order to improve traffic states estimation from Lagrangian data, we propose to deal with macroscopic traffic flow models of the GSOM family. As these models combine the simplicity of the LWR model with the dynamics of driver specific attributes, we are able to recapture more specific phenomenon with a higher accuracy. While methods of data assimilation have been only developed for first order models up to now, this work presents a new algorithm to reconstruct traffic states from both Eulerian and Lagrangian data. We take advantage of a very recent article [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF] in which a variational principle has been proved for models of the GSOM family.

Organization of the paper

The rest of this paper is structured as follows. Section 2 presents more in detail the GSOM models and sheds a specific light on the LWR model which is widely used in traffic engineering. The variational principles for the GSOM models including LWR model are briefly recalled in Section 3. Section 4 is devoted to the presentation of the main elements of our computational method. Finally, Section 5 proposes a numerical example of our method.

GSOM traffic flow models 2.1 Formulation of GSOM models

In [START_REF] Lebacque | Second order traffic flow modeling: supply-demand analysis of the inhomogeneous Riemann problem and of boundary conditions[END_REF][START_REF] Lebacque | Generic second order traffic flow modelling[END_REF], the authors introduce a general class of macroscopic traffic flow models called the Generic Second Order Models (GSOM) family. Any model of the GSOM family can be stated in conservation form as follows

     ∂ t ρ + ∂ x (ρv) = 0 Conservation of vehicles, ∂ t (ρI) + ∂ x (ρvI) = ρϕ(I)
Dynamics of the driver attribute

I, v = I(ρ, I) Fundamental diagram, (2.1)
where ρ stands for the density of vehicles, v for the flow speed (equal to the mean spatial velocity of vehicles), x and t for position and time. The variable I is a specific driver attribute which can represent for example the driver aggressiveness, the driver destination or the vehicle class. The flow-density fundamental diagram (FD) is defined by

F : (ρ, I) → ρI(ρ, I).
Notice moreover that it was shown in [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF] that the notions of Supply and Demand functions defined in [START_REF] Lebacque | The Godunov scheme and what it means for first order traffic flow models[END_REF] for the classical LWR model could be extended to the GSOM family. The GSOM models admit two kinds of waves:

• Kinematic waves or 1-waves as in the seminal LWR model: a wave propagates density variations at speed ν = ∂ ρ I(ρ, I) while the driver attribute I is continuous across such a wave.

• Contact discontinuities or 2-waves: a wave propagates variations of driver attribute I at speed ν = I(ρ, I) while the flow speed is constant across such a wave.

Examples

The GSOM family recovers a wide range of existing models:

• The LWR model [START_REF] Lighthill | On kinetic waves. II. Theory of Traffic Flows on Long Crowded Roads[END_REF][START_REF] Richards | Shock Waves on the Highway[END_REF] itself is simply a GSOM model with no specific driver attribute, expressed as follows

∂ t ρ + ∂ x (ρv) = 0 Conservation of vehicles, v = I(ρ, x) Fundamental diagram. (2.
2)

The fundamental diagram (FD) for the LWR model states that traffic flow is always at an equilibrium state. It is commonly assumed that the flow is an increasing function of density between zero (corresponding to an empty section) and a critical density and then the flow decreases until the maximal density (corresponding to a bumper-to-bumper situation). However the FD shape is always a subject of debates (see for instance [START_REF] Fan | A comparison of data-fitted first order traffic models and their second order generalizations via trajectory and sensor data[END_REF]) and there exists a wide variety in the literature encompassing concave and triangular flow functions (see Figure 1 and also Chapter 3 of [START_REF] Garavello | Traffic flow on networks[END_REF] for additional examples). • The LWR model with bounded acceleration proposed in [START_REF] Lebacque | A two phase extension of the LWR model based on the boundedness of traffic acceleration[END_REF][START_REF] Lebacque | Two-phase bounded-acceleration traffic flow model: analytical solutions and applications[END_REF][START_REF] Leclercq | Bounded acceleration close to fixed and moving bottlenecks[END_REF] is also a GSOM model in which the propagated driver attribute is simply the speed of vehicles.

• The ARZ model (standing for Aw, Rascle [START_REF] Aw | Resurrection of "Second Order" Models of traffic flow[END_REF] and Zhang [START_REF] Zhang | A non-equilibrium traffic model devoid of gas like behavior[END_REF]) for which the driver attribute is taken as

I = v -V e (ρ) that is I(ρ, I) = I + V e (ρ).
• The Generalized ARZ model proposed in [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF] that can be also seen as a particular case of the model described in [START_REF] Zhang | A conserved higher-order anisotropic traffic flow model: description of equilibrium and non-equilibrium flows[END_REF]. These models introduce an interaction mechanism between two different FDs for distinguish equilibrium and non-equilibrium states.

• Multi-commodity models (multi-class, multi-lanes) of Jin and Zhang [START_REF] Jin | Multicommodity kinematic wave simulation model for network traffic flow[END_REF], Bagnerini and Rascle [START_REF] Bagnerini | A multiclass homogenized hyperbolic model of traffic flow[END_REF] or Herty, Kirchner, Moutari and Rascle [START_REF] Herty | Multicommodity flows on road networks[END_REF]. It encompasses also the model of Klar, Greenberg and Rascle [START_REF] Klar | Congestion on multilane highways[END_REF].

• The Colombo 1-phase model deduced in [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF] from the 2-phase model of Colombo [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF]. In this case, the driver attribute I is a scalar which is non-trivial in congested situation. In fluid area, the model follows the classical LWR model.

• The stochastic GSOM model of Khoshyaran and Lebacque [START_REF] Khoshyaran | A stochastic macroscopic traffic model devoid of diffusion[END_REF]. The driver attribute I is a random variable depending on the vehicle index N and on the random event ω such that I = I(N, t, ω). The random perturbations do not affect the vehicle dynamics but affect the driver perception and its behavior.

The interested reader is referred to [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF] and references therein for more details on examples.

Variational principles in traffic flow modeling

In traffic flow literature, the variational formulation was first conjectured by Newell in [START_REF] Newell | A simplified theory of kinematic waves in highway traffic, (I) General theory, (II) Queueing at freeway bottlenecks, (III) Multi-destination flows[END_REF] for the LWR first-order traffic flow model. It was then properly established and generalized by Daganzo in [START_REF] Daganzo | A variational formulation of kinematic waves: Basic theory and complex boundary conditions[END_REF][START_REF] Daganzo | A variational formulation of kinematic waves: Solution method[END_REF][START_REF] Daganzo | A variational formulation of kinematic waves: Bottleneck properties and examples[END_REF][START_REF] Daganzo | On the variational theory of traffic flow: well-posedness, duality and applications, Networks and Heterogeneous Media[END_REF]. In this section, we first present the variational formulation of the GSOM family lead by the ideas developed for the LWR model. Then the numerical methods for solving these variational formulations are introduced.

Variational formulation of the GSOM family

Guided by the ideas developed for the LWR model, variational formulations have been recently developed and proved for models of the GSOM family in both Eulerian [START_REF] Li | The Variational Formulation of a Non-equilibrium Traffic Flow Model: Theory and Implications[END_REF] and Lagrangian frameworks [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF].

Lagrangian setting

In [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF], the authors prove the existence of a variational principle for the GSOM models family expressed under its Lagrangian form, introducing r := 1/ρ the spacing and N the vehicle label

     ∂ t r + ∂ N v = 0 Conservation of vehicles, ∂ t I = ϕ(N, I, t)
Dynamics of the driver attribute I, v = V(r, I) = I(1/r, I) Fundamental diagram.

(3.3) Remark 3.1 (Vehicle labeling) Note that vehicles are labeled in the classical traffic sense, according to their passing time through a reference spatial point. Thus, the label axis is in the opposite sense that the spatial axis.

Considering the position X (N, t) = t -∞ v(N, τ )dτ , we obviously have that v = ∂ t X and r = -∂ N X . (3.4)
Hence the system (3.3) could be written as a Hamilton-Jacobi equation satisfied by X :

∂ t X -W(N, -∂ N X , t) = 0, (3.5)
where W denotes the speed-spacing fundamental diagram of vehicle N at time t. It is defined such that W(N, r, t) := V(r, I(N, t)) = I(1/r, I(N, t)),

where the driver attribute I(N, t) solves the following nonlinear first order ODE ∂ t I(N, t) = ϕ(N, I, t), I(N, 0) = i 0 (N ), for any N.

By classical results on optimal control problems i.e. dynamic programming on Hamilton-Jacobi-Bellman equations, we can check that:

X (N T , T ) = min u,(N 0 ,t 0 ) T t 0 M(N, u, t)dt + ξ(N 0 , t 0 ), Ṅ = u N (t 0 ) = N 0 , N (T ) = N T (N 0 , t 0 ) ∈ J (3.6)
where M is the Legendre-Fenchel transform of W according to its second variable (see Figure 2 for an illustration of the Legendre-Fenchel transform), defined such that

M(N, u, t) = sup r∈R {W(N, r, t) -ur}. (3.7)
M(N, p, t) pq W(N, q, t) Notice that the variational formulation of GSOM models (3.6) is only available if and only if r → W(N, r, t) is concave such that the Legendre-Fenchel transform of M gives back the function W.

W(N, r, t) q r p p
In the previous formulation (3.6), J is the locus of initial and boundary conditions and ξ(N 0 , t 0 ) is the initial data of the position of vehicle N 0 at time t 0 for any (N 0 , t 0 ) ∈ J . The optimal trajectories are given by the characteristics (Eulerian and Lagrangian characteristics are equal) described by the following coupled ODEs

Ṅ = ∂ r W(N, r, t), ṙ = -∂ N W(N, r, t), (3.8)
where u = Ṅ is the command of the optimal control problem (3.6). The system (3.8) is a simple system of ODEs in the (N, r) plane.

Eulerian setting

In [START_REF] Li | The Variational Formulation of a Non-equilibrium Traffic Flow Model: Theory and Implications[END_REF], the authors prove the existence of a variational formulation of a class of models from the GSOM family. Unlike the paper of Lebacque and Khoshyaran [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF], these models are expressed from the Eulerian viewpoint as a system of two conservation laws

∂ t ρ + ∂ x f (ρ, s) = 0 Conservation of vehicles, ∂ t s + ∂ x g(ρ, s) = 0
Dynamics around the equilibrium.

(3.9)

The first conservation law is obviously satisfied by the vehicle density ρ and the corresponding flux f is the product of density times speed f (ρ, v) = ρv. The second conservation law is satisfied by a variable s (a non-equilibrium measure) which can be replaced by s = ρI yielding the GSOM family without source terms ϕ(I) = 0. The corresponding flux g needs to satisfy some specific conditions to ensure that the wave speed is less or equal to traffic speed.

The paper [START_REF] Li | The Variational Formulation of a Non-equilibrium Traffic Flow Model: Theory and Implications[END_REF] highlights that taking ϕ(I) = 0 is equivalent to say that the driver attribute I is invariant along vehicle trajectories which seems to be corroborated for instance by the work of Duret et al. [START_REF] Duret | Estimating individual speed-spacing relationship and assessing ability of Newell's car-following model to reproduce trajectories[END_REF] on the NGSIM I-80 trajectories data-set (in congested situation).

For the models (3.9), the authors show in [START_REF] Li | The Variational Formulation of a Non-equilibrium Traffic Flow Model: Theory and Implications[END_REF] a variational formulation based on Lax-Hopf like formulas for both N ρ and N s defined as the cumulative quantity of respectively ρ and s.

When taking into account a non trivial relaxation term (ϕ(I) = 0), the authors assume that the problem reduces to solving a LWR model in large time because the relaxation term induces an exponential decay in time of the difference between (3.9) and a regular LWR model (2.2).

Remark 3.2 (Variational formulations of the LWR model) The variational theory for the LWR relies on the three dimensional representation of traffic flow or the so-called Moskowitz function [START_REF] Moskowitz | Notes on freeway capacity[END_REF][START_REF] Makigami | Three-dimensional representation of traffic flow[END_REF] yielding Hamilton-Jacobi formulations in Eulerian and Lagrangian coordinates (see also [START_REF] Laval | The Hamilton-Jacobi partial differential equation and the three representations of traffic flow[END_REF] and references therein). Under suitable assumptions, the solutions to the HJ equations could be determined thanks to a variational formulation, known as Lax-Hopf formula in the mathematics literature [START_REF] Evans | Partial Differential Equations: Second Edition[END_REF].

There exist different methods to obtain the variational formulation of these problems, including calculus of variations, dynamic programming [START_REF] Daganzo | A variational formulation of kinematic waves: Basic theory and complex boundary conditions[END_REF][START_REF] Daganzo | A variational formulation of kinematic waves: Solution method[END_REF] and viability theory [START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part I: Theory[END_REF][START_REF] Han | Lagrangian-based hydrodynamic model: Freeway traffic estimation[END_REF].

Notice that from a mathematical point of view it was rigorously established in [START_REF] Han | Lagrangian-based hydrodynamic model: Freeway traffic estimation[END_REF] that the (viscosity) solutions of HJ equations in Eulerian and Lagrangian coordinates are strictly equivalent.

Computational methods for GSOM models

Lax-Hopf algorithms in LWR case

Up to our best knowledge, existing numerical procedures for data assimilation use the seminal paper [START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part I: Theory[END_REF] which proposes a semi-explicit form of the solution to Hamilton-Jacobi equations with concave flow-density fundamental diagrams and with any piecewise affine (PWA) initial and boundary conditions. They introduce a generalized Lax-Hopf formula and the inf-morphism property [START_REF] Aubin | Dirichlet problems for some Hamilton-Jacobi equations with inequality constraints[END_REF][START_REF] Aubin | Viability theory: new directions[END_REF] to compute the solution by taking the infimum of all solutions associated with simpler partial problems. It follows computational methods that are commonly referenced as "grid free" schemes [START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part II: Computational methods[END_REF][START_REF] Mazaré | Analytical and grid-free solutions to the Lighthill-Whitham-Richards traffic flow model[END_REF]. Unlike dynamic programming methods [START_REF] Daganzo | A variational formulation of kinematic waves: Basic theory and complex boundary conditions[END_REF], Lax-Hopf algorithms have been proved to be exact in general cases of concave flow-density fundamental diagrams (i.e. for convex Hamiltonians) [START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part I: Theory[END_REF]. Some elements of comparison are given in [START_REF] Mazaré | Analytical and grid-free solutions to the Lighthill-Whitham-Richards traffic flow model[END_REF] between dynamic programming, Lax-Hopf algorithm and more classical methods such as Godunov scheme [START_REF] Daganzo | The cell transmission model, part II: network traffic[END_REF][START_REF] Lebacque | The Godunov scheme and what it means for first order traffic flow models[END_REF] and wave tracking algorithm [START_REF] Henn | A wave-based resolution scheme for the hydrodynamic LWR traffic flow model[END_REF][START_REF] Holden | Front tracking for hyperbolic conservation laws[END_REF]. The Lax-Hopf algorithm has been extended to the Lagrangian version of the LWR problem in [START_REF] Han | Lagrangian-based hydrodynamic model: Freeway traffic estimation[END_REF] for a triangular flow-density FD and for piecewise affine initial and boundary conditions.

Review of computational methods for GSOM models

There already exists some works on computational methods for models of the GSOM family [START_REF] Lebacque | Riemann problem resolution and Godunov scheme for the Aw-Rascle-Zhang model[END_REF][START_REF] Li | The Variational Formulation of a Non-equilibrium Traffic Flow Model: Theory and Implications[END_REF][START_REF] Qiu | Exact and grid-free solutions to the Lighthill-WhithamRichards traffic flow model with bounded acceleration for a class of fundamental diagrams[END_REF]. However these existing methods are developed in the Eulerian framework which does not seem to be the best one to deal with Lagrangian data. Moreover the paper [START_REF] Qiu | Exact and grid-free solutions to the Lighthill-WhithamRichards traffic flow model with bounded acceleration for a class of fundamental diagrams[END_REF] deals with the LWR model with bounded acceleration which is a very specific GSOM model. The algorithm is very similar to the ones developed in [START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part II: Computational methods[END_REF][START_REF] Mazaré | Analytical and grid-free solutions to the Lighthill-Whitham-Richards traffic flow model[END_REF] but it is not applicable to general models of the GSOM family.

In [START_REF] Li | The Variational Formulation of a Non-equilibrium Traffic Flow Model: Theory and Implications[END_REF] the authors present a numerical scheme that is based on a discretization of Eulerian time and space domain. The idea of the scheme (close to dynamic programming) is then to update the cumulative quantities

N ρ := +∞ x ρ(y, t)dy and N s := +∞ x s(y, t)dy,
on the vertexes of a mesh (whatever its shape if this mesh is dense enough) following the optimal paths in an iterative way because of the coupling of the equations to solve. The values of ρ (resp. s) are deduced from a first order approximation of the derivative of N ρ (resp. N s ). Notice that this computational method is not exact in general.

A finite difference scheme

A numerical scheme based on a uniform grid of the Lagrangian domain with steps ∆t and ∆n, was presented in [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF].

We define for any n, t ≥ 0 X t n := X (n∆n, tδt). As usual with finite difference scheme, we need to introduce a Courant-Friedrichs-Lewy (CFL) condition to ensure the scheme to be monotone and convergent ∆n ∆t ≥ sup

N,r,t |∂ r W (N, r, t)| .
Then we consider the classical first order finite difference scheme as follows

X t+1 n = X t n + ∆t W n∆n, X t n-1 -X t n ∆n , t∆t . (3.10)
By construction the above scheme can be interpreted as the seminal Godunov (finite volume) scheme (3.11) for the Lagrangian formulation of the GSOM model (see [START_REF] Costeseque | A convergent scheme for Hamilton-Jacobi equations on a junction: application to traffic[END_REF])

r t+1 n = r t n + ∆t ∆n W t n-1 -W t n , (3.11)
where we have defined the discrete spacing and respectively the numerical speed as follows

       r t n := X t n-1 -X t n ∆n , W t n := W n∆n, r t n , t∆t .
The upstream and downstream boundary conditions for the finite difference scheme (3.10) are fully described in [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF]. They match the Bardos-LeRoux-Nédélec or Dubois-LeFloch boundary conditions and it was already shown that such conditions are equivalent to prescribe supply and demand conditions (see [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF] and references therein).

Up to our best knowledge, there is no other existing work on computational methods for models of the GSOM family in the Lagrangian framework. The method of finite differences does not enjoy the semi-analytical expression of the solution of the HJ equation (3.5) thanks to Lax-Hopf formula (4.12). In the next section, we provide a semi-analytical algorithm to compute the solution based on that representation formula.

A "grid free" scheme

We are now interested in numerical methods to efficiently recover the solution of the Hamilton-Jacobi problem (3.5). The computation of numerical solutions of the HJ equation has already attracted an important interest in the mathematical community. The majority of numerical schemes which were proposed to solve HJ equations are based on finite difference methods, Semi-Lagrangian schemes and discontinuous Galerkin methods (the interested reader is referred to [START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part II: Computational methods[END_REF][START_REF] Costeseque | A convergent scheme for Hamilton-Jacobi equations on a junction: application to traffic[END_REF] and references therein).

Computational strategy

The idea to compute the solution of the Eulerian GSOM model (3.9) is to recast it under its Lagrangian form (3.3). As it was shown in the previous section, the position X of vehicle N at time t solves the Hamilton-Jacobi equation (3.5). As explained in [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF], the Initial and Boundary Value Problem (3.5) admits a quite simple representation formula (3.6), very similar to the Hopf-Lax formula presented for the LWR model. The Hopf-Lax formula (3.6) can be simplified because it is well-known from optimal control theory, that the optimal trajectories for which the minimum is attained are the characteristics (3.8). Hence we have to solve the following system of coupled ODEs (3.8). Then the generalized Hopf-Lax formulation (3.6) can be recast as follows

X (N T , T ) = min (N 0 ,r 0 ,t 0 ) T t 0 M(N, ∂ r W(N, r, t), t)dt + ξ(N 0 , t 0 ), Ṅ (t) = ∂ r W(N, r, t) ṙ(t) = -∂ N W(N, r, t) N (t 0 ) = N 0 , r(t 0 ) = r 0 , N (T ) = N T (N 0 , r 0 , t 0 ) ∈ K (4.12)
where K is the set of initial/boundary values obtained by combining the initial values of J with the initial/boundary value r 0 deduced from the initial values ξ(N 0 , t 0 ). According to the principle of inf-morphism property [START_REF] Aubin | Dirichlet problems for some Hamilton-Jacobi equations with inequality constraints[END_REF][START_REF] Aubin | Viability theory: new directions[END_REF][START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part I: Theory[END_REF][START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part II: Computational methods[END_REF], if the initial/boundary condition data ξ are located on a union (non necessarily disjoint) of sets

K = l K l ,
it suffices to solve partial problems on each set K l and to compute the minimum of the solutions to these sub-problems

X (N T , T ) = min l X l (N T , T ), (4.13)
with

X l (N T , T ) := min (N 0 ,r 0 ,t 0 ) T t 0 M(N l , ∂ r W(N l , r l , t), t)dt + ξ(N 0 , t 0 ). Ṅl (t) = ∂ r W(N l , r l , t) ṙl (t) = -∂ N W(N l , r l , t) N l (t 0 ) = N 0 , r l (t 0 ) = r 0 , N l (T ) = N T (N 0 , r 0 , t 0 ) ∈ K l (4.14)
In the remaining of this article, we will apply the inf-morphism property by considering initial and boundary conditions which are piecewise affine (PWA). Thus we will calculate the solution generated by each piece using (4.14) and then apply (4.13) in order to obtain the solution of (3.6).

Algorithm for piecewise affine value conditions

Hereafter, we study separately the different elements which contribute to the value of the solution of (3.5) in its Lagrangian setting. We distinguish

• the initial condition at time t = t 0 describing the initial position of vehicles ξ(N, t 0 ) for any considered N ,

• the "upstream" boundary condition that is the trajectory of the first vehicle N = N 0 traveling on the section ξ(N 0 , t) for any considered t,

• and internal boundary conditions given for instance by cumulative vehicle counts at fixed location X = x 0 .

Note that Lagrangian data which are individual vehicle trajectories are considered as a particular case of "upstream" boundary condition for a given N ≥ N 0 (see Remark 4.5).

While the Lax-Hopf algorithm can handle infinite horizon problems either in the Eulerian or in the Lagrangian framework, we restrict ourselves to finite values for a convenient numerical implementation. Without loss of generality, we assume that N ∈ [N 0 , N max ] with N max < +∞ and similarly t ∈ [t 0 , t max ] with t max < +∞.

We recall that according to Remark 3.1, in all the following figures, the N -axis and x-axis are increasing in opposite directions. Hence, the highest labels actually match vehicles that are located further upstream.

Initial conditions

In this case, at t = t 0 , the positions ξ(n, t 0 ) of vehicles n are given. We have that r 0 (N ) = -∂ N ξ(N, t 0 ), for any N.

The initial conditions for the characteristics are the couples (N, r 0 (N )). First, we need to discretize the set of vehicle labels into intervals [n p , n p+1 ] of length ∆n, for p = 1, ..., P , in such a way that the dynamics ϕ of the driver attribute I can be approximated by

ϕ p in the interval [n p , n p+1 ] ϕ(N, I, t) = ϕ p (I, t), for any N ∈ [n p , n p+1 ].
If the discrete step ∆n is small enough, we could also assume that for any N ∈ [n p , n p+1 ] the initial data are piecewise constant

I(N, t 0 ) = I 0,p , r(N, t 0 ) = r 0,p .
We will now first define properly the initial condition and then we calculate the solution generated by the p th component (p = 1, ..., P ).

Definition 4.1 (PWA initial condition)

Let t 0 ≥ 0 be fixed. Then the p th component of the initial condition (p = 1, ..., P ) is given by

X ini (N, t 0 ) = r 0,p (N -n p ) + α p , for any N ∈ [n p , n p+1 ].
To ensure continuity of the initial data X ini on [N 0 , N max ], we require that

α p = p-1 l=1
r 0,l ∆n, for any p = 1, ..., P.

In a first step, we want to compute the characteristics generated by the initial conditions (I(N, t 0 ) = i 0,p and r(N, t 0 ) = r 0,p given on the interval [n p , n p+1 ] at time t 0 ). In the general case, we would have

I(N, t) = I p (t) for any N ∈ [n p , n p+1 ].
In the following, we consider a characteristic denoted by N (t) for t ∈ [t 0 , t max ] and we distinguish two cases according to the location of this characteristic. To considerably simplify the presentation, we will restrict ourselves to the case of systems with no relaxation that is the dynamics of driver attributes ϕ p (I, t) = 0 at any time t ≥ t 0 and for any interval [n p , n p+1 ] with p = 1, ..., P .

If we assume that the dynamics is null, then it obviously follows that İp = 0 which leads to

I p (t) = I 0,p for any t ∈ Ω p . (4.16)
In traffic modeling, this choice is relevant if the driver attribute I represents for example some origin-destination (OD) information or if it characterizes the vehicle kind. The driver attribute is conserved at any time and along any characteristic inside each "strip" [n p , n p+1 ] × [t 0 , t max ] of width ∆n. By the way, the simplification (4.16) allows to consider characteristics which are straight lines while there are not in the general case.

At the edge N = n p .

The edge N = n p separates two different traffic states that exhibit two different speed-spacing fundamental diagrams. We need to take care of the kinematic waves that can appear at this locus. As illustrated on Figure 3, we set the characteristic speed related to the spacing r 0,p as follows ν p := ∂ r V (r 0,p , I 0,p ) , and we also define the speed of the "refracted" characteristic wave through the discontinuity of I such that

ν * p := ∂ r V r * 0,p , I 0,p . (4.17)
This characteristic speed is related to a new spacing value r * 0,p defined as follows

r * 0,p =    V -1 (V (r 0,p-1 , I 0,p-1 ) , I 0,p ) if V (r 0,p-1 , I 0,p-1 ) < sup r∈R V (r, I 0,p ) , +∞ else.
(4.18) that maps a given (initial) spacing r 0,p to its projection r q 0,p (which is the spacing value transported along the corresponding characteristic) in a strip [n q , n q+1 ] × [t 0 , t max ], with any q ≥ p and p ∈ {1, ..., P }. The transported value can be computed thanks to a recursive composition of function τ on consecutive strips r p 0,q := τ (r 0,p , q) = τ (τ (. . . τ (r 0,p , p + 1) , . . . , q -1) , q) , where r p+1 0,p := τ (r 0,p , p + 1) = r * 0,p+1 can be easily computed according to (4.18), for any p ∈ {1, ..., P }. Indeed, the speed is conserved through a discontinuity of I i.e.

n p I 0,p r 0,p r 0,p-1 I 0,p-1 ν p-1 ν * p N n p+1 n p-1 ν p
V r q 0,p , I 0,q = V r q+1 0,p , I 0,q+1 , for any q ≥ p.

Obviously we also have that r p 0,p := τ (r 0,p , p) = r 0,p .

Regarding the values of the initial spacing r 0,p and r 0,p-1 , we can distinguish two cases that could occur starting to an edge N = n p :

• either ν * p > ν p and in this case there will be a shock wave. The characteristics carrying respectively the initial states r 0,p-1 and r 0,p will cross each other, offering a partial superposition of both characteristics domains (see (a) on Figure 4). The solution is then simply computed thanks to the inf-morphism property. Hence we only consider the two characteristics (indicated as 1 and 2 on Figures 6, 7 and9) that encompass the domain of influence generated by the initial spacing r 0,p .

• or ν * p < ν p and in this case there will be a rarefaction fan. The characteristics waves carrying the initial states r 0,p-1 and r 0,p diverge (see (b) on Figure 4). The solution on the area between both extreme characteristics may not be computed and even the infmorphism property will not be able to recover it. That is why we need to consider the first two characteristics (as above) and a third wave corresponding to the first characteristic matching with the initial state r 0,p-1 (labeled 0 on Figures 6, 7 and9). Notice that by convention, for a given p = 1, ..., P , we add the rarefaction wave 0 separating the two states at N = n p (if any) to the solution generated by the initial condition on [n p , n p+1 ].

Remark 4.3 (Degenerate case for characteristics through a discontinuity of I)

When passing through a discontinuity of I (assume (a) and (b) the states on both sides of the discontinuity as illustrated on Figure 5), the characteristics speed may be changed. Nevertheless it is well known (Rankine-Hugoniot jump condition) that through such a discontinuity, the traffic speed should be unmodified

v (a) := V(r (a) , I (a) ) = V(r (b) , I (b) ) =: v (b) .
In case of v (a) > v max ) to catch up with the downstream vehicles (a). This induces the apparition of a vacuum area between both groups of vehicles. There is a rarefaction wave as vehicles (b) accelerate to attain their top speed and a shock wave following the last vehicle of type (a). Notice that the vacuum states (a ∞ ) and (b ∞ ) defined on Figure 5 coincide for r = +∞ (or equivalently at ρ = 0). It is noteworthy that the vacuum area is not visible in Lagrangian viewpoint (see Figure 5 (b)) while it matches the rarefaction fan and characteristics of type (a) in Eulerian viewpoint (see the area delimited by characteristic with speed v (b) max and dashdotted characteristic with speed v (a) on Figure 5 (d)).

The contribution X p of the initial condition defined on [n p , n p+1 ] × {t 0 } for any p ∈ {1, ..., P } is then computed as follows:

(i) Initialize the partial solution X p at +∞ on the whole computational domain [N 0 , N max ]× [t 0 , t max ].

(ii) Determine the number of characteristics to compute (two or three) according to what could occur at the edge N = n p .

(iii) Compute the equation N (t) of each characteristic while t ≤ t max and N ≤ N max .

(iv) Calculate the (exact) solution X p all along each computed characteristic generated by the interval [n p , n p+1 ] × {t 0 }, namely characteristics 1 and 2 , and 0 whenever it appears (see Figure 6). To this aim, we use the generalized Lax-Hopf formula (4.14) which gives that Ẋp (t) = M N (t), Ṅ (t), t .

r W (a ∞ ) (b ∞ ) (a) (b) v (b) v (a) v (b) max (b) (a) N t (a ∞ ) = (b ∞ ) (a) Speed-spacing fundamental diagram (b) Some characteristics in Lagrangian viewpoint ρ (b) (a ∞ ) = (b ∞ ) (a) v (b) max v (a) F 1/v (b) max t x ( a ∞ ) = ( b ∞ ) (a) 1/v (a) (b) (c) Flow-density fundamental diagram (d) Some characteristics in Eulerian viewpoint
The interested reader is referred to [START_REF] Brenier | Un algorithme rapide pour le calcul de transformées de Legendre-Fenchel discrètes[END_REF][START_REF] Corrias | Fast Legendre-Fenchel transform and applications to Hamilton-Jacobi equations and conservation laws[END_REF][START_REF] Lucet | Faster than the fast Legendre transform, the linear-time Legendre transform[END_REF] for additional information about the Legendre-Fenchel transform and fast algorithms for its numerical computation.

(v) Compute the exact value at any point within the characteristics strip (delimited by characteristics 1 and 2 on Figure 6) using the fact that for any point (N, t) belonging to the characteristic strip, the position at this point can be deduced by a simple translation of the position on characteristic 1 (see Figure 6). Indeed we have for any N ∈ [n q , n q+1 ] with q ≥ p and any t ∈ t

(2)

N , t (1)
N X p (N, t) = X p N, t (1)
N + t t (1) N V (r(N, τ ), I(N, τ )) dτ = X p N, t (1)
N + V r q 0,p , I 0,q t -t (1) N , (4.20)
where we recall that ϕ p (I, t) = 0 for any couple (I, t) and r q 0,p is computed according to (4.19). The time t (vi) In the case of a rarefaction fan, evaluate the value of X p at each point within the influence domain of the considered initial condition (illustrated on Figure 6) by an interpolation technique based on triangular meshes (the value at each triangle vertex is exact). for t = t 0 to t max do {Iteration on time} 10:

Compute N (j) (t) using (4.15) {Computation of the characteristic equation} 11:

Compute X p (N (j) (t), t) using (4.14) {Computation of the position along the characteristic} for t = t for N = N p to N max do {Iteration on greater labels} 21:

for t = t end for 25: end if Ensure:

X p (., .) over [N 0 , N max ] × [t 0 , t max ] r 0,p 0 I 0,p N 0 N n p n p+1 N max t t max

"Upstream" boundary conditions

In the Lagrangian setting, boundary conditions describe floating vehicle conditions in Eulerian setting. Indeed such a condition is equivalent to consider the trajectory t → ξ(N 0 , t) of a given vehicle N 0 . Then we have that v 0 (t) := ∂ t ξ(N 0 , t) = V(r 0 (N 0 , t), I(N 0 , t)), for any t ≥ t 0 , which can be solved and yields a unique solution (because V is strictly increasing) r 0 (t) = V -1 (v 0 (t), I(N 0 , t)), for any t ≥ t 0 .

The initial conditions for the characteristics are the couples (N 0 , r 0 (t)).

As for initial condition, we introduce a discrete time step ∆t and we consider a full discretization of the time domain [t 0 , t max ] into segments [t q , t q+1 ] of length ∆t with q = 1, ..., Q. We still consider that I(N, t) is piecewise constant w.r.t. N and constant w.r.t. time t (ϕ = 0) such that I(N 0 , t) = I 0,p=1 . If ∆t is small enough, we can consider that for any q ∈ {1, ..., Q}, we have v 0 (t) =: v 0,q , r 0 (t) = V -1 (v 0,q , I 0,p=1 ) =: r 0,q , for any t ∈ [t q , t q+1 ].

Definition 4.4 (PWA upstream boundary condition)

Let N 0 ∈ Z be fixed. Then the q th component of the upstream boundary condition (q = 1, ..., Q) is given by X up (N 0 , t) = v 0,q (t -t q ) + β q , for any t ∈ [t q , t q+1 ].

To ensure continuity of the upstream boundary data X ini on [t 0 , t max ], we require that

β q = q-1 l=1
v 0,l ∆t, for any q = 1, ..., Q.

The contributions of upstream boundary conditions are computed in a similar way than those of initial condition described in the previous subsection. The domain of influence of upstream boundary conditions is simply separated into what happens on each interval ]t q , t q+1 [ and at the edge t q for any q ∈ {1, ..., Q} (see Figure 7).

It is worth noting that we have a rarefaction wave at t = t q for any q ∈ {1, ..., Q} if and only if v 0,q-1 < v 0,q (or equivalently r 0,q-1 < r 0,q ). Note that in this case, the rarefaction wave denoted 0 separating the two states at t = t q (see Figure 7) is added by convention, to the solution computed from the boundary condition on [t q , t q+1 ].

N 0 t 0 N N max t t max 2 1 0 r 0,q t q+1 t q

Internal boundary conditions

Note that internal boundary condition is here understood in the Lagrangian framework. It does not correspond to internal vehicle trajectories which can also be incorporated in the algorithm (see the previous subsection and Remark 4.5). Notice that internal boundary condition is called mixed condition in [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF] to avoid any misunderstanding.

We assume that the data comes from vehicles N (t) located at a point ξ(N (t), t) at time t. Thus we get

d dt ξ(N (t), t) = -Ṅ (t)r 0 (t) + W(N (t), r 0 (t), t). (4.21)
As not all data are compatible (in the sense that they need to comply to the traffic dynamics laws), we need to introduce compatibility conditions. From Legendre-Fenchel transform, we have sup

r∈R -Ṅ(t)r + W(N (t), r, t) = M(N (t), Ṅ (t), t),
hence the equation yielding r 0 (t) (4.21) admits a solution only if

d dt ξ(N (t), t) ≤ M(N (t), Ṅ (t), t). (4.22)
Let us describe the specific situation of ξ(t) = 0 and Ṅ (t) ≥ 0 for any time t. This case occurs when the data originates at a fixed measurement point x 0 like a fixed detector data that is ξ(N (t), t) = x 0 for all time t. Then N (t) represents the cumulative flow at point x 0 . In this case the above compatibility condition (4.22) is satisfied if moreover we assume that Ṅ ≤ F max (I). It is quite natural to impose such a condition because Ṅ (t) is equal to the instantaneous traffic flow which is evidently bounded by F max (I) (see Figure 8). Moreover Algorithm 2 Pseudo-code for the computation of X q on the computational domain [N 0 , N max ] × [t 0 , t max ], under the upstream boundary condition X up . Require: N 0 , N max , t 0 , t max , ∆t, t q and r 0,q {Input label domain, time domain and initial condition} 1: X q ← +∞ {Initialization of the position function to infinity} 2: if r 0,q-1 > r 0,q then {Determination of the number of characteristic curves to compute} for t = t 0 to t max do {Iteration on time} 9:

Compute N (j) (t) using (4.15) {Computation of the characteristic equation} 10:

Compute X p (N (j) (t), t) using (4.14) {Computation of the position along the characteristic} 11:

end for 12: end for 13: for N = N p to N max do {Iteration on greater labels} 14:

for t = t

N to t (2)
N do {Iteration on time} for N = N p to N max do {Iteration on greater labels} 20:

for t = t (0) N to t (1)
N do {Iteration on time}

21:

Compute X (N, t) using interpolation on the three closest exact values {Computation of the position within the rarefaction fan} end for 24: end if Ensure: X q (., .) over [N 0 , N max ] × [t 0 , t max ] such a condition ensures that the Lagrangian stays non-negative which is coherent with its interpretation as a instantaneous speed cost. As illustrated on figure 8, there are two solutions (except in the very particular case of Ṅ = F max ), one under-critical and the other over-critical according to classical traffic definitions. Recall that the speed-spacing fundamental diagram W depends on the vehicle N . If we propagate both pairs of characteristics, then the inf-morphism property will automatically select a single solution which matches the over-critical situation (because the speed is lower than in the under-critical situation). It may happen that the congested solution is not the good one. Then to avoid any mistake, we assume that the fixed detector gives the cumulative flow and also the traffic flow speed. In this way, we can select one traffic state only.

W(N, r, t) In our case, we are interested in including some Eulerian data coming from classical fixed sensors like inductive loop detectors. To achieve this aim, Eulerian data become internal boundary condition into the Lagrangian framework. We assume that a Eulerian sensor located at a fixed position x 0 gives us the incremental cumulative vehicle count which is then interpolated in a piecewise affine function N (t) for t ∈ [t 0 , t max ]. We then define by f 0,p the value of Ṅ (t) for any t ∈ [t p , t p+1 ] and we set tp the time such that N ( tp ) = N p for any p ∈ {1, ..., P }.

ν under r under 0,p r over 0,p ν over r Ṅ > 0 0 ρ max F max Ṅ 1/r over 0,p F low, F 1
Definition 4.6 (PWA internal boundary condition) Let x 0 ≥ 0 be fixed. Then the p th component of the internal boundary condition (p = 1, ..., P ) is given by

X int (n, t) = x 0 , for n = f 0,p t -tp + γ p and t ∈ [ tp , tp+1 ],
where we define tp for any p as the time such that n = n p . To ensure continuity of the trajectory t → n(t) on which the internal boundary data X int is prescribed, we require that

γ p = p-1 l=1
f 0,l ( tl+1 -tl ), for any p = 1, ..., P.

We assume that N is piecewise affine on each discrete segments [n p , n p+1 ] that is equivalent to say that Ṅ is piecewise constant (i.e. constant in each strip [n p , n p+1 ] × [t 0 , t max ] for p = 1, ..., P ). It is easy to deal with this case in the algorithm because the computational steps are similar to both previous cases for initial and upstream boundary conditions. The only difference resides in the fact that we have a characteristic strip that matches either the under-critical or the over-critical traffic state. It is simple to verify that in the first (undercritical) case, the p th component generates characteristics which are emitted with speeds ν under > Ṅ while it is the inverse in the other case (see Figure 9). From a traffic point of view, it is relevant with observations stating that a deceleration wave spreads over vehicles located further upstream at a slower speed than an acceleration wave. Hence, from a computational viewpoint, if we are in the under-critical (resp. over-critical) case, then we will use the same strategy than Algorithm 1 (resp. Algorithm 2).

r 0,p 1 0 2 n p n p+1 N 0 N N max t t max t 0 r 0,p 2 0 1 n p n p+1 N 0 N N max t t max t 0

Recapitulation of the overall algorithm

The different elements of the computations for piecewise affine value conditions are recapitulated in the following pseudo-code (see Algorithm 3).

Remark 4.7 This pseudo-code underlines the property of the Lax-Hopf algorithm to treat in a parallel way the different value conditions (initial, upstream boundary and internal boundary).

5 Numerical example

Instantiation

In order to simplify the computations and to ease the presentation of the following example, we consider that the driver attribute I is piecewise constant with respect to N at initial time t = t 0 . It means that there exist some platoons of vehicles which share the same driver attribute. It could be for example some vehicles of the same kind (cars and trucks), or vehicles that go to the same destination, or vehicles that have the same desired maximal speeds. Actually, in this numerical example, we have considered ∆N = 1 so the attribute I is constant for each vehicle. As ϕ = 0, this attribute is conserved at any time t ≥ t 0 (see Figure 11 (left)). We consider for this numerical example a Colombo 1-phase model [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF] given by

W(r, N, t) :=      V - β r if r ≥ r crit (I), (I + q * r) 1 - 1 rR else,
where β, q * , V and R are given parameters of the model. More precisely, q * stands for the theoretical maximal flow, while V and R describe respectively the maximal speed and density. The critical spacing r crit (I) separating the free and congested cases is computed as the inverse of the critical density

ρ crit (I) := 1 2 β - I R V + q * R -I - V + q * R -I 2 -4q * β - I R .
Algorithm 3 Pseudo-code implementation for the Lax-Hopf based computation of the position function X on the computational domain [N 0 , N max ] × [t 0 , t max ] prescribed by the user. Require: N 0 , N max , t 0 and t max {Input label domain, time domain} 1: X ← +∞ {Initialization of the position function to infinity} 2: for p = 1 to P do {Iteration on the initial conditions} for T = t 0 to t max do {Iteration on time} 4:

for N T = N p to N max do {Iteration on greater labels} 5:

compute X ini p (N T , T ) using Algorithm 1 {Component induced by the initial condition r 0,p }

The numerical values considered in this example are given below

           V = 25 m.s -1 V = 22 m.s -1 R = 1/5 veh.m -1 R = 1/30 veh.m -1 and    β = V - V R = 90 m 2 .s -1 q * = 1 veh.s -1
The flow-density and speed-spacing fundamental diagrams are given on Figure 10. Notice that these numerical values are satisfactory for a two-lanes road section. The driver attribute I varies between 0 and I max = 5 which is designed such that the flow-density fundamental diagram F(ρ, I) stays concave with respect to ρ.

g 1 :=        g 1 (r crit (I)) if β p < r crit (I),
g 2 :=                g 2 (r min ) if q * < p and I R(p -q * ) < r min , g 2 (r crit (I)) if q * ≥ p or if q * < p and I R(p -q * ) > r crit (I), g 2 β p if q * < p and r min ≤ I R(p -q * ) ≤ r crit (I).
The Lagrangian M is then computed for our numerical instantiation, as illustrated on Figure 11 (right). We consider piecewise affine initial conditions that is we prescribe each vehicle position X(t, n) at time t = t 0 for n ∈ [N 0 , N max ] according to a normal distribution of spacing between r min = 5 m and a 50 m "maximal" spacing. The initial condition for this numerical test are plotted on Figure 12. Moreover we prescribe piecewise affine upstream boundary condition that is the trajectory X(t, n) of the first vehicle N = N 0 for t ∈ [t 0 , t max ]. This upstream boundary condition is illustrated on Figure 13. We also have included two (fictive) vehicle trajectories for N = 25 and N = 46 as Lagrangian internal boundary conditions. They are plotted on Figure 14. For these trajectories, we have considered randomly distributed speed values that are piecewise constant. Note that the speed distribution is chosen such that the mean value is equal to 18 m.s -1 in a first time and then to 10 m.s -1 . In that way, one can verify the trajectories are piecewise affine and averaging on time, vehicles are globally decelerating.

Numerical result and interpretation

The solution is thus given by applying the computation algorithm described in Algorithm 3. The simulations deal with 50 vehicles on a 2 minutes run. We also compare the solution with initial and upstream boundary conditions and the solution with initial, upstream and internal boundary conditions. The numerical solutions X are respectively plotted on Figures 15 (a In Figures 15 (c) and (d), we have also plotted the speeds of vehicles which are computed as the discrete gradient of the solution X (N, t) according to time t and for any vehicle N . One can observe the shock and rarefaction waves generated from initial and boundary conditions. Note that the internal boundary conditions are not necessarily defined for every t ∈ [t 0 , t max ]. It is the reason why there is an area up to N = 46 and t = 110 that matches a high speed rarefaction fan since the solution is no more constrained by the internal boundary condition. One can note that the integration of two slower vehicles (N = 25 and N = 46) in the traffic flow implies the modification of the solution. The values of X for N ≥ 25 and then N ≥ 46 are notably decreased. Indeed slower vehicles impose lower speeds to the following vehicles as it could be noticed on vehicle trajectories plotted on Figure 16. Thus, data assimilation allows to modify the previous over-estimation we made by considering only initial and upstream boundary conditions. In the numerical example above, we do not consider real data assimilation. The next step for numerical results should be the assimilation of real data in a model of the GSOM family. As it was described in the previous section, vehicle trajectories and fixed detector measurements can be easily incorporated in the algorithm. However we need to identify a driver attribute I and to know how it evolves on the computational domain. That is the tough point because it is hard to measure or to evaluate. Notice that an example of speed data assimilation with the ARZ model which is a model of the GSOM family can be found in [START_REF] Lebacque | Riemann problem resolution and Godunov scheme for the Aw-Rascle-Zhang model[END_REF].

Concerning the identification of I, a method could be deduced from [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF], exploiting data-fitted fundamental diagrams. In the case of a scalar attribute I, we can use for instance speed and spacing measurements to evaluate the value of I.

Discussion and conclusion

In this paper, we are interested in the assimilation of different kinds of data into a generic class of macroscopic traffic flow models for improving traffic state estimation. The data come from mixed sources including Lagrangian vehicle trajectories and Eulerian cumulative vehicle counts obtained from fixed detectors. We then consider macroscopic traffic flow models of the Generic Second Order Modeling (GSOM) family in the Lagrangian system of coordinates which seems to be the most adapted framework for dealing with such data. We describe a computational method providing solutions to these models under piecewise affine initial and boundary conditions. The numerical method is based on the variational theory which has been extended to the GSOM family very recently [START_REF] Lebacque | A variationnal formulation for higher order macroscopic traffic flow model of the GSOM family[END_REF]. Solution can be computed thanks to Lax-Hopf like formula and a generalization of the inf-morphism property [START_REF] Claudel | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part I: Theory[END_REF]. It is possible to considerably reduce the number of integral curves by following only the characteristics which are the optimal trajectories. Another computational trick is to solve separately many partial problems (by discretizing initial and boundary conditions) instead of the general problem.

Extensions of this work could include more general assumptions on the dynamics ϕ = 0 of driver attributes, always with the assumption of piecewise affine value conditions. The computational benefits of the variational method over the finite difference method which is also presented in the paper need to be illustrated on several comparisons. We expect that the variational method has a lower computational cost (due to parallel computing e.g.) and a globally higher accuracy compared with the other method.

A still open but interesting question is the application of the algorithm to data assimilation on road networks. The challenging point is to deal with fundamental diagrams which will depend on the drivers attribute but also on the spatial position. To our best knowledge Lax-Hopf formula seems not to be so tractable in that framework because the formula does not simplify so much for space dependent Hamiltonians. One considered possibility is to use the extended concepts of supply and demand functions to models of the GSOM family [START_REF] Lebacque | An intersection model based on the GSOM model[END_REF]. This subject is being currently investigated by the authors.

0FFigure 1 :

 1 Figure 1: Illustrations of some flow functions F for the LWR model: Greenshields (left), triangular (center) and exponential (right).

Figure 2 :

 2 Figure 2: Legendre transform M of function W.

  On the interval N ∈]n p , n p+1 [. Consider Ω p := {t | N (t) ∈]n p , n p+1 [} for any p = 1, ..., P . We then have to solve the following nonlinear first order ODE İp (t) = ϕ p (I p (t)), for t ∈ Ω p , I p (t 0 ) = I 0,p . On the interval ]n p , n p+1 [, we have that I(N, t) = I p (t) is independent of N i.e. ∂ N I = 0. Then from (3.8), we deduce ṙ = -∂ N W(N, r, t) = -∂ I V(N, I(N, t)) ∂ N I(N, t) = 0. Thus r(N, t) on any characteristic in Ω p . Hence it simply suffices to solve the second ODE in (3.8) that is Ṅ = ∂ r W(N, r, t) to get the equation of the characteristic on the interval ]n p , n p+1 [. We finally get N (t) = t u ∂ r W(N, r, τ )dτ + N (u), for any t, u ∈ Ω p . (4.15)

Figure 3 :

 3 Figure 3: Schematic view of what happens for the characteristic wave generated from the edge N = n p and passing through a discontinuity of I.

Remark 4 . 2 (

 42 Transported spacing value along a characteristic)We extend the idea of what happens at each discontinuity of I by introducing a function τ : (r 0,p , q) → r q

Figure 4 :

 4 Figure 4: Characteristics through a discontinuity of I of r.

  , I (b) ), it is obvious that the equality between these two speeds cannot be complied. In an Eulerian setting, it means that the group of upstream vehicles (b) cannot accelerate enough (even increasing their speeds to v (b)

Figure 5 :

 5 Figure 5: Critical "vacuum" case appearing from special condition values.

N

  , i = 1, 2 corresponds to the time when the characteristic i crosses the line {N } × [t 0 , t max ].

  end for 14: for N = N p to N max do {Iteration on greater labels} 15:

{Iteration on time} 16 :

 16 Compute X (N, t) using(4.20) {Computation of the position within the characteristics strip} 17: end for 18: end for 19: if n charac = 3 then {Case of the rarefaction fan} 20:

{Iteration on time} 22 :

 22 Compute X (N, t) using interpolation on the three closest exact values {Computation of the position within the rarefaction fan}

Figure 6 :

 6 Figure 6: Domain of influence of the initial data [n p , n p+1 ] × {t 0 }.

Figure 7 :

 7 Figure 7: Domain of influence of the upstream boundary data {N = N 0 } × [t q , t q+1 ].

Remark 4 . 5 (

 45 Extension to any individual trajectory condition) Any Lagrangian data giving the trajectory of a vehicle N * with N * ∈ [N 0 , N max ] will strictly generate the same type of calculations than the ones described for the upstream boundary condition. It is then easy to deal with data coming from mobile sensors moving within the traffic stream.

  end if 7: for j = 1 to n charac do {Iteration on the characteristics} 8:

15 :

 15 Compute X (N, t) using(4.20) {Computation of the position within the characteristics strip} 16: end for 17: end for 18: if n charac = 3 then {Case of the rarefaction fan} 19:

Figure 8 :

 8 Figure 8: Existence of two solutions corresponding to a condition ξ(t) = 0 and Ṅ (t) ≥ 0.

Figure 9 :

 9 Figure 9: Domain of influence of an internal boundary data (in red).

  veh/km) Flow F (veh/h) Fundamental diagram F(ρ,I)

Figure 10 :g 1

 101 Figure 10: Speed-spacing fundamental diagram V(r, I) (left) and flow-density fundamental diagram F(ρ, I) (right).

r

  min ≤r≤r crit (I)

Figure 11 :

 11 Figure 11: Initial values of driver attribute I(N, t 0 ) (left) and Lagrangian function M(N, u, t) (right).

Figure 12 :

 12 Figure 12: Initial conditions for the GSOM PDE at t = t 0 .

  )

Figure 13 :Figure 14 :

 1314 Figure 13: Boundary conditions for the GSOM PDE at N = N 0 .

Figure 15 :

 15 Figure 15: Numerical solution for GSOM PDE obtained on the computation domain [N 0 , N max ] × [t 0 , t max ].

Figure 16 :

 16 Figure 16: Vehicle trajectories obtained for integers of the GSOM numerical solution X(t, n).

  Moskowitz surface X(t, n) (b) Speeds

Figure 17 :

 17 Figure 17: Numerical solution for GSOM PDE obtained on the computation domain [N 0 , N max ] × [t 0 , t max ] with Lagrangian and Eulerian data assimilation.

  Algorithm 1 Pseudo-code for the computation of X p on the computational domain [N 0 , N max ] × [t 0 , t max ], under the initial condition X ini . Require: N 0 , N max , t 0 , t max , ∆N , N p and r 0,p

			{Input label domain, time domain and
		initial condition}
	1: X p ← +∞	{Initialization of the position function to infinity}
	2: Compute µ * p and r * 0,p according to (4.17)-(4.18) 3: if µ * p > µ p then {Determination of the number of characteristic curves to compute}
	4:	n charac = 2
	5: else
	6:	n charac = 3
	7: end if
	8: for j = 1 to n charac do	{Iteration on the characteristics}
	9:

Aknowledgements

The first author would like to thank Jia Li, PhD candidate at UC Davis, California, USA, for deeper insight of his work. This work was partially supported by the ANR (Agence Nationale de la Recherche) through HJnet project ANR-12-BS01-0008-01.

6:

if X ini p (N T , T ) < X then 7:

X ← X ini p (N T , T ) {Update the position function} end for 11: end for 12: for q = 1 to Q do {Iteration on upstream boundary conditions} for T = t q to t max do {Iteration on greater times} 14:

for N T = N 0 to N max do {Iteration on labels} 15:

compute X up q (N T , T ) using Algorithm 2 {Component induced by the upstream boundary condition X N 0 , t q } 16:

if X up q (N T , T ) < X then 17:

X ← X up q (N T , T ) {Update the position function} for T = t p to t max do {Iteration on greater times} 24:

for N T = N p to N max do {Iteration on greater labels}