Guillaume Costeseque

Jean-Patrick Lebacque

Asma Khelifi

Lagrangian GSOM traffic flow models on junctions

Keywords: Traffic flow, junction, GSOM models, Lagrangian, numerical scheme

This paper is concerned with the macroscopic modeling and simulation of traffic flow on junctions. More precisely, we deal with a generic class of second order models, known in the literature as the GSOM family. While classical approaches focus on the Eulerian point-of-view, here we recast the model using its Lagrangian coordinates and we treat the junction as a specific discontinuity in Lagrangian framework. We propose a complete numerical methodology based on a finite difference scheme for solving such a model and we provide a numerical example.

INTRODUCTION

Motivation

In this paper, we are motivated by road network modeling, thanks to macroscopic traffic flow models. First order traffic flow models have been used for quite a long time for modeling traffic flows on networks (see [START_REF] Garavello | Traffic flow on networks[END_REF]; [START_REF] Lebacque | First order macroscopic traffic flow models for networks in the context of dynamic assignment[END_REF] for instance). In particular, the seminal LWR model standing for [START_REF] Lighthill | On kinematic waves II. A theory of traffic flow on long crowded roads[END_REF]; [START_REF] Richards | Shock waves on the highway[END_REF] has been widely used. However, first order models do not allow to recapture accurately specific and meaningful traffic flow phenomena. Thus we focus on the Generic Second Order Models (GSOM) family which encompasses a large variety of higher order traffic flow models. GSOM models have been already well studied on homogeneous sections but they have attracted little attention for their implementation on junctions, as it is discussed in Section 3. However, junctions are the main source of congestion for traffic streams on a network.

In this paper, we want to develop a junction model which is compatible with microscopic and macroscopic descriptions, and satisfies classical constraints coming from engineering, as for instance the invariance principle discussed in [START_REF] Lebacque | First-order macroscopic traffic flow models: Intersection modeling, network modeling[END_REF]; [START_REF] Tampère | A generic class of first order node models for dynamic macroscopic simulation of traffic flows[END_REF]; Costeseque and Lebacque (2014a). The microscopic representation of traffic flow is particularly suited for traffic management methods, while staying compatible with a macroscopic representation allowing global evaluation. The key idea for conciliating both microscopic and ⋆ This work was partially supported by the ANR (Agence Nationale de la Recherche) through HJnet project ANR-12-BS01-0008-01. macroscopic representations is to recast the macroscopic model under its Lagrangian coordinates. Indeed the Lagrangian framework focuses directly on the particles and incidentally it allows to keep track of individual behaviors (see for instance [START_REF] Leclercq | The Lagrangian coordinates and what it means for first order traffic flow models[END_REF] in the case of first order LWR model).

Organization of the paper

The article is organized as follows. In Section 2, the generic class of second order macroscopic traffic flow models called GSOM family, is introduced. In Section 3, we review and discuss the literature for solving GSOM models posed on junctions. Our aim is to show that the Lagrangian framework is well-suited for designing the solution to GSOM problems even if incorporating discontinuities. The complete numerical methodology is described in Section 4 and a numerical example is given incidentally in Section 5. Finally, we provide some conclusions on this work and give some insights on future research in Section 6.

GSOM FAMILY

Formulation of GSOM models

