Oliver Kolb
email: kolb@uni-mannheim.de

Simone Göttlich
email: goettlich@uni-mannheim.de

Paola Goatin
email: paola.goatin@inria.fr

Capacity drop and traffic control for a second order traffic model

Keywords: AMS Classification. 90B20, 35L65 Keywords. Traffic flow, second order model, on-ramp coupling, numerical simulations, optimal control

In this paper, we illustrate how second order traffic flow models, in our case the Aw-Rascle equations, can be used to reproduce empirical observations such as the capacity drop at merges and solve related optimal control problems. To this aim, we propose a model for on-ramp junctions and derive suitable coupling conditions. These are associated to the first order Godunov scheme to numerically study the well-known capacity drop effect, where the outflow of the system is significantly below the expected maximum. Control issues such as speed and ramp meter control are also addressed in a first discretize-then optimize framework.

Introduction

The description of vehicular traffic flow based on systems of conservation or balance laws has been proposed by many researchers during the last decades, see for example [START_REF] Garavello | Traffic flow on networks[END_REF][START_REF] Treiber | Data, models and simulation[END_REF] and the references therein. In contrast to first order traffic models consisting only of one scalar equation for the traffic density, second order models are characterized by a second equation for the evolution of traffic mean velocity. Among these, Aw-Rascle [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] proposed a model that overcomes the drawbacks of previous models, and is capable to capture some traffic characteristics linked to the anticipation behaviour of drivers. Greenberg [START_REF] Greenberg | Extensions and amplifications of a traffic model of Aw and Rascle[END_REF] presented an extension of the Aw-Rascle model including a relaxation term towards a preferred (or equilibrium) velocity. Since then many contributions have appeared on theoretical and numerical investigations [START_REF] Aw | Derivation of continuum traffic flow models from microscopic follow-the-leader models[END_REF][START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF][START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF][START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF], networked problems [START_REF] Garavello | Traffic flow on a road network using the Aw-Rascle model[END_REF][START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF][START_REF] Herty | Optimization criteria for modelling intersections of vehicular traffic flow[END_REF][START_REF] Herty | Coupling conditions for a class of second-order models for traffic flow[END_REF][START_REF] Siebel | Balanced vehicular traffic at a bottleneck[END_REF], and more recently phenomenological and data driven approaches [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF][START_REF] Jin | A kinematic wave theory of capacity drop[END_REF][START_REF] Leclercq | Capacity drops at merges: New analytical investigations[END_REF][START_REF] Piccoli | Second-order models and traffic data from mobile sensors[END_REF][START_REF] Srivastava | Empirical observations of capacity drop in freeway merges with ramp control and integration in a first-order model[END_REF].

In this work, we present results concerning the coupling of interlinked road networks to onramps. It is known that in such situations capacity drops occur for high traffic densities, see e.g. [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF][START_REF] Jin | A kinematic wave theory of capacity drop[END_REF][START_REF] Leclercq | Capacity drops at merges: New analytical investigations[END_REF][START_REF] Parzani | Second-order model and capacity drop at merge[END_REF][START_REF] Srivastava | Empirical observations of capacity drop in freeway merges with ramp control and integration in a first-order model[END_REF]. Capacity drop means that the measured outflow of the system is smaller than the expected maximal flow if the system is well utilized. The difference might be up to 10% or even higher [START_REF] Jin | A kinematic wave theory of capacity drop[END_REF][START_REF] Srivastava | Empirical observations of capacity drop in freeway merges with ramp control and integration in a first-order model[END_REF]. This can be explained by a less efficient driving style when exiting a congested zone (upstream the on-ramp), which reduces the downstream flow compared to free-flow conditions [START_REF] Treiber | Data, models and simulation[END_REF]. The coupling of on-ramps to roads is similar to the modeling of a merge [START_REF] Jin | On the distribution schemes for determining flows through a merge[END_REF]. However, in contrast to the majority of already existing approaches, here we assume that the on-ramp is given by an ordinary differential equation, see [START_REF] Delle Monache | A PDE-ODE model for a junction with ramp buffer[END_REF][START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF]. For second order models, the challenge here is to find appropriate conditions to ensure the conservation of mass and momentum flow. Once these coupling conditions are defined, they can be integrated in a finite-volume type numerical scheme to compute the evolution of traffic conditions on the network. We use a first-order Godunov scheme combined with a time splitting to resolve the relaxation term. In particular, we observe that the model is able to reproduce the capacity drop effect for increasing inflows, similar to the study made in [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF] and in contrast to [START_REF] Parzani | Second-order model and capacity drop at merge[END_REF], where the phenomenon is only observed for a short time period. Furthermore, our approach gives the opportunity to consider questions of optimal control in a first discretize-then optimize framework. Typical control issues arising in the on-ramp context are speed limit and ramp metering. These kind of traffic flow control problems have been mainly studied in the context of first-order models based on the Lighthill-Whitham-Richards (LWR) equations [START_REF] Lighthill | On Kinematic Waves. II. A Theory of Traffic Flow on Long Crowded Roads[END_REF], see e.g. [START_REF] Delle Monache | Traffic regulation via controlled speed limit[END_REF][START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF][START_REF] Hegyi | Optimal coordination of variable speed limits to suppress shock waves[END_REF][START_REF] Reilly | Adjoint-based optimization on a network of discretized scalar conservation laws with applications to coordinated ramp metering[END_REF]. To the best of our knowledge, speed limit and ramp meter control has been not considered for the Aw-Rascle model so far. We discuss the numerical results obtained from the adjoint calculus and compare them to the LWR equations.

The outline is as follows: In Section 2 we present the Aw-Rascle model with relaxation term and the coupling conditions at on-ramps. Our discretization scheme will be introduced in Section 3. Numerical results concerning the capacity drop and optimal control issues can be found in Section 4. Comparisons to the solution of the LWR model are also given for all experiments.

Modeling

We introduce the modeling equations given by the Aw-Rascle (AR) model [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] including a relaxation term as originally introduced in [START_REF] Greenberg | Extensions and amplifications of a traffic model of Aw and Rascle[END_REF]. The model consists of a 2 × 2 system of conservation laws for the density and a sort of generalized momentum derived from the (anticipate) acceleration equation.

In this work, our focus is on the coupling of the Aw-Rascle model to on-ramps whose dynamics are modeled by ordinary differential equations. We will derive appropriate coupling conditions that allow for the definition of boundary conditions at junctions. As we will see in the numerical experiments, our modeling approach covers the capacity drop effect and allows to solve optimal control problems such as speed limit and ramp metering, cf. subsection 4.3.

Road dynamics

We briefly recall the Aw-Rascle traffic flow model [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] and explain how it can be extended to the context of networks. Traffic states are described by the density ρ i (x, t) and the mean speed of vehicles v i (x, t) on each road i at position x and time t.

Given some initial state ρ i (x, 0), v i (x, 0) on each road i, the dynamics for x ∈ (0, L i ) and t ∈ (0, T ) are described by [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF][START_REF] Greenberg | Extensions and amplifications of a traffic model of Aw and Rascle[END_REF]

