A front tracking method for a strongly coupled PDE-ODE system with moving density constraints in traffic flow

In this paper we introduce a numerical method for tracking a bus trajectory on a road network. The mathematical model taken into consideration is a strongly coupled PDE-ODE system: the PDE is a scalar hyperbolic conservation law describing the traﬃc ﬂow while the ODE, that describes the bus trajectory, needs to be intended in a Carath´eodory sense. The moving constraint is given by an inequality on the ﬂux which accounts for the bottleneck created by the bus on the road. The ﬁnite volume algorithm uses a locally non-uniform moving mesh which tracks the bus position. Some numerical tests are shown to describe the behavior of the solution.


Introduction
The first macroscopic model for traffic flow dates back to the 1950s when Lighthill and Whitham [19] and, independently, Richards [20], proposed a fluid dynamics model to describe traffic flow on an infinite single road, using a non-linear hyperbolic partial differential equation (PDE). This model is commonly referred to as the LWR model. The Cauchy problem has then been extended to initial boundary value problem in [2]. More recently several authors proposed models that track a single vehicle moving in the vehicular traffic. In these models the single vehicle trajectory is described with an ordinary differential equation (ODE), see [4,8,11,18] and references therein. Indipendently, in the transport engeneering community, different numerical methods that approximate solutions for problems with moving bottlenecks have been developed, see [12,13]. In these works the moving constraints are replaced by a sequence of fixed ones and the discontinuity is applied at the upstream cell interface with respect to the bottleneck position. Moreover, they only deal with with triangular flux diagrams.
In this article we consider the model introduced in [15] to model the effect of urban transport systems, such as buses, in a road network. From an analytical point of view, we deal with a hyperbolic conservation law which describes the evolution of the main traffic (LWR model), an ODE which describes the bus trajectory and an inequality constraint which models the bottleneck effect created by the presence of a bus on a road. Existence of solutions to this problem for general BV data was 1 proved in [14]. This model can be seen as a generalization to moving constraints of the problem consisting in a scalar hyperbolic conservation law with a (fixed in space) flux constraint, introduced and studied in [1,9,10]. The article presents a numerical method for computing the solutions of this strongly coupled PDE-ODE systems. The results are obtained by combining a tracking algorithm in Lagrangian coordinates which uses a locally nonuniform mesh as in [21] and a tracking algorithm which reconstructs the bus position through its interaction with the density waves as in [8].
The article is organized as follows. Section 2 contains some preliminary notations and definitions. Some theoretical background and the existence theorem for Cauchy problems are introduced. In Section 3 we present a finite volume scheme with moving mesh and the tracking algorithm for the solution of the ODE. In Section 4 we present some numerical tests which show the effectiveness of our approximation.