In Lebacque et al. (2005[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. Assume that ρ(t, x) stands for the density of vehicles at location x ∈ R and time t ≥ 0. Any model of the GSOM family can be stated in conservation form as follows ∂ ∂t

ρ ρI + ∂ ∂x ρv ρvI = 0 ρϕ(I)
, on (0, +∞)×R, (1)

with the initial conditions ρ(0, x) = ρ 0 (x), on R,

I(0, x) = I 0 (x), on R.
The speed is defined by v := I(ρ, I) where I denotes the speed-density Fundamental Diagram (FD) which depends on the choice of the GSOM model. It is usually assumed to be non-increasing with respect to its first variable. The flow-density FD is then defined by

F : (ρ, I) → ρI(ρ, I).
(2)

The variable I(N (t, x)) ∈ R m , for some m ∈ N, represents an attribute which is specific to the vehicle N located at position x and time t. This attribute (which is later denoted I(t, x) by abuse of notation) can represent, for instance, the driver aggressiveness, the driver destination or the vehicle class. The function ϕ : I → ϕ(I) accounts for the dynamics of the attribute I.

The GSOM family recovers a wide range of existing models among which there is the LWR model of [START_REF] Lighthill | On kinematic waves II. A theory of traffic flow on long crowded roads[END_REF]; [START_REF] Richards | Shock waves on the highway[END_REF], a GSOM model with no specific driver attribute (I = κ for any κ ∈ R),

∂ t ρ + ∂ x (ρv) = 0 on (0, +∞) × R, ρ(0, x) = ρ 0 (x), on R (3)
and the ARZ model [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] and [START_REF] Zhang | A non-equilibrium traffic model devoid of gaslike behavior[END_REF]) for which the driver attribute is taken as the difference between actual speed and the equilibrium speed V e (ρ) depending only on the density, that is

I = v -V e (ρ).
The speed-density FD boils down to I(ρ, I) = I + V e (ρ).

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

A note on Supply-Demand functions

It is worth noticing 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 (3) and expanded to the case of the LWR model on junctions in [START_REF] Lebacque | First order macroscopic traffic flow models for networks in the context of dynamic assignment[END_REF], could be also extended to the GSOM family, as it was shown in [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF]. These functions built on the flowdensity FD are essential to build monotone finite volume schemes for solving the hyperbolic system (1). Supply and demand functions are also particularly relevant for traffic flow modeling through junctions.

If F denotes the flow-density FD defined in (2) (by extension of notation, we will consider F(., .; x) as the flowdensity FD at position x ∈ R), then equilibrium demand and supply functions are defined as follows

   ∆ (ρ, I; x) := max 0≤k≤ρ F k, I; x -, Σ (ρ, I; x) := max k≥ρ F k, I; x + (4)
where x + (resp. x -) denotes the upper (resp. lower) limit.

In the case of GSOM models, the extension of local traffic supply and demand definitions is far from being straightforward. Indeed, it has been pointed out in [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF] that downstream supply depends on the upstream driver attribute (which has not already passed through the considered point). The upstream demand and downstream supply at a point x and time t are defined such that δ(x, t) := ∆ ρ(x -, t), I(x -, t); x , σ(x, t) := Σ ρ(x, t), I(x + , t); x , (5) with ρ(x, t) := I -1 I ρ(x + , t), I(x + , t); x + , I(x -, t); x + . (6)

Then the passing flow at location x, time t is given by the celebrated Min formula q(t, x) = min {δ(t, x), σ(t, x)} .

I I(., I -) I(., I + ) v := I(ρ + , I + ) ρ := I -1 (v, I -) ρ + ρ
Fig. 1. The downstream supply depends on the upstream driver attribute. Here, we have defined

I + := I(x + , t), I -:= I(x -, t) and ρ + := ρ(x + , t).
Let us introduce the modified demand for the GSOM models as follows Ξ (ρ, I 1 , I 2 ; x) := Σ I -1 I ρ, I 2 ; x + , I 1 ; x + , I 2 ; x , (7) such that ( 5)-( 6) boil down to σ(x, t) := Ξ ρ(x, t), I(x + , t), I(x -, t); x .

Lagrangian setting of the GSOM family

The common expression of GSOM models in Eulerian coordinates (t, x) is given by (1). However, it is well-known that Lagrangian framework (t, n) is particularly convenient for dealing with flows of particles and it is especially true in traffic flow modeling (see [START_REF] Leclercq | The Lagrangian coordinates and what it means for first order traffic flow models[END_REF]; van [START_REF] Van Wageningen-Kessels | Discontinuities in the Lagrangian formulation of the kinematic wave model[END_REF] and references therein).

Assume that N (t, x) ∈ R describes the label of vehicle at position x and time t. We set the spacing r := 1 ρ and the speed-spacing FD as follows

V(r, I) := I 1 r , I for any (r, I) ∈ (0, +∞) × R m . (8) If we set the following change of coordinates    N (t, x) := +∞ x ρ(t, ξ)dξ, T := t such that ∂ N = -r∂ x , ∂ T = ∂ t + v∂ x
where v denotes the speed of particles, then the GSOM model ( 1) can be recast in Lagrangian form as follows

∂ T r + ∂ N V(r, I) = 0 on (0, +∞) × R, ∂ T I = ϕ(I) on (0, +∞) × R, (9)
with initial conditions

r(0, n) = r 0 (n), on R, I(0, n) = I 0 (n), on R.
One can notice also that we have r = -∂ N X where X (T, N ) denotes the position of particle N at time T which solves the following Hamilton-Jacobi equation

∂ T X = V (-∂ N X , I) , ∂ T I = ϕ(I).
(10)

BRIEF REVIEW OF THE LITERATURE

There already exist a few works on the Lagrangian modeling of junctions based on GSOM models; see for instance [START_REF] Moutari | A hybrid Lagrangian model based on the Aw-Rascle traffic flow model[END_REF]; van [START_REF] Van Wageningen-Kessels | Discontinuities in the Lagrangian formulation of the kinematic wave model[END_REF]; [START_REF] Khoshyaran | Lagrangian modelling of intersections for the GSOM generic macroscopic traffic flow model[END_REF]. However, some of these works are based on very specific examples extracted from the GSOM family. But it is not straightforward to extend the numerical methodologies presented in these papers to the generic GSOM model ( 9).

In [START_REF] Khoshyaran | Lagrangian modelling of intersections for the GSOM generic macroscopic traffic flow model[END_REF], the authors consider the Godunov scheme applied to (9). The authors extend this particle discretization to networks, addressing the problem of junction modeling through a supply-demand approach. The authors make the choice to introduce an internal state model (see [START_REF] Lebacque | An intersection model based on the GSOM model[END_REF]; Khoshyaran and Lebacque ( 2009)) and assume that the particles share the same attribute once they have passed. By the way, the authors deal also with densities and flows which is not particularly convenient with GSOM models in the Lagrangian framework. Hopefully, dealing with spacing instead of density will ease the resolution of the model.

While boundary conditions can be treated within the framework of supply-demand flows methodology (see Lebacque et al. (2005[START_REF] Lebacque | Generic second order traffic flow modelling[END_REF] and [START_REF] Khoshyaran | Lagrangian modelling of intersections for the GSOM generic macroscopic traffic flow model[END_REF]), expressions of upstream and downstream boundary conditions into Lagrangian coordinates can be obtained in the framework of variational approach for GSOM models (see [START_REF] Lebacque | A variational formulation for higher order macroscopic traffic flow models of the GSOM family[END_REF]). It will be developed in next section.

METHODOLOGY FOR THE LAGRANGIAN MODELING OF JUNCTIONS

In this section, we describe the numerical scheme adapted for the generic Lagrangian GSOM model ( 9) posed on a junction. We partially follow [START_REF] Lebacque | A variational formulation for higher order macroscopic traffic flow models of the GSOM family[END_REF] in which the authors describe boundary conditions for the Hamilton-Jacobi equation ( 10) associated to the Lagrangian GSOM model ( 9).

We set a junction as the union of N I incoming and N O outgoing branches that intersect at a unique point called the junction point (or the node in the traffic literature). We also define I := 1, N I (resp. J := N I +1, N I +N O ) the set of incoming (resp. outgoing) links.

Lagrangian discretized model

Let us introduce ∆t and ∆N the time and particle steps respectively. We set I t n := I(t∆t, n∆N ), for any t ∈ N and any n ∈ Z.

We have the choice between two Finite Difference schemes: either we deal with the discrete particle spacing defined as r t n := r(t∆t, n∆N ), for any (t, n) ∈ N × Z,

where r solves the Lagrangian GSOM model ( 9), or we consider the discrete particle position that reads X t n := X (t∆t, n∆N ), for any (t, n) ∈ N × Z, where (X , I) solves the Hamilton-Jacobi problem (10) associated to (9). Notice that the spatial extension of particle n ∈ Z is [X t n , X t n-1 [. In the first case, we would have to solve the following numerical scheme for ( 9)

       r t+1 n := r t n + ∆t ∆N V t n-1 -V t n , V t n := V r t n , I t n , I t+1 n = I t n + ∆tϕ I t n (11)
while in the second case, i.e. for (10) the appropriate numerical scheme is defined as follows

       X t+1 n = X t n + ∆tV t n , V t n := V X t n-1 -X t n ∆N , I t n , I t+1 n = I t n + ∆tϕ I t n ( 12
)
Notice that both approaches are very similar and give back the same results. Indeed, one can remark that if we have

r t n := X t n-1 -X t n ∆N , ( 13
)
then ( 11) is simply deduced from (12). By the way, knowing the spacing at each numerical steps, it is easy to compute the position of all particles, given the leader particle trajectory as a boundary condition.

It is worth noting that both schemes are first order schemes. The first one (11) can be interpreted as the seminal Godunov scheme (see [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF]) applied with demand and supply and the second discrete model ( 12) is an explicit Euler scheme.

In order to the numerical schemes ( 11) and ( 12) be monotone, time and label discrete steps need to satisfy a Courant-Friedrichs-Lewy (CFL) condition given by ∆N ∆t ≥ sup

N,r,t |∂ r V(r, I(t, N ))| . ( 14
)
For deducing a particle (Lagrangian) discretization of a traffic flow model on a junction, it is necessary to take into consideration different elements:

(i) the link model, which is given by either ( 11) or ( 12); (ii) the upstream (resp. downstream) boundary conditions for any incoming (resp. outgoing) link; (iii) the internal junction model, say the way particles are assigned from incoming road i ∈ I to outgoing road j ∈ J and eventually the internal dynamics of the junction point; (iv) link-junction and junction-link interfaces.

These constituting elements are addressed in what follows.

Downstream boundary condition

We assume that we are located at the downstream boundary of a given outgoing road j ∈ J . Assume the exit point S located at x S . The downstream boundary data at x S is given by the downstream supply σ which is discretized as σ t at time step t. Let n be the last particle located on the link (or at least a fraction η∆N of it is still on the link, with 0 < η ≤ 1). See Figure 2. We define the spacing associated to particle n as

r t n ∆N x S n -1 X t n r * r * r crit (I)
r t n := x S -X t n η∆N .
The fraction η is instantiated at the first time step t n following the exit of particle (n -1), as follows

η = x S -X tn n r tn n ∆N
. Now, we have to distinguish two cases:

• either V(r t n , I t n ) ≤ σ t r t n : in this case, the downstream supply is sufficient to accommodate the demand on the link. The spacing is conserved.

• or V(r t n , I t n ) > σ t r t n : in this case, the demand on the link cannot be fully satisfied since the downstream supply limits the outflow. Then, we have to solve

V(r t n , I t n ) = σ t r t n
and we choose the smallest value i.e. r t n = r * (see Figure 2). It means that we select the solution corresponding to the congested phase.

Then, we still update the position of particle n as usual, using (12). We also need to update the fraction η if the particle has not totally exited the link i.e. if X t+1 n < x S . The updated fraction is computed as follows

η ← η - ∆t r t n ∆N V r t n , I t n .

Upstream boundary conditions

We assume that we are located at the upstream boundary of a given incoming road i ∈ I. Consider the entry point E located at position x E . The boundary data is constituted by the upstream demand δ which is discretized as δ t for time step t. Let n be the last particle entered in the link. The next particle (n + 1) is still part of the demand. See Figure 3. Assume that particle (n + 1) will enter in the link at time (t + ε)∆t.

δ t n n + 1 r t n+1 ∆N r * r * r crit (I) V(., I) r x E X t n
Fig. 3. Illustration of upstream boundary condition.

We want to adapt to the upstream boundary condition, the methodology developed for the downstream boundary condition (which is detailed in Section 4.2) by introducing a proportion η of the particle that has already entered the link. The problem we face is that, unlike for the downstream where we know exactly the position of the last particle which has exited the link, we do not know precisely the position of the next particle which will enter the link. The situation at upstream is not exactly the inverse of what happens at the downstream that is why the algorithm is not so simple. Thus, we have to position the next particle that will enter the link.

More precisely, if one consider that the last particle that has entirely entered the link at time t is labeled n, we introduce a fraction η t of the particle (n + 1) which has already got into the link at time t. If we denote by q t the effective flow at the upstream entry and at time t then the proportion η t+1 at time (t + 1) is given by

η t+1 = q t (1 -ε)∆t
∆N where we assume that the particle n has entirely entered the link at time t n := (t + ε)∆t (see Figure 4). Then we have to compute the flow q t . If one consider a fictitious "junction" model just upstream the entry point, in which particles are stored before being injected into the link whenever it is possible, then we can deduce a stock model which is similar to an internal junction model. If F t is the number of particles stored inside the fictitious junction, then the evolution of the stock is given by

ε t+1 n+1 ∆t q t (t -1)∆t t n (n + 1) x E (n) x ∆t x E X t n ηr t n+1 ∆N n n + 1 σ t δ t t (t + 1)∆t t∆t t n+1
F t+1 = F t + (δ t -q t )∆t, (15)
where δ t is the (cumulative) demand and q t is the effective flow of particles which enters the link. Notice that the particle is generated at time (t + ε)∆t if and only if

F t + (δ t -q t )∆t = ∆N.
Then, with a simple test, we can distinguish two cases:

• if F t > 0, then there is a (vertical) queue just upstream the entry point and we get

q t = min σ t , Q max (I t n+1 ), F t ∆t + δ t ,
where Q max (I t n ) is the maximal flow obtained for the flow-density FD corresponding to the attribute I t n . • if F t = 0, then there is no queue and the flow is simply given by the minimum between the (local) upstream demand δ t which is given and the (local) downstream demand σ t , say

q t = min σ t , δ t .
We recall that the demand is defined according to (7),

say σ t = Ξ 1 r t n , I t n+1 , I t n ; x E .
In summary, the algorithm is composed as follows

(1) assume that we know the flow q t-1 passing through the entry point at time (t -1)∆t, (2) we update the fraction η t n+1 of particle (n + 1) which has already entered the link at time t such that η t n+1 = q t-1 ε t n ∆t ∆N where ε t n ∆t = t∆tt n and t n is the exact date at which the rear of particle (n) enters the link at x E .

(3) we also compute the spacing at time t according to particle (n + 1)

r t n+1 = X t n -x E η t n+1 ∆N
, and the exact position of particle (n + 1) and time t X t n+1 = X t nr t n+1 . (4) then we can compute the trajectory of particle (n+ 1) for following time steps as follows X t+1 n+1 = X t n+1 + ∆tV r t n+1 , I t n+1 , and we distinguish two cases:

• if X t+1 n+1 ≤ x E , then we go back to the first step and we itemize in time.

• if X t+1 n+1 > x E , then (the rear of) particle (n + 1) has entirely entered the link and we compute the exact time of its entry t n+1 as follows

t n+1 = t + (1 -ε t+1 n+1 ∆t, with ε t+1 n+1 = X t+1 n+1 -x E X t+1 n+1 -X t n+1 .
Then we itemize by considering next particle (n+ 2) (if it has been generated) and so on.

The algorithm in this case is more complex than in the case of "pure" Lagrangian described in [START_REF] Lebacque | A variational formulation for higher order macroscopic traffic flow models of the GSOM family[END_REF] but we can manage the exact arrival time of particles in the upstream buffer. Thus, the methodology can be directly applied to treat any junctionlink interface as we will see in what follows.

Internal state junction model

We consider a point-wise junction model with an internal state (first introduced in Khoshyaran and Lebacque ( 2009)) that is used as a buffer between incoming and outgoing branches of the junction. We recall that this buffer has internal dynamics and we can define an internal supply which depends on the number of stored particles.

In Eulerian framework, the internal state has some specific attributes such as    N z (t), total number of particles in the junction, N z,j (t), number of particles going on (j), I z (t), driver attribute in the junction.

Notice that the link-junction (resp. junction-link) interface is treated as a downstream (resp. upstream) boundary condition. Thus, we apply the algorithms described above, considering the local supply (resp. demand) of the buffers inside the junction point which are defined according to the number of stored particles.

There exists different strategies to deal with the assignment of particles through the junction. We can assume that we know

• either the assignment coefficients (α i,j ) i,j say the proportion of particles coming from any road i ∈ I that want to exit on road j ∈ J . Thus, a particle n ∈ Z entering from (i) may exit on road (j) with a probability α i,j ; • or the outgoing branch on which the particle n ∈ Z will exit; • or the origin-destination (OD) information (on a complex network) for any particle n ∈ Z. Coupled with an assignment model, one can compute the path of each particle.

Moreover we can distinguish (at least) two different cases for describing the internal dynamics of the junction. Indeed, one can consider that once particles have entered the junction, whatever are their origins, they are immediately assigned to the buffer corresponding to their wished outgoing branch j ∈ J . But it is also possible to consider that inside the junction point, any particle has a non-trivial travel time before to join their exit, which can be affected by the total number of particle inside the junction point or by the "physical" conflicts that can appear between the internal lanes of the junction point.

NUMERICAL EXAMPLE

Instantiation

Let us consider a junction with two incoming and two outgoing roads and the Colombo 1-phase model (see [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF]). It is noteworthy that with this choice, the speed-spacing FD V(•, •) is non-decreasing w.r.t. its second argument. The distribution of particle attribute I(•, •) is displayed on Figure 5. We consider initial and boundary conditions which correspond respectively to the initial positions of particles on the considered network, the upstream demands on each incoming link and the downstream supplies on the outgoing links. They are not displayed here by lack of available space. We assume that the junction point has a finite storage capacity.

Numerical result

The numerical solution to (10) is obtained thanks to (12). We then obtain the trajectories of each particle. For a better graphical representation, we consider the discrete traffic density ρ(•, •) defined as the inverse of the discrete spacing r t n for any t ∈ N, n ∈ Z (see also ( 13)). The density for the incoming link is plotted in Figure 6. The reader can notice that our numerical method can accurately recapture the shock wave due to the congestion and then the rarefaction wave, due to the decrease of the upstream demand, that mitigates the traffic jam later on. However, one can remark the numerical viscosity.

CONCLUSION AND FUTURE RESEARCH

In this paper, we have discussed a totally new numerical method to deal with the family of GSOM models posed on a junction. The generic GSOM model is recast in the Lagrangian framework and we have a careful look at the boundaries conditions for links and junctions. Notice that in our scheme vehicles are discretized into packets of ∆N particles. Hence, our scheme can be seen as a microscopic car-following model for the particular choice of ∆N = 1.

Recent models like [START_REF] Bressan | Continuous riemann solvers for traffic flow at a junction[END_REF] and [START_REF] Bressan | Conservation law models for traffic flow on a network of roads[END_REF] can be fully recast into the framework described in our article and solved using our algorithm. Indeed, the attribute is given by the assignment coefficients which are hopefully advected with the traffic flow (if users do not change their minds).

By the way, we highlight below some interesting research directions. The discrete model ( 12) can be replaced by an implicit scheme or more complex time integration schemes (see for instance [START_REF] Treiber | Comparing numerical integration schemes for time-continuous car-following models[END_REF]). Such numerical schemes can be justified mainly if we consider a source term at the r.h.s. in (9, say ϕ(I) = 0. In the particular case of ϕ = 0, explicit Euler scheme is very satisfying.

Another direction of research would be to numerically compare our method and the variational approach (see [START_REF] Costeseque | A variational formulation for higher order macroscopic traffic flow models: numerical investigation[END_REF]) adapted for junction modeling, which has not been done right now.

Fig. 2 .

 2 Fig. 2. Illustration of downstream boundary condition.

Fig. 4 .

 4 Fig. 4. Illustration of new upstream boundary condition.

  Fig. 5. Particle attribute values I(•, •).

Fig. 6 .

 6 Fig. 6. Numerical density values ρ(•, •).