∂ t ρ i + ∂ x (ρ i v i ) = 0, (1)
∂ t (ρ i w i ) + ∂ x (ρ i v i w i ) = -ρ i v i -V i (ρ i ) δ , (2)
w i = v i + p i (ρ i ), (3)
or in conservative form with

y i = ρ i w i ∂ t ρ i y i + ∂ x y i -ρ i p i (ρ i ) y i -ρ i p i (ρ i ) yi ρi = 0 -(yi-ρipi(ρi))-ρiVi(ρi) δ =gi(ρi,yi) , (4)
where p i (ρ) is a known pressure function satisfying p i (ρ) > 0 and ρp i (ρ) + 2p i (ρ) > 0 for all ρ.

The latter condition ensures that the curve {w i (ρ, v) = v + p i (ρ) ! = c} for any constant c > 0 is strictly concave and passes through the origin of the (ρ, ρv)-plane. Further, there exists a sonic point σ i (c) maximizing the flux ρv along the curve {w i (ρ, v) = c}. The relaxation term in the velocity equation includes the preferential speed function V i (ρ) which the drivers tend to adopt and a relaxation parameter δ.

We consider the preferential velocity depending on the density

V i (ρ) = v max i 1 - ρ ρ max i (5)
and the pressure function (as in [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF][START_REF] Parzani | Second-order model and capacity drop at merge[END_REF]) [START_REF] Delle Monache | A PDE-ODE model for a junction with ramp buffer[END_REF] equipped with maximal density ρ max i > 0, maximal velocity v max i > 0, reference velocity v ref i > 0 and γ i > 0. Later, in Section 4.3, we will also consider time dependent velocities v max i (t) and v ref i (t) so that actually p i (ρ) = p i (ρ, t), V i (ρ) = V i (ρ, t) and g i (ρ, y) = g i (ρ, y, t).

p i (ρ) = v ref i γ i ρ ρ max i γi
Note that above and in the following, we use w i = w i (x, t) as space and time dependent state variable but also as function, e.g. in the form w i (ρ, v, t) = v + p i (ρ, t). Analogously, we will use the notation v i (U, t) = wp i (ρ, t) for the velocity of a state U = (ρ, ρw). Similar to first order traffic models, we define the demand and supply functions for each road i as follows: For a given constant c (corresponding to a fixed value of w) we have

D i (ρ, c, t) = c -p i (ρ, t) ρ if ρ ≤ σ i (c, t), c -p i (σ i (c, t), t) σ i (c, t) if ρ ≥ σ i (c, t), (7)
S i (ρ, c, t) = c -p i (σ i (c, t), t) σ i (c, t) if ρ ≤ σ i (c, t), c -p i (ρ, t) ρ if ρ ≥ σ i (c, t), (8)
where

σ i (c, t) = ρ max i c γ i v ref i (t) (1 + γ i ) 1 γ i (9)
is the sonic point on the curve {w i (ρ, v, t) = v + p i (ρ, t) ! = c} in the (ρ, ρv)-plane. An illustration of the considered demand and supply functions is given in Figure 1. Supply and demand functions are needed to formulate the coupling of different roads or on-ramps.

Coupling and boundary conditions

In the following, we present the serial coupling between single roads as well as the coupling with an on-ramp. The considered coupling and boundary conditions can be given in terms of mass flow q = ρv and "momentum flow" qw. The computation of the actual states at a junction is not necessary (cf. the discretization in Section 3). Different to other coupling conditions [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF][START_REF] Herty | Optimization criteria for modelling intersections of vehicular traffic flow[END_REF][START_REF] Jin | On the distribution schemes for determining flows through a merge[END_REF][START_REF] Parzani | Second-order model and capacity drop at merge[END_REF][START_REF] Siebel | Balanced vehicular traffic at a bottleneck[END_REF], we assume that in the situation of a merge the priority parameter P is independent of the demand of the ingoing roads. This guarantees the consistency of the corresponding Riemann solver, see [9, Definition 4.2.2].

1-to-1 junction

q 1 q 2 Figure 2: 1-to-1 junction.
We use index 1 for the ingoing road, index 2 for the outgoing road, and consider the data U i = (ρ i , ρ i w i ) at the adjacent boundaries of roads 1 and 2, respectively. Note that we omit the time dependency of the states U i and the fluxes q i in the following for a better reading, whereas we still indicate the explicit time dependency of the demand and supply functions.

Flow maximization at the junction over all admissible states leads to

q 1 = q 2 = q = min D 1 ρ 1 , w 1 , t , S 2 ρ2 , w 1 , t , (10)
where ρ2 is either obtained by the intersection of the curves

{v 2 (U, t) = v 2 (U 2 , t)} and {w 2 (U, t) = w 1 }, (11)
or ρ2 = 0. Note that for the considered pressure functions this problem can be solved explicitly:

ρ2 = ρ max 2 max γ 2 v ref 2 (t) w 1 -v 2 (U 2 , t) , 0 1 γ 2 . ( 12
)
Then, q 1 = q 2 = q determines the mass flow out of road 1/into road 2. Further, with w = w 1 at the junction, the momentum flow is given by q w.

On-ramp at a 1-to-1 junction We consider a 1-to-1 junction with an on-ramp. Here, we use index 1 for the ingoing road, index 2 for the on-ramp, index 3 for the outgoing road, and consider the data U i = (ρ i , ρ i w i ) on roads 1 and 3, respectively, and the queue length l 2 at the on-ramp.

q 1 q 3 q 2
For the demand D 2 (l, t) at the on-ramp we apply

D 2 (l, t) = u 2 (t) f max 2 if l > 0, min{f in 2 (t), f max 2 } if l = 0, ( 13
)
where l is the length of the queue, u 2 (t) ∈ [0, 1] is the (time-dependent) metering rate, f in 2 (t) the "inflow" of cars arriving at the on-ramp and f max 2 the maximum flow onto the main road. To get a unique solution, we assign the priority parameter P to road 1, and apply

q 1 = min D 1 ρ 1 , w 1 , t , max P S 3 ρ3 , w 1 , t , S 3 ρ3 , w 1 , t -D 2 l 2 , t , (14)
q 2 = min D 2 l 2 , t , max (1 -P )S 3 ρ3 , w 1 , t , S 3 ρ3 , w 1 , t -D 1 ρ 1 , w 1 , t , (15)
q 3 = q 1 + q 2 , ( 16
)
where ρ3 is either obtained by the intersection of the curves

{v 3 (U, t) = v 3 (U 3 , t)} and {w 3 (U, t) = w 1 }, (17)
or ρ3 = 0. Again, the density ρ3 can be computed explicitly:

ρ3 = ρ max 3 max γ 3 v ref 3 (t) w 1 -v 3 (U 3 , t) , 0 1 γ 3 . ( 18
)
As before, the momentum flow is computed by multiplication of q 1 and q 3 with w = w 1 , while the formal momentum flow wq 2 from the on-ramp does not influence the state of the on-ramp and is therefore not computed. The evolution of the queue at the on-ramp is given by