Mathematical Model
We consider a slow moving large vehicle (e.g. a bus) that reduces the road capacity and generates a moving bottlenck for the surrounding traffic flow. This can be modeled with a PDE-ODE coupled system consisting in a scalar conservation law representing the traffic flow with a density constraint and an ODE describing the slower vehicle trajectory: x ∈ R, ρ(t, y(t)) ≤ αρ max , t ∈ R + , y(t) = ω(ρ(t, y(t)+)), t ∈ R + , y(0) = y 0 .
(1) Above ρ = ρ(t, x) ∈ [0, ρ max ] is the scalar conserved quantity denoting the mean traffic density and ρ max is the maximal density allowed on the road. The flux function f (ρ) : [0, ρ max ] → R + is a strictly concave function such that f (0) = f (ρ max ) = 0. It is given by the following flux-density relation where v is a smooth decreasing function denoting the mean traffic speed that in this article is set to be v(ρ) = V (1 − ρ ρmax ), with V the maximal velocity allowed on the road. y = y(t) represents the slower vehicle position, which moves with a speed ω(ρ(t, y(t)+)). The bus velocity w therefore depends on the downstream traffic density and is given by that it, the slower vehicle moves with a constant speed V b < V as long as it is not slowed down by the downstream traffic conditions. If this is the case, it will move at the same speed of the main traffic. The coefficient α ∈ ]0, 1[ gives the reduction rate of the road capacity due to the presence of the large vehicle. For simplicity, in the following we assume that ρ max = V = 1 so that the model becomes x ∈ R, ρ(t, y(t)) ≤ α, t ∈ R + , y(t) = ω(ρ(t, y(t)+)), t ∈ R + , y(0) = y 0 . (3) The above model was introduced in [15], in an engineering framework, to study the effect of urban transport systems in a road network. Its analytical properties were addressed in [14].

The Riemann Problem with moving density constraint
Consider (3) with the particular choice We want to define the corresponding Riemann solver with moving density constraint. To this end, we rewrite the equations in the bus reference frame (setting Figure 1, right, and we get under the constraint ρ(t, 0) ≤ α.

Remark 1 The above definition is well-posed even if the classical Riemann solution
and hence

Remark 2
The density constraint ρ(t, y(t)) ≤ α does not appear explicitly in Definition 2.1, and in the following Definition 2.2. It is handled by the corresponding condition on the flux The corresponding density on the reduced roadway at x = y(t) is found taking the solution to the equation closer to ρ y . = ρ(t, y(t))), see Figure 1.

Cauchy Problem: existence of solutions
We briefly recall the known results for this type of problem (for details see [14]). We start giving our definition of solution.
moreover, ρ satisfies the Kružhkov entropy conditions [17] on , for every k ∈ [0, 1] and for all ϕ ∈ C 1 c (R 2 ; R + ) and ϕ(t, y(t)) = 0, t > 0, Remark that the above traces exist because ρ(t, ·) ∈ BV(R) for all t ∈ R + . The complete proof can be found in [14]. It consists in constructing a sequence of approximate solutions via the wave-front tracking method, and prove its convergence and then check that the limit functions satisfy conditions (9a)-(9d) of Definition 2.2.

Numerical scheme
Our aim is to present a numerical method to compute solutions to strongly coupled constrained PDE-ODE problems. Since the solutions of the Riemann problem are known explicitely, it is natural to develop a Godunov-type method. The standard Godunov method, in principle, could be applied, however, the results produced are not correct, since it will not reproduce all the characteristics of the solutions and it fails to capture the presence of the non-classical shock. This can be overcame by applying a front tracking-capturing method which uses a Lagrangian algorithm in which the interface is tracked, such as the one in [21], together with a numerical method that tracks at each time step the slower vehicle trajectory, as in [8].
3.1 Godunov-type scheme for hyperbolic PDEs with constraint First, we briefly recall the classical Godunov scheme and then we show how it is modified to fulfill our needs. We use the following notation: x n j are the grid points at time t n with n ∈ N and j ∈ Z. A computational cell is given by [ is the cell size at time t n . The Godunov scheme [16] is based on exact solutions to Riemann problems. The main idea of this method is to approximate the initial datum by a piecewise constant function, then the corresponding Riemann problems are solved exactly and a global solution is simply obtained by piecing them together. Finally, one takes the mean on the cell and proceeds by iteration. Given U (t, x), the cell average of U in the cell j and at time t n is defined as Then the Godunov scheme consists of two main steps: 1. Solve the Riemann problem at each cell interface x n j+ 1 2 with initial data (U n j , U n j+1 ). 2. Compute the cell averages at time t n+1 in each computational cell and obtain U n+1 j .
We remark that waves in two neighbouring cells do not intersect before t n+1 = t n + ∆t n if the following CFL condition holds: where λ n j+ 1 2 is the wave speed of the Riemann problem solution at the interface x j+ 1 2 at time t n . The classical Godunov scheme can be expressed in conservative form as where F (U, V ) = f (R(U, V )(0)) is the corresponding numerical flux. Since our aim is to track the trajectory of the bus using a Lagrangian algorithm a moving mesh has to be used. In particular, we develop an algorithm which follows at each time step the bus trajectory and modifies the mesh when and the constraint is enforced. This factor needs to be considered when considering a Godunov scheme. In particular, if the inequality (13) is satisfied, then the solution of the Riemann solver is not the classical one and hence, the Godunov scheme cannot be applied as it is. We are going to shift grid points locally and, as a consequence, we will have a locally non-uniform mesh due to a cell interface moving with the bus trajectory. We will use the superscript new to indicate the quantities that are modified at time t n with the grid. Assume that at time t n , y n is the bus position and y n ∈ ]x n m− 1 2 , x n m+ 1 2 ], for some m. When (13) holds, the algorithm for the adaptive mesh reads as follows: Each time the constraint is enforced the bus position follows the non-classical shock trajectory: y n+1 = x n+1 The other cell interfaces are kept unchanged. An explicit formula for the scheme can be derived in the following way. Consider for example the finite volume cell C m−1 in Figure 4 (abcd). Integrate the conservation law over the finite volume: For semplicity, we introduceF which is given bỹ so thatF corrisponds to the F α computed in (7) when the constraint is active and to the calssical Godunov flux when the constraint is not enforced. After the mesh has been resized and adjusted as described earlier in this section, we update the cell averages for all cells with the following conservative formula:

Numerical method for the ODE
We expose here how to solve numerically the ODE, keeping track of the bus position. At each time t n we determine the position y n of the driver by studying interactions between the bus trajectory and the density waves within a fixed cell. We distinguish the two cases: • Inequality (13) is satisfied. Then the bus moves with fixed velocity V b and we update the bus position y n+1 = y n + V b ∆t n .
• Inequality (13)  [. In both cases, we check if the wave starting at the cell interface is a shock or a rarefaction and compute the time of interaction between the wave and the bus trajectory. In the case of the rarefaction the initial and final time of interaction is computed and the position of the bus is updated by solving explicitely an ordinary differential equation. According to the new position of the bus, the cell index is updated.

Numerical algorithm
In this section we describe in detail the algorithm which is composed of the following steps: Compute the densities averages at time t n+1 using formula (19). Compute the bus position:

Numerical results
For illustration, we choose a concave fundamental diagram with the following flux function: In this case the Godunov numerical flux is given by with ρ cr = 0.5 the density at which the unique maximun of the flux function is attained, f (ρ cr ) = f max . In this section we present two numerical tests performed with the scheme previously described. We deal with a road of length 1 parametrized by the interval [0, 1]. In both the simulations we fix V b = 0.3, α = 0.6.
1. Case I: We consider the following initial data The solution is given by two classical shocks separated by a non-classical discontinuity, as illustrated in Figures 5 and 6.   2. Case II: We consider the following initial data ρ l (0, x) = 0.8, ρ r (0, x) = 0.53, y 0 = 0.5.
The values of the initial conditions create a rarefaction wave followed by a non-classical and a classical shocks on the density, as illustrated in

Conclusions
This article introduces a numerical method for a strongly coupled PDE-ODE problem with moving density constraint which can model the presence of a moving bottleck on the roads. The PDE describes the evolution of the main traffic in time while the ODE describes the bus trajectory. The theoretical framework of the model was set up. The density is computed using a Godunov-type scheme with a locally nonuniform mesh. Then the position of the bus is recostructed determining the effects of the interactions with density waves as in [8]. Some numerical tests are presented to show the effectiveness of the scheme.