d l 2 (t) dt = f in 2 (t) -q 2 , (19)
where we typically start with empty queues, i.e., l 2 (0) = 0.

On-ramp at origin

q 1 q 2
Figure 4: On-ramp at origin.

We consider an on-ramp at a vertex with only one outgoing road. We use index 1 for the on-ramp with queue length l 1 , and index 2 for the road with data U 2 = (ρ 2 , ρ 2 w 2 ).

Similar to above we apply

D 1 (l, t) = u 1 (t) f max 1 if l > 0, min{f in 1 (t), f max 1 } if l = 0 ( 20
)
for the demand at the on-ramp, where l is the length of the queue, u 1 (t) ∈ [0, 1] is the (timedependent) metering rate, f in 1 (t) the inflow of cars arriving at the on-ramp and f max 1 the maximum flow onto the road.

In this case there is no value for w to evaluate the supply function of the road, which is needed to determine the actual flux onto the road q in a similar fashion as in the case of a simple 1-to-1 junction, i.e., similar to [START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF]. Therefore, we consider an auxiliary left state U 1 mimicking the desired inflow of the on-ramp. We assume that the velocity of the auxiliary state is at equilibrium with respect to the desired velocity of the road, i.e., we search for a state

U 1 = ρ, ρ w such that ρV 2 (ρ, t) = D 1 (l, t), (21)
where w = w 2 (ρ, V 2 (ρ, t), t). Assuming that the flux D 1 (l, t) ≤ f max 1 can be realized with the velocity function V 2 , we get the solutions

ρ ± = ρ max 2 2 ± (ρ max 2 ) 2 4 - ρ max 2 D 1 (l, t) v max 2 (t) . ( 22
)
Here we choose ρ -, since this choice also fulfils

D 2 (ρ -, w, t) = D 1 (l, t). (23)
Finally, we apply

q 1 = q 2 = q = min D 1 (l 1 , t), S 2 ρ2 , w, t (24)
with w = w 2 (ρ -, V 2 (ρ -, t), t), and ρ2 is either obtained by the intersection of the curves

{v 2 (U, t) = v 2 (U 2 , t)} and {w 2 (U, t) = w}, (25)
or ρ2 = 0. The explicit representation is the same as [START_REF] Greenberg | Extensions and amplifications of a traffic model of Aw and Rascle[END_REF]. Similar to above, the evolution of the queue at the on-ramp is given by

d l 1 (t) dt = f in 1 (t) -q 1 . (26)
Outflow conditions

q 1 q 1 Figure 5: Outflow at a vertex.
At nodes with only one ingoing road (with index 1) we consider absorbing boundary conditions up to a given maximum flow rate f out 1 (t):

q 1 = min D 1 ρ 1 , w 1 , t , f out 1 (t) . ( 27
)
The momentum flow is given by q 1 w 1 .

Discretization

For the numerical solution of the described model, a finite number of time points t n = n∆t is considered, where ∆t = T /Nt. Further, each road i is divided into Nx i cells of equal size ∆x i = L i /Nx i . To fulfill the CFL condition, we will claim (at least) ∆t

v max i ≤ ∆x i in all numerical examples.

Roads

The underlying balance law (4) is discretized with a fractional step method: We use a first order Godunov scheme for the flux term and an implicit Euler method for the relaxation term. Due to the fundamental idea of the Godunov scheme, the computation of the actual states at junctions is not necessary and the coupling will be given in terms of fluxes only. Note that in the following U i = (ρ i , ρ i w i ) denotes again the state.

On each road i the initial conditions (for j ∈ {1, . . . , Nx i }) are given by the cell averages

U 0 i,j-0.5 = 1 ∆x i j∆xi (j-1)∆xi U i (x, 0) dx . (28)
Then, for n ∈ {0, . . . , Nt -1} and j ∈ {1, . . . , Nx i } we apply

Ũ n+1 i,j-0.5 = U n i,j-0.5 - ∆t ∆x i F n i,j -F n i,j-1 , (transport) (29)
U n+1 i,j-0.5 = Ũ n+1 i,j-0.5 + ∆t g(U n+1 i,j-0.5 ).

(relaxation term)

The flux terms F n i,0 and F n i,Nxi are given via coupling/boundary conditions (see next subsection) and for j ∈ {1, . . . , Nx i -1} we have

F n i,j = q n i,j
w n i,j-0.5 q n i,j with q n i,j = min D i (ρ n i,j-0.5 , w n i,j-0.5 , t n ), S i (ρ n i,j , w n i,j-0.5 , t n ) , (31) where ρn i,j is either obtained by the intersection of the curves

{v i (U, t n ) = v n i,j+0.5 } and {w i (U, t n ) = w n i,j-0.5 }, (32)
or ρn i,j = 0. Note that this treatment is consistent with the treatment of a 1-to-1 junction where two roads with the same parameters meet.

Coupling and boundary conditions

We complete the numerical scheme by the description of the coupling and boundary conditions according to the model equations from Section 2.2. As already mentioned, we only need the coupling fluxes here and not the actual states at the junctions.

1-to-1 junction

As before, we use index 1 for the ingoing road and index 2 for the outgoing road. According to [START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF] we apply for n ∈ {0, . . . , Nt -1} qn = min D 1 ρ n 1,Nx1-0.5 , w n 1,Nx1-0.5 , t n , S 2 ρn 2,0 , w n 1,Nx1-0.5 , t n ,

where ρn 2,0 is either obtained by the intersection of the curves

{v 2 (U, t n ) = v n 2,0.5 } and {w 2 (U, t n ) = w n 1,Nx1-0.5 }, (34)
or ρn 2,0 = 0. Due to the choice of the pressure function (6), ρn 2,0 can be directly computed according to [START_REF] Greenberg | Extensions and amplifications of a traffic model of Aw and Rascle[END_REF]. Finally, we get

F n 1,Nx1 = F n 2,0 = qn w n 1,Nxi-0.5 qn . ( 35
)
On-ramp at a 1-to-1 junction

We use index 1 for the ingoing road, index 2 for the on-ramp and index 3 for the outgoing road.

According to [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF] we apply

D2 (l, t) = u 2 (t) min f in 2 (t) + l ∆t , f max 2 (36)
for the demand at the on-ramp. Then, using the abbreviations

d n 1 = D 1 ρ n 1,Nx1-0.5 , w n 1,Nx1-0.5 , t n , ( 37
)
d n 2 = D2 (l n 2 , t n ), (38)
s n 3 = S 3 ρn 3,0 , w n 1,Nx1-0.5 , t n , (39)
we get from ( 14) and ( 15)

q n 1 = min {d 1 , max {P s 3 , s 3 -d 2 }} , (40)
q n 2 = min {d 2 , max {(1 -P )s 3 , s 3 -d 1 }} , (41)
q n 3 = q n 1 + q n 2 , (42)
where ρn 3,0 is either obtained by the intersection of the curves

{v 3 (U, t n ) = v n 3,0.5 } and {w 3 (U, t n ) = w n 1,Nx1-0.5 }, (43)
or ρn 3,0 = 0. Again, in our case a direct computation is possible due to [START_REF] Jin | A kinematic wave theory of capacity drop[END_REF]. Finally, we have

F n 1,Nx1 = q n 1 w n 1,Nxi-0.5 q n 1 , F n 3,0 = q n 3 w n 1,Nxi-0.5 q n 3 ( 44
)
and apply

l n+1 2 = l n 2 + ∆t f in 2 (t n ) -q n 2 (45)
for the update of the queue length via the explicit Euler scheme. Note that we always achieve a non-negative queue length by our choice of the demand function (36). Alternatively, the time step ∆t could be adapted, cf. [START_REF] Delle Monache | A PDE-ODE model for a junction with ramp buffer[END_REF].

On-ramp at origin

We use index 1 for the on-ramp and index 2 for the road. Similar to above we apply

D1 (l, t) = u 1 (t) min f in 1 (t) + l ∆t , f max 1 ( 46
)
for the demand at the on-ramp. Next, following the idea introduced in Section 2.2, we compute

ρ n 1,-= ρ max 2 2 - (ρ max 2 ) 2 4 - ρ max 2 D1 (l n 1 , t n ) v max 2 (t n ) ( 47
)
and qn = min D1 (l n 1 , t n ), S 2 ρn 2,0 , wn , t n (48)
with wn = w 2 (ρ n 1,-, V 2 (ρ n 1,-, t n ), t n ), and ρn 2,0 is either obtained by the intersection of the curves

{v 2 (U, t) = v n 2,0.5 } and {w 2 (U, t n ) = wn }, (49)
or ρn 2,0 = 0. Again, ρn 2,0 can be computed explicitly according to [START_REF] Greenberg | Extensions and amplifications of a traffic model of Aw and Rascle[END_REF]. Finally, we use

F n 2,0 = qn wn qn ( 50
)
and apply

l n+1 1 = l n 1 + ∆t f in 1 (t n ) -qn (51)
for the update of the queue length. Again, non-negativity of the queue length is ensured by our choice of the demand function (46).

Outflow conditions

At nodes with only one ingoing road (with index 1) we apply according to ( 27)

F n 1,Nx1 = q n 1 w n 1,Nx1-0.5 q n 1 (52) with q n 1 = min D 1 ρ n 1,Nx1-0.5 , w n 1,Nx1-0.5 , t n , f out 1 (t n ) . (53)
In the next section, we present a numerical study to show the characteristics of our modeling approach. We also comment on already existing examples from the literature, pointing out differences and similarities.

Numerical results

Our numerical results emphasize on three different scenarios that our approach is able to deal with from a modeling and computational point of view. First, we show that, for a single road, the Aw-Rascle (AR) model with relaxation term numerically converges to the Lighthill-Whitham-Richards (LWR) model as expected [START_REF] Greenberg | Extensions and amplifications of a traffic model of Aw and Rascle[END_REF]. The second example considers the capacity drop for a 1-to-1 junction with on-ramp. We will see that this phenomenon is not covered by the classical LWR model. The third example tackles speed control and coordinated ramp metering control issues. It will turn out that the control strategies considered individually do not lead to the best possible travel times, while the combination of both lead to satisfying results.

Aw-Rascle towards LWR

We consider a single road as depicted in Figure 6 with length L = 1 km, ρ max = 200 cars km , v max = v ref = 100 km h and ρ 0 = 80 cars km as initial density. For the AR model we further use γ = 2, v 0 = 60 km h (in equilibrium state) and decreasing values for the relaxation parameter δ starting from a usual value of δ = 0.005 hours (18 seconds, see e.g. [START_REF] Hegyi | Optimal coordination of variable speed limits to suppress shock waves[END_REF]). We consider a simulation time of T = 36 seconds and f in (t) = 4000 + 800 cos(4π t T ) for 0 ≤ t ≤ T 2 , 4800 otherwise as inflow profile at the origin "in" (in cars h ). As discretization parameters, we apply ∆x = 10 meters and ∆t = 0.18 seconds Figure 7 shows density and velocity profiles on the whole road after 36 seconds. As expected, the solution of the Aw-Rascle model with relaxation term tends towards the LWR model for decreasing relaxation parameter δ. Moreover, decreasing densities are correlated with increasing velocities. Note that the application of larger values for the parameter γ in the pressure term leads to higher peaks in the simulation results (not shown here).

Capacity drop

Now we turn our attention to study the capacity drop effect, cf. [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF]. We will see that increasing the inflow into the network depicted in Figure 8 will at a certain point lead to a decreasing outflow. Thus, a significant difference between the effective and maximal possible outflow occurs in the Aw-Rascle model. Conversely, this effect is not captured by the classical LWR model. This is due to the different modeling of the supply function. In the congested regime, the value of the supply function is at its maximum in the LWR model and therefore no capacity drop occurs. However, in the AR model, the supply function additionally depends on the "incoming" value of w that allows for a different behavior.

in on-ramp out road1 road2

Figure 8: Two roads with an on-ramp in between.

We consider two roads of 1 km length with an on-ramp in between, see Figure 8. On both roads we use ρ max = 180 cars km , v max = v ref = 100 km h , γ = 2, δ = 0.005 h and an initial density of 50 cars km . At the origin "in" we consider a constant (desired) inflow f in 1 = 3500 cars h . At the on-ramp we consider different (desired) inflows, starting from f in 2 = 500 cars h up to 2500 cars h and down to 500 cars h again. As priority parameter at the on-ramp we use P = 0.5. Table 1 and Figure 9 show the simulation results for the AR and the LWR model (with ∆x = 100 metres and ∆t = 1.8 seconds). The first two columns report the desired and the actual inflow at the on-ramp. The following three values are the resulting density, velocity and the value of w just upstream the on-ramp (end of the first road). The last two columns show the total outflow at the end of the second road. We begin with the interpretation of the results obtained from the AR model: Obviously, up to an inflow of 1000 cars h at the on-ramp, the total inflow (of 4500 cars h ) is within the capacity of the given road. When the inflow at the on-ramp further increases, the resulting total flux cannot be realized at some point. From then on, the value of w at the end of the first road directly influences the total flux entering the second road (compare equations ( 14) and [START_REF] Herty | Optimization criteria for modelling intersections of vehicular traffic flow[END_REF] and the definition of the supply function ( 8)). As a consequence, the outflow at the end of the second road for the cases f in 2 ∈ {1500 cars h , 2000 cars h , 2500 cars h } is lower than the outflow for the cases f in 2 ∈ {500 cars h , 1000 cars h }. Due to the choice of the priority parameter P = 0.5, the effect of a decreasing outflow while the desired inflow increases stagnates as soon as the fluxes from the first road and the on-ramp are equal, which is the case for f in 2 ∈ {2000 cars h , 2500 cars h }. Interestingly, even when the desired inflow at the on-ramp is lowered again, the original state for the same values of f in 1 and f in 2 is not reached. This is due to the fact that the outflow in the congested situation is below the accumulated desired inflow. Accordingly, the queue at the origin permanently increases in the final situation even though the capacity of the road could handle the desired inflows in the free flow situation.

Note that compared to the existing literature the capacity drop is investigated for different situations. In [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF] no relaxation term is considered but a similar result is obtained. The construction principle for the capacity drop is the same (decreasing value of w for increasing flux), but in our investigations the possible equilibrium states are restricted due to the relaxation term so that higher densities on the incoming road lead to lower velocities (in difference to the equilibria 3 and 5 in [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF]Table 3]). In [START_REF] Parzani | Second-order model and capacity drop at merge[END_REF], where a different junction model and also a different numerical scheme is used, the capacity drop effect for the model with relaxation term is only observed for short time intervals, whereas we identify a permanent drop, cf. Figure 9. Finally, considering the same scenario with the LWR model (see [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF] for the details), one observes that the outflow increases with the accumulated desired inflows until it reaches the maximum capacity of 4500 cars h = ρ max 2 • v max 2 (see again Table 1 and Figure 9). According to the chosen priority P = 0.5, the actual inflows at the origin and at the on-ramp are 2250 cars h in that situation. Any further increase of the desired inflow at the on-ramp does not have any effect on the situation on the roads. Further, unlike the AR model, a decrease of the desired inflow f in 2 below 1000 cars h leads back to the original situation (as soon as the queues have emptied).

Coordinated speed control and ramp metering

For optimization purposes, we are finally interested in controlling maximal speed limits as well as on-ramp inflows to minimize the total travel time. As already introduced in equations ( 13) and [START_REF] Kolb | Mathematical optimization of water networks[END_REF], the ramp metering rate u i (t) is used to control the demand at the on-ramp. Similar to our recent work [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF], we will use

v max i = v max i (t) ∈ [v low i , v high i
] as another time-dependent (piecewise constant) control variable so that actually V i (ρ) = V i (ρ, t) and g i (ρ, y) = g i (ρ, y, t) for the relaxation term. Regarding the influence of the pressure term, we will consider two variants:

1. v ref i (t) = v max i (t), i.e.
, the pressure term directly depends on the (controllable) speed limit and therefore p i (ρ) = p i (ρ, t);

2. v ref i (t) = v high i , i.e.
, the pressure term is independent of the current speed limit in the sense that ∂pi ∂v max i = 0.

For our investigations, we consider the road network depicted in Figure 10 (adapted from [START_REF] Hegyi | Optimal coordination of variable speed limits to suppress shock waves[END_REF]) with the road parameters of Table 2. The exponents in the pressure function are γ i = γ = 2 for all roads and δ = 0.005 h for the relaxation parameter. The priority parameter P at the on-ramp equals 0.5 and f max = 2000 cars h . At the origin "in" we consider f max = 4000 cars h .

in mid mid2 on-ramp out road1 road2 road3 road4

Figure 10: Road network with an on-ramp at the node "on-ramp". We consider a time horizon of T = 3.0 hours and the boundary conditions shown in Figure 11. For the given scenario, we are interested in minimizing the total travel time

roads i T 0 Li 0 ρ i (x, t) dx dt + on-ramps j T 0 l j (t) dt , (54)
given an upper bound of 100 vehicles in the queue of the on-ramp. To solve this optimization task, we apply a first-discretize-then-optimize approach and adjoint calculus, see also [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF]. Thus, for given control decisions (speed limits and metering rates), the discretization scheme described in Section 3 is always used to evaluate the objective function (54) (using the trapezoidal rule for quadrature). Further, sensitivity/gradient information with respect to the control decisions is computed based on the same discretization (we refer to [START_REF] Kolb | Simulation and Optimization of Gas and Water Supply Networks[END_REF][START_REF] Kolb | Mathematical optimization of water networks[END_REF] for more details). Finally, the SQP solver DONLP2 [START_REF] Spellucci | A new technique for inconsistent QP problems in the SQP method[END_REF][START_REF]An SQP method for general nonlinear programs using only equality constrained subproblems[END_REF] is used for the optimization of the control decisions. The applied discretization parameters are ∆x = 250 meters and ∆t = 7.2 seconds. Table 3 shows the total travel times for the different models with and without optimization. The resulting queues for the AR model with scaling of the pressure function ( ∂pi ∂v max i = 0) are shown in Figure 12. Without ramp metering and speed control the queue at the on-ramp stays empty whereas more than 300 cars accumulate in the queue at the origin. When only ramp metering is considered, the queue at the origin is reduced to zero, while up to 100 cars accumulate in the queue at the on-ramp. Obviously, the prescribed upper bound for the queue at the on-ramp is active, impeding further improvement of the objective function. When ramp metering and speed control are used, the upper bound of 100 cars is even never touched. In the case of speed control only, the queue at the on-ramp stays empty, whereas the queue at the origin cannot be kept empty during the complete time horizon.

Note that the results for the model without scaling of the pressure function ( ∂pi ∂v max i = 0) are quite similar considering the queues and therefore are not plotted. Nevertheless the possible improvement with respect to total travel time is much larger in the case of speed control only (872.6 versus 1122.8). Concerning the LWR model, both queues stay empty during the whole time horizon already in the uncontrolled case and the optimization procedure terminates directly in the other cases -leading to the same result. The optimality of the uncontrolled case with respect to total travel time results from the fact that the classical LWR model does not capture the capacity drop effect (cf. Section 4.2). Figure 13 exemplarily shows where the improvement in total travel time comes from in the case of coordinated ramp metering and speed control in the model with scaling of the pressure function: The outflow of the system (plot on the right) in the optimized scenario is permanently above the outflow of the uncontrolled system until it drops to the low inflow level after two hours. Figure 14 shows the applied (optimal) controls. As shown in Table 3, the uncontrolled case is already optimal for the LWR model, i.e., no ramp metering or speed control is required. Accordingly, we take a look at the influence of the relaxation parameter δ within the second order model, which regulates the acceleration towards the desired velocity V (ρ, t) of the LWR model.

Table 4 shows the simulation and optimization results for a wide choice of δ (on all roads). Since smaller values of δ lead to larger values for the acceleration towards the desired velocity, one would expect that the results for small values of δ get close to the results of the LWR model (as in Section 4.1). Surprisingly, it is the other way around for the considered scenario. 725.9 725.9 725.9

First, one observes larger total travel times in the case of smaller values for δ. These result from the faster reaction of the velocity on the increase in density behind the on-ramp at the beginning of the scenario, which in turn accelerates the increase in density compared to the other parameter choices (see Figure 15 where the results for the case ∂pi ∂v max i = 0 are plotted). For very large values of δ there is almost no reaction of the velocity on the density increase and no congestion. Qualitatively, the corresponding results are quite similar to the results of the LWR model (see again Figure 15).

Conclusion

In this work, we have set up coupling conditions for the Aw-Rascle model linked to on-ramps described by ordinary differential equations. Applying a first order Godunov scheme, the presented approach allows for numerical investigations of the capacity drop effect and optimal control problems.

Future work will include extensions to on-ramps that are described by first order traffic models such as the LWR model. We will derive appropriate coupling conditions and focus on theoretical properties.

Figure 1 :

 1 Figure 1: Demand and supply functions on a fixed road for ρ max = 1, v ref = 2, γ = 2 and the given values for c.

Figure 3 : 1 -

 31 Figure 3: 1-to-1 junction with on-ramp.

Figure 6 :

 6 Figure 6: A single road.

δ = 5 • 5 • 10 - 5 Figure 7 :

 551057 Figure 7: Density (top) and velocity (bottom) after 36 seconds for different choices of the relaxation parameter δ and for the LWR model.

Figure 9 :

 9 Figure 9: Actual outflow depending on the desired inflow at the on-ramp for the AR and the LWR model.

Figure 11 :

 11 Figure 11: Inflow profiles for the network in Figure 10.

  queue at on-ramp[cars]

Figure 12 :

 12 Figure 12: Queue at the origin "in" and the on-ramp with and without optimization.

Figure 13 :

 13 Figure 13: Flow at the origin "in" of the network (left) and at the node "out" (right) with and without optimization.

Figure 14 :

 14 Figure 14: Optimal control of v max i (t) on road2 (top left) and road3 (top right) and u(t) at the on-ramp (bottom).

δ = 5 • 5 • 10 - 1

 55101 10 -3 (original) δ = 5 • 10 -5 (smallest) δ =

Figure 15 :

 15 Figure 15: Density (top) and velocity (bottom) behind the on-ramp in the uncontrolled case for different choices of the relaxation parameter δ and for the LWR model.

Table 1 :

 1

	Capacity drop effect

Table 2 :

 2 Properties of the roads in Figure10road length[km] ρ max [ cars km ] v low [ km h ] v high [

	km h ] initial density [ cars km ]

Table 3 :

 3 Optimization results for the network in Figure10AR, ∂pi

		∂v max i	= 0 AR, ∂pi ∂v max i	= 0 LWR
	no control	1871.7	1871.7	834.9
	ramp metering only	1325.3	1325.3	834.9
	speed control only	1122.8	872.6	834.9
	both control types	814.5	818.4	834.9

Table 4 :

 4 Total travel time for different choices of the relaxation parameter δ δ no control opt. control, ∂pi

	∂v max i	= 0 opt. control, ∂pi ∂v max i	= 0


Maria Laura

Delle Monache

Paola Goatin

A numerical scheme for moving bottlenecks in traffic flow

Keywords: Scalar conservation laws with local moving constraints, Traffic flow modeling, PDE-ODE coupling, Conservative finite volume schemes

We consider a strongly coupled PDE-ODE system that describes the moving bottlenecks created by several buses on a road. The model consists of a scalar conservation law modeling the main traffic evolution and a series of ODEs accounting for the trajectories of the buses, which depend on the surrounding traffic density. The moving bottlenecks are expressed by inequality constraints on the flux. We generalize a conservative finite volume scheme for the constrained hyperbolic PDE using a reconstruction technique and a tracking algorithm for the ODEs. Numerical simulations illustrate the impact of the buses on traffic flow. This is joint work with C. Chalons

Introduction

We consider a circular bus route of unit length, where we model traffic flow in the presence of N buses, which are responsible for moving bottlenecks on the road. This situation can be described by a coupled macroscopic-microscopic model as introduced in [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF]. The dynamics of the cars density and the slower vehicles reads

∂ t ρ + ∂ x f (ρ) = 0, (t, x) ∈ R + × [0, 1], (1a)
ρ(0, x) = ρ 0 (x), x ∈ [0, 1], (1b)
f (ρ(t, y i (t))) -ẏi (t)ρ(t,

y i (t)) ≤ F α,i t ∈ R + , i = 1, ..., N, (1c)
ẏi (t) = ω(ρ(t, y i (t)+)),

t ∈ R + , i = 1, ..., N, (1d)
y i (0) = y 0,i , i = 1, ..., N, (1e)
ρ(t, 0+) = ρ(t, 1-) t ∈ R + , (1f)
where ω(ρ(t, y i (t)+)) denotes the right-hand limit of ω(ρ(t, x)) at x = y i (t) and

F α,i . = αR 4V (V -ẏi (t)) 2 , i = 1, . . . , N. 1
Above, ρ = ρ(t, x) ∈ [0, R] is the scalar conserved quantity and represents the traffic density, whose maximum attainable value is R. The flux function f : [0, R] → R + is assumed to be strictly concave and such that f (0) = f (R) = 0.

In this paper, we will take f (ρ) = ρv(ρ), where v(ρ) = V (1 -ρ/R) is the mean traffic speed, V being the maximal velocity allowed on the road. The time-dependent variable y i denotes the position of the i-th bus. If the traffic is not too congested, each bus moves at its maximal speed V b < V . When the surrounding traffic density is higher, the bus adapts its velocity accordingly. In particular, it is not possible for a bus to overtake cars. Besides, we remark that the case of buses with different maximal speeds requires a specific treatment and is not covered in this paper. In particular, we do not consider the possibility for a bus to reach and bypass another bus running ahead. From a mathematical point of view, the velocity of the bus can be described by the following traffic density dependent function (see Figure 1) The bus acts as a moving obstacle, thus hindering the traffic flow as expressed by the flux constraint (1c). There, α ∈ ]0, 1[ is the reduction rate of the road capacity due to the presence of the bus, see [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled pde-ode problem[END_REF][START_REF] Delle Monache | A front tracking method for a strongly coupled PDE-ODE system with moving density constraints in traffic flow[END_REF][START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF] for details on the derivation of the constraint (1c). Note that the inequality is always satisfied if ẏi (t) = v(ρ), since the left-hand side is 0. Moreover, it is well defined even if ρ has a jump at y i (t) because of the Rankine-Hugoniot conditions. Denoting by S 1 = R/Z the 1-dimensional torus, the notion of solution is given as follows.

ω(ρ) = V b if ρ ≤ ρ * . = R(1 -V b /V ), v(ρ) otherwise. ρ R ρ * ω(ρ) Vb v(ρ)

Bus speed

Cars speed

Definition 1.1 A (N+1)-tuple (ρ, y 1 , . . . , y N ) ∈ C 0 R + ; L 1 ∩ BV(S 1 ; [0, R]) × (W 1,1 (R + ; S 1 )) N is a solution to (1) if 1. ρ is a weak solution of (1a), (1b), (1f), i.e. for all ϕ ∈ C 1 c (R × S 1 ; R) R + S 1 (ρ∂ t ϕ + f (ρ)∂ x ϕ) dx dt + S 1
ρ 0 (x)ϕ(0, x) dx = 0 ; moreover, ρ satisfies Kružhkov entropy conditions [START_REF] Kružhkov | First order quasilinear equations with several independent variables[END_REF] on

R + × S 1 \ N i=1
{(t, y i (t)) : t ∈ R + }, i.e. for every k ∈ [0, R] and for all ϕ ∈ C 1 c (R × S 1 ; R + ) and ϕ(t, y i (t)) = 0, t > 0, i = 1, . . . , N ,

R + S 1 (|ρ -k|∂ t ϕ + sgn(ρ -k) (f (ρ) -f (k)) ∂ x ϕ) dx dt + S 1
|ρ 0 -k|ϕ(0, x) dx ≥ 0 ;

2. for i = 1, . . . , N , y i is a Carathéodory solution of (1d), (1e), i.e. for a.e.

t ∈ R + y i (t) = y 0,i + t 0 ω(ρ(s, y i (s)+)) ds ; 3. for i = 1, . . . , N , the constraint (1c) is satisfied, in the sense that for a.e. t ∈ R + lim x→yi(t)± (f (ρ) -ω(ρ)ρ) (t, x) ≤ F α,i .
The existence of solutions for the general Cauchy problem with BV initial data is proved in [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF]. Remark that, when the constraint (1c) is enforced, a discontinuity appears, which violates the classical Kružhkov entropy condition along the bus trajectory and will be referred to as non-classical shock. For related studies on conservation laws with fixed local constraints, the reader is referred to [START_REF] Andreianov | Finite volume schemes for locally constrained conservation laws[END_REF][START_REF] Colombo | A well posed conservation law with a variable unilateral constraint[END_REF]. Model (1) was first introduced and studied in an engineering setting [START_REF] Giorgi | Prise en compte des transports en commun de surface dans la modélisation macroscopique de l'écoulement du trafic[END_REF][START_REF] Lebacque | Introducing buses into firstorder macroscopic traffic flow models[END_REF] (see also [START_REF] Daganzo | Moving bottlenecks: A numerical method that converges in flows[END_REF][START_REF] Daganzo | On the numerical treatement of moving bottlenecks[END_REF] for a numerical treatment). A rigorous mathematical framework and an existence result were later proposed in [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF]. Other systems modeling moving bottlenecks on roads are proposed in the applied mathematics literature: in [START_REF] Gasser | Vehicular traffic flow dynamics on a bus route[END_REF][START_REF] Lattanzio | Moving bottlenecks in car traffic flow: a PDE-ODE coupled model[END_REF] the flux constraint is modeled by a mollifier, thus preserving the smoothness of the flux function; a 2 × 2 system modeling traffic flow coupled with a second order ODE is presented in [START_REF] Borsche | On the coupling of systems of hyperbolic conservation laws with ordinary differential equations[END_REF]. In these works, standard numerical methods can be applied. The non classical shocks appearing in (1), instead, require a specific numerical treatment. A finite volume algorithm based on a locally non-uniform moving mesh which tracks the bus position was proposed in [START_REF] Delle Monache | A front tracking method for a strongly coupled PDE-ODE system with moving density constraints in traffic flow[END_REF]. In the present paper we generalise to the case of N buses a conservative scheme introduced in [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled pde-ode problem[END_REF] (see also [START_REF] Aguillon | A Non-dissipative Reconstruction Scheme for the Compressible Euler Equations[END_REF][START_REF] Aguillon | Capturing nonclassical shocks in nonlinear elastodynamic with a conservative finite volume scheme[END_REF][START_REF] Aguillon | Nondiffusive conservative schemes based on approximate Riemann solvers for Lagrangian gas dynamics[END_REF][START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF] for papers using a similar technique in a different setting), which precisely captures all discontinuities. Besides, the buses trajectories are computed by a tracking algorithm which accurately solves (1d), (1e), adapted from [START_REF] Bretti | A tracking algorithm for car paths on road networks[END_REF].

The Riemann Problem

We recall here the definition of the Riemann Solver for problem [START_REF] Aguillon | A Non-dissipative Reconstruction Scheme for the Compressible Euler Equations[END_REF] given in [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF] for a single bus (the index i will be dropped in this section). We denote by RS the standard (i.e., without the constraint (1c)) Riemann solver for

∂ t ρ + ∂ x f (ρ) = 0, (t, x) ∈ R + × R, ρ(0, x) = ρ 0 (x), x ∈ R,
with initial datum

ρ 0 (x) = ρ L if x < 0, ρ R if x > 0,
i.e., the (right continuous) map (t, x) → RS(ρ L , ρ R )(x/t) given by the standard weak entropy solution. Moreover, let y 0 = 0 be the initial position of the bus and ρα and ρα , with ρα ≤ ρα , be the intersection points of the flux function f (ρ) with the line f α (ρ α ) + V b (ρ -ρ α ) (see Figure 2(a)).

Definition 1.2

The constrained Riemann solver RS α is the map

RS α : [0, R] 2 → L 1 loc (R; [0, R]) defined as follows. 1. If f (RS(ρ L , ρ R )(V b )) > F α + V b RS(ρ L , ρ R )(V b ), then RS α (ρ L , ρ R )(x/t) = RS(ρ L , ρα )(x/t) if x < V b t, RS(ρ α , ρ R )(x/t) if x ≥ V b t, and
y(t) = V b t. 2. If V b RS(ρ L , ρ R )(V b ) ≤ f (RS(ρ L , ρ R )(V b )) ≤ F α +V b RS(ρ L , ρ R )(V b ), then RS α (ρ L , ρ R ) = RS(ρ L , ρ R ) and y(t) = V b t. 3. If f (RS(ρ L , ρ R )(V b )) < V b RS(ρ L , ρ R )(V b ), i.e., v(RS(ρ L , ρ R )(V b )) < V b then RS α (ρ L , ρ R ) = RS(ρ L , ρ R ) and y(t) = v(ρ R )t.
Note that, when the constraint is enforced (point 1. in the above definition), a non-classical shock arises between ρα and ρα , which satisfies the Rankine-Hugoniot condition but violates the Lax entropy condition.

2 Numerical scheme 2.1 Numerical treatment for the conservation law Approximate solutions for (1) can be computed extending the algorithm proposed in [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled pde-ode problem[END_REF] for a single bus. This scheme allows to accurately capture the non-classical shock on a fixed mesh and at the same time to track the buses positions. Let ∆x and ∆t be the fixed space and time discretizations such that we can set x j+ 1 2 = j∆x for j ∈ Z and t n = n∆t for n ∈ N. Moreover, we denote by y n i the i-th bus position at time t n , so that

y n i ∈ C mi = [x mi-1 2 , x mi+ 1 2
) for some m i ∈ Z, i = 1, . . . , N . The idea is that, whenever a non-classical shock appears, the information held by the initial data are not sufficient to correctly generate the exact solution. Therefore, it is necessary to input the left and right traces of the non-classical discontinuity in the algorithm. Whenever condition

f (ρ n mi ) > F α + V b ρ n mi , (3)
is satisfied, a non-classical shock is expected to appear at time t n near the i-th bus position. If at time t n a non-classical shock is created at y n i ∈ C mi , for some i = 1, . . . , N , then the sole information ρ n m,i , which is present in the cell C mi using a classical Godunov approach, is not enough to capture it. Hence, in the spirit of [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF], we propose to introduce in the cells C mi the left and right states ρ n mi,l = ρα and ρ n mi,r = ρα of the non-classical discontinuity which is expected to be present in the constrained Riemann solution associated with ρ L = ρ n m-1,i and ρ R = ρ n mi+1 in case inequality (3) is satisfied. Since the presence of the non-classical shock is due to the presence of the bus, we require that this reconstructed discontinuity is located inside the cell C mi at the position xmi = x mi-1 2 + d n mi ∆x, where d n mi ∈ [0, 1] is given by

d n mi = ρα -ρ n mi ρα -ρα , (4)
which implies that the proposed in-cell reconstruction is fully conservative, since

d n mi ρ n mi,l + (1 -d n mi )ρ n mi,r = ρ n mi .
Then, the numerical flux

F n mi+ 1 2 at x mi+ 1 2 is defined by ∆tF n mi+ 1 2 = min (∆t mi+ 1 2 , ∆t)f (ρ n mi,r ) + max (∆t -∆t mi+ 1 2 , 0)f (ρ n mi,l ) (5)
where ∆t

mi+ 1 2 = 1 -d n mi V b
∆x represents the time needed by the reconstructed discontinuity to reach the interface x = x mi+ 1 2 . We refer the reader to [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF][START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled pde-ode problem[END_REF] for more details.

This method allows to precisely capture the non-classical shocks introducing only one point of numerical diffusion, while guaranteeing mass conservation. To avoid instabilities due to numerical diffusion near classical shocks, which can make the scheme fail to capture non-classical solutions, we propose to adopt this strategy to reconstruct also classical shocks. We proceed as follows: given the left and right traces ρ n L and ρ n R of the shock, its position xj = x j-1 2 +d n j ∆x ∈ C j is computed from

d n j = ρ n R -ρ n j ρ n R -ρ n L ∈ [0, 1], (6)
x n mi-1

x n mi+1

x n mi and the numerical flux becomes

y n i ρ n mi-1 ρ n mi+1 ρ n mi,r ρ n mi,l
• if λ(ρ n L , ρ n R ) > 0, ∆tF n j+ 1 2 = min (∆t j+ 1 2 , ∆t)f (ρ n R ) + max (∆t -∆t j+ 1 2 , 0)f (ρ n L ), (7)
with ∆t

j+ 1 2 = 1 -d n j λ(ρ n L , ρ n R ) ∆x; • if λ(ρ n L , ρ n R ) < 0, ∆tF n j-1 2 = min (∆t j-1 2 , ∆t)f (ρ n L ) + max (∆t -∆t j-1 2 , 0)f (ρ n R ), (8)
with ∆t j-

1 2 = d n j -λ(ρ n L , ρ n R ) ∆x; • if λ(ρ n L , ρ n R ) = 0, F n j-1 2 = f (ρ n L ) and F n j+ 1 2 = f (ρ n R ).
Above, λ(ρ n L , ρ n R ) is the speed of the shock given by the Rankine-Hugoniot condition. This additional reconstruction allows to handle correctly interactions between classical shocks and non-classical waves.

Remark 1 Whenever a non-classical shock interacts with a classical one in the definition of the numerical flux, the priority is given to the reconstruction of the non-classical discontinuity.

Numerical treatment for the ODE

The buses trajectories are tracked using an algorithm based on the scheme introduced in [START_REF] Bretti | A tracking algorithm for car paths on road networks[END_REF], which accurately tracks the bus trajectory inside each computational cell. In particular, the following procedure is adopted:

• If v(ρ n mi ) ≥ V b
, then the i-th bus moves at constant speed V b . In this case, the bus trajectory is computed using the explicit Euler scheme y n+1 i = y n i + V b ∆t.

• If v(ρ n mi ) < V b , the i-th bus moves at the same speed of the surrounding traffic. In this case, we use the tracking algorithm in [START_REF] Bretti | A tracking algorithm for car paths on road networks[END_REF]. It computes the interactions between the bus trajectory and the waves produced by the Rieamnn solver at cells interfaces. A different numerical treatment is required whether the interacting wave is a rarefaction or a shock.

The complete algorithm is summarized in Algorithm 1.

Numerical results

In this section, we present two numerical tests performed with the scheme previously described to show the influence of the buses on traffic flow, and vice-versa. In both simulations, we set

V b = 0.3, α = 0.3, V = R = 1.
Test 1. We consider ρ 0 (x) = 0.4, which represent traffic in free flow equilibrium, and N = 3 buses initially located at y 0,1 = 0.2, y 0,2 = 0.4, y 0,3 = 0.6. As shown in Figure 4, the presence of buses produces permanent stop-and-go waves, thus worsening traffic conditions, with creation of queues and the passage from a regime of free flow to a regime of congested flow.

Test 2. We consider the initial datum ρ 0 (x) = 0.099 if 0 < x < 0.5, 0.99 if 0.5 < x < 1, consisting of both a free flow and a congested state. Moreover, N = 2 buses are initially located at y 0,1 = 0.45, y 0,2 = 0.5. As shown in Figure 5, the initial distance y 0,2 -y 0,1 = 0.05 between the buses gets considerably reduced due to the interaction with the shock wave. Afterwards, the buses move at the same speed, keeping their distance and creating a single moving bottleneck.

Conclusions

In this paper, we extend the numerical method introduced in [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled pde-ode problem[END_REF] to the case of several buses with the same speed law. It consists of a reconstruction based finite volume scheme that allows to capture precisely non-classical shocks for a coupled PDE-ODE problem with moving constraints. Classical and non-classical shocks are propagated exactly on a fixed mesh. We propose two numerical tests that illustrate the performances of the algorithm and are significant from the modeling point of view.

Figure 1 :

 1 Figure 1: Bus and cars speeds.

Figure 2 :

 2 Figure 2: Flux functions for ẏi = V b .

Figure 3 :

 3 Figure 3: Reconstruction of a non-classical shock.

Figure 4 :

 4 Figure 4: Test 1: Time evolution of cars density and buses positions: transition from free flow to stop-and-go waves.

Figure 5 :

 5 Figure 5: Test 2: Time evolution of cars density and buses positions: reduction of distance between two following buses.

Acknowledgments

This research was supported by the European Research Council under the European Union's Seventh Framework Program (FP/2007-2013) / ERC Grant Agreement n. 257661. This work was finalized while the second author was visiting Inria Sophia Antipolis -Méditerranée supported by Inria Associated Team "ORESTE -Optimal REroute Strategies for Traffic managEment".

Algorithm 1 Algorithm for the conservative scheme and the bus tracking Input data: Initial and boundary condition for the PDE and the ODE, m i cells indexes of the buses positions y n i , i = 1, . . . , N . Compute the densities at time t n+1 from the density values at time t n using the Godunov scheme. for j=1,. . . , J do for i=1,. . . ,N do Find j such that ρ n j < ρ n j+1 if j = m i , and

Compute the densities averages at time t n+1 using the conservative scheme

Compute the buses positions for i=1,. . . ,N do