Numerical discretization of Hamilton-Jacobi equations on networks

We discuss a numerical discretization of Hamilton--Jacobi equations on networks. The latter arise for example as reformulation 
of the Lighthill--Whitham--Richards traffic flow model. We present coupling conditions for the Hamilton--Jacobi equations 
and derive a suitable numerical algorithm. Numerical computations of travel times in a round-about are given.


1.
Introduction. Traffic flow models, especially on road networks, have been intensively studied in the mathematical [3,6,13,16,20] as well as in the engineering community [1,7,9,10,14,18,19,27,28,30] during the last years. We are interested in first-order macroscopic models based on partial differential equations for the traffic density [24,26] with the prototype being the Lighthill-Whitham-Richards (LWR) model. When considering a traffic network the crucial point is the coupling at a traffic junction leading to coupling conditions. Well-posedness results for those conditions have been obtained for example in [15]. When considering car trajectories the Hamilton-Jacobi (HJ) reformulation of the LWR model can be used. The discussion of traffic model in HJ form has been introduced in engineering literature as for example [8,10,30]. Therein [29] the solution is also known as Moskowitz function. Recently, there has been an intense discussion on analytical properties of this equation in the context of traffic flow networks. For example, in [2] the HJ formulation has been used to deduce optimal starting times for a congested road among other results. In [22] general coupling conditions for systems of HJ equations have been studied analytically. A particular application of the well-posedness result therein is the coupling condition for traffic flow using a fixed ratio for incoming and outgoing traffic flows. For question of reconstruction of parameters of traffic models the HJ approach has been succesfully discussed and applied in [1,28].
We are interested in a numerical scheme combining the existing coupling conditions for LWR models [15,21] stated in the density and flow variables with a suitable numerical method for the HJ equation. The advantage of the numerical computation of the HJ model is to directly represent the trajectories of particular traffic 686 SIMONE GÖTTLICH, UTE ZIEGLER AND MICHAEL HERTY members and to compute the duration of specific journeys through the network. Further, we introduce a new coupling condition for merging junctions replacing the common right of way parameter [5,6] by a priority rule. For the numerical method we extend the numerical algorithm of [23]. Other numerical approaches have been discussed in [1,10,30].
2. Traffic flow network model. We give a brief review on the LWR traffic flow model [26,31] for road networks and consider the coupling conditions for several specific junctions.
Here, ρ : (x, t) → ρ(x, t) ∈ [0, ρ max ] ⊂ R + denotes the density of cars, x ∈ [0, L] ⊂ R + describes the location on the road, L is the length of the road (possibly being ∞) and t ∈ R + is time. As in [15,21] we assume that the flux function (also called fundamental diagramm) f (ρ) is concave with a unique maximum at the designated point ρ * ∈ [0, ρ max ]. This allows to define the function τ : [0, ρ max ] → [0, ρ max ] mapping the density to a distinct density value with equal flux: We denote by f −1 where v(ρ) = v max ρ max (ρ max − ρ), is the velocity which cars are assumed to have depending on the actual traffic density, v max is the maximal allowed velocity for the road and ρ max is the maximal traffic density (bumper-to-bumper density) [3]. Another broadly used function with constant velocity λ is the triangular flow function [1,10,28] (and in the context of telecommunication networks [12]): A road network is a directed graph G(E, V ), where E denotes the set of edges which represent the roads and V the set of vertices which represent the traffic intersections or junctions. The length of each road i, leading from one junction to the next, is given by L i ∈ R + . The roads are unidirectional. Lanes for different directions can be described by separate edges. Incoming roads of each junction v are denoted by δ in v and outgoing by δ out v . The density at the junction for an incoming edge i will be denoted byρ i (t) = ρ i (L i −, t) and the density for an outgoing edge j isρ j (t) = ρ j (0+, t). As coupling condition the conservation of mass is imposed: For simplicity we use the following notation for the flow at junctions:γ i := f (ρ i ) for incoming edges andγ j := f (ρ j ) for outgoing edges. At dispersing junctions (junctions with more than one outgoing road), we assume to have given (possibly time-dependent) distribution parameters 0 ≤ α ij ≤ 1 indicating the percentage of the cars coming from road i going to road j. Obviously, i∈δ out v α ij (t) = 1. The condition (4) does not guarantee a unique traffic densities at the junction. As in [3,4,15,21] the following construction is used to obtain unique traffic densitiesρ i ,ρ j at the junction given the respective traffic flowsγ i andγ j . For sake of clarity, we skip the time variable in the following notation. The densities at the junction areρ i andρ j , respectively. The admissible sets for those densities depend on the initial densities at the junction (ρ i (L i ) and ρ i (0), respectively) and are then for incoming roadŝ and for outgoing roads Givenγ i (orγ j ) the uniquely defined densitiesρ i (orρ j ) within the admissible sets areρ respectively. We refer to Figure 1 for an example.
and γ max In the following paragraphs, we explicitly state coupling conditions in order to obtain unique traffic flowsγ i andγ j . Those, by the above discussion, then lead to unique traffic densities which in turn can be used as boundary conditions for (1).

SIMONE GÖTTLICH, UTE ZIEGLER AND MICHAEL HERTY
Coupling conditions for traffic fluxes have been proposed by [3,15]. Here, we focus on priority rules and round-abouts.
Two connected roads. In a bottleneck situation, the capacity of the traffic load may decrease at a certain point. This can be viewed as a network model with two connected roads as depicted in Figure 2 such that max ρ f j (ρ) < max ρ f i (ρ) : where γ max 1 and γ max 2 are given by (9) and (10), respectively. The unique solution isγ Dispersing junction. We assume the distribution rate α ij for i = 1, j ∈ {2, 3} at the junction to be previously known. As in [9] we assume, that all the flow at the junction is restricted as soon as one of the outgoing roads is not able to absorb the designated flow. This corresponds to a first-in-first-out rule of cars and is a realistic assumption, since a car waiting at the junction blockes all the traffic behind until it continues. A mathematical model respecting this property stated in terms of the flowsγ i andγ j at the vertex is: This linear programming problem is solved explicitly: where γ max 1 is given by (9) and γ max 2 and γ max 3 are given by (10). The resulting boundary densities at the junction are again given by (7) and (8).  Merging junction. In contrast to the previous discussion problem (13) may have multiple solutions in the case of a merging junction depicted in Figure 4(a). Therefore, additional conditions have to be imposed as constraints to (13). In [3,5,9,20] propose a right of way parameter q ∈]0, 1[ prescribing the proportion of flow coming from road 1 and 2. Here, we formulate a priority rule, where the traffic of the main road always is prioritized over the traffic of a side road. As soon as road 3 reached a state of dense traffic, cars from road 1 are preferred. As a mathematical formulation we propose to replace (13) by (15) where the prioritization of the flow coming from road 1 is obtained by using a weighting parameter ω > 1 in the objective function: Lemma 2.1. There exists a unique solution of (15) given bȳ where γ max 1 and γ max 2 are given by equation (9) and γ max 3 is given by equation (10).
The proof of Lemma 2.1 can be found in Appendix A. Round-about. We consider a round-about depicted in Figure 5(a). For simplicity we assume it is composed of four junctions with two incoming and two outgoing roads each. Additionally, the central ring of the round-about has priority over the connecting roads, e.g., road 1 is prioritized over road 2 in Figure 5(c). Furthermore, we assume road 2, no car is going to road 3, but all go to road 4. Hence, the distribution parameters for the dispersing junction are α 2,3 = 0 and α 2,4 = 1.
Finally, the traffic distribution from road 1 to road 3 and 4 is also prescribed by α 1,3 and α 1,4 , respectively. Therefore, the coupling condition for one part of the round-about is  Lemma 2.2. Let α 2,3 = 0 and α 2,4 = 1. If α 1,3 α 1,4 = 0, then, problem (17) has a unique solution. the solution is given bŷ The proof of Lemma 2.2 can be found in Appendix B. Note, that in [3,4,6] the considered distribution parameters α have to be strictly larger than zero and strictly smaller than 1. The proof of Lemma 2.2 especially considers the case, where α 2,3 = 0 and α 2,4 = 1.
3. Numerical scheme for the Hamilton-Jacobi reformulation. As in [8,10,28,29], the traffic network model can be reformulated as HJ equation. In this section, we briefly recall the relation between LWR and HJ, and derive a numerical scheme [23]. A HJ equation with Hamiltonian f is given by If we consider roads on which vehicles cannot overtake, it is possible to number them according to the order, they pass a certain point of the road. In [4,27,29,30] a continuous function is considered, where the space-time trajectory of each car is given by its the curves of cumulative counts. In detail, if we start counting with the foremost car at time t = t 0 we get and for a general point in time t, the car number at (x, t) is given by where the value of the left boundary is given by Consequently, the curve given by (19) the continuity equation used in the LWR-model (1) holds. Differentiation of (19) with respect to x yields: Consequently, if we find an M that satisfies (19), ρ := M x also satisfies (1). On traffic problems we have ρ ≥ 0, hence M is monotonically increasing in x.
Extension to the network case. For the network model of HJ equation we provide the additional index e indicating the current road. We have e ∈ {1, ..., |E|}. Then, the complete model reads right boundary condition . (21) The values ofρ e andρ e are obtained by the coupling conditions. They therefore depend on the adjacent edges of the node. Within the coupling conditions the car density ρ e (x, t) for x = 0 and x = L has to be evaluated. Therefore, in the numerical algorithm we first reconstruct at every time step t n the density ρ(·, t n ). Using the coupling conditions we then compute theρ e (t n+1 ) andρ(t n+1 ) to finally obtain the boundary conditions for M e .
3.1. Numerical scheme for a single road. Consider a single road first. We introduce a spatial grid i ∈ {0, . . . , n x }, x i = ∆x · i, n x = L ∆x , and a temporal grid t ∈ {0, . . . , n t }, where the gridsize ∆t is set according to the CFL-condition for piecewise differentiable functions: The superindex denotes the time step and the subindex the spatial point. We set ρ t i := ρ(∆x · i, t · ∆t) and M t j = M ((j − 1 2 )∆x, t · ∆t) where j = {0, n x + 1}. Note that the grid points for the discretization of M (x, t) are shifted by ∆x 2 c compared to the grid for ρ(x, t). For the discretization of HJ we use the central first order scheme [23]. The time evolution of M at the inner grid points 0 < j < n x + 1 is computed as follows: and due to the CFL condition a n j ≥ max The coupling is done in terms of densities (at least at the boundary of the domain). Therefore we need to reconstruct the derivative of M at the junction. This is done by finite differences and we define the density ρ t i by where (j − 1 2 )∆x = i∆x. We have to add an additional index e for all e ∈ E to M and ρ, respectively, when considering the network formulation. However, for the sake of readability, we will drop this index whenever the context is clear.

3.2.
Discretization of the boundary condition. Due to dispersion effects of the discretization scheme [25], the densities at the junction might not always be captured correctly. Therefore we need to introduce suitable ghost-cells added on both ends of each road. The wave fronts travel along the roads until they reach the next junction, providing the coupling condition with information about the actual density on the road. We let the waves run through these artificial cells, but take the value at the road boundary to compute the coupling condition ( Figure 6). This finally leads to the correct density information at the junction, see below for more details.
This method only works, when the number of ghost-cells is large enough to absorb the dispersion of the wave front. In the sequel we will show, that two ghost-cells on each side of the road are sufficient. A correlation between spatial and temporal grid and a particular choice of parameter a n j is required to ensure this result. Lemma 3.1. If the parameter a n j of the HJ-Scheme (23) is a n j := max and the grid size ∆t is set to the maximal possible value satisfying the CFL-condition (22), then the HJ scheme (23) is equivalent to the Lax-Friedrich scheme: Proof. Scheme (23) and equation (24) allow for the following calculation: which is precisely the Lax-Friedrich scheme (26).

Lemma 3.2.
Assume that f is given by equation (3) and assume the density at time t is of Heaviside type. i.e., ∃x ∈ [0, L] such that ρ(x, t) = l, x ≤x r, x >x.
Then, for the Lax-Friedrich scheme and the dispersion over time of the wave front will not exceed two grid cells.

HAMILTON-JACOBI ON NETWORKS 693
Proof. The Lax-Friedrich scheme yields The spatial grid is given such that the discontinuity of the initial condition is located between grid point i and gridpoint i + 1. Hence, the density values at the discontinuity at timestep t are given by: , . . . , r).
The scheme preserves the density values inside the constant regions, because (28) for an arbitrary space grid point j. Hence, it is sufficient to consider the density evolution next to the discontinuity. For this purpose we distinguish several cases: Applying (28) to ρ t , we get Hence, we get a sharp forward traveling front without any dispersion.
Applying (28) to ρ t , we get Hence, ρt +1 again fulfills the assumptions imposed to ρt, with the shape shifted by one space step either to the left or to the right. Therefore, by applying the Lax-Friedrich scheme iteratively over time, the dispersion will never become greater than two cells.
Hence, the resulting wave front is moving backwards carrying along two middle density values ρ * .

3.3.
Details of the numerical algorithm in the network case. The complete numerical scheme for solving HJ equations on road networks is described in Algorithm 1. Some steps are particularly illustrated in Figure 6. A crucial point in the computation are the coupling conditions as depicted in Figure 6(c). As denoted in line 16 of Algorithm 1, equations (9) to (18) are used. The detailed procedure is the following: Consider a junction v with at most two incoming roads (∈ δ in v ) and at most two outgoing roads (∈ δ out v ). The leftmost gridpoints of the incoming roads and the rightmost gridpoints of the outgoing roads in terms of the density ρ have been computed for timestep t + 1, see Figure 6(b). Hence, the values for ρ t+1 e,nx ∀e ∈ δ in v and ρ t+1 e,0 ∀e ∈ δ out v are given corresponding to ρ e (L) and ρ e (0) in the continuous case. Now, we use equations (9) and (10) to obtain the maximal possible flow γ max e for all roads e connected to the junction. Depending on the junction type we compute the coupling flowsγ e ∀e ∈ δ in v and γ e ∀e ∈ δ out v using equations (12), (14), (16) or (18), respectively. The density boundary valuesρ e ∀e ∈ δ in v andρ e ∀e ∈ δ out v are uniquely given by (7) and (8). An illustration of this procedure is given in Figure 7. Some further remarks on the algorithm 1 in order: line 2. Note that the Godunov scheme [17] does not need any ghostcells to compute the coupling condition. The presented scheme introduces numerical diffusion such that the ghost cells need to be sufficiently large. Its size has been discussed in the previous Lemma. The length of the ghost cells is chosen equal to the size of the interior cells. In order to have the same speed of propagation those cells do not enter the computation of the length of the road.

4.
Computational results for round-abouts. We apply the merging and dispersing junction model to a small network consisting of eight roads, see Figure 8. The network describes a small traffic round-about examined in [3]. We use the same setting as in [3], where the flow is given by f (ρ) = ρ(1 − ρ) and initial as well In [3] this test case is compared for different right-of-way parameters q ∈]0, 1[, determining the proportion of cars coming from each road at merging junctions. The priority rule used in this paper corresponds to the case q = 0. Figure 8 shows the traffic density exemplarily for four roads at four different points in time. Since the boundary condition is constant, the density evolution reaches an equilibrium and does not change for t > 5. The traffic at the inner circle has priority, therefore, the round-about does not get blocked. This is qualitatively the same behaviour as in [3], when a small parameter q is used. Since our model uses strict priorities, the equilibrium state is reached faster compared with [3]. We have used the proposed method for simulations and reconstructed the density values as in (24). The thin lines show the result obtained by using the Godunov scheme, which we use as a benchmark. For road 4 at time t = 2 it can clearly be seen, that by Godunov a sharper resolution of the shock wave is obtained. However, the actual density levels are equivalent for both schemes.
Next, we consider a larger round-about composed of four junctions with two incoming and two outgoing roads as in Figure 9(a).
According to the enumeration, roads 1 to 4 are leading towards the inner circle which is composed of roads 5 to 8. Roads 9-12 point out of the round-about. Usually, drivers cannot drive as fast inside the inner circle as on the other roads. For this reason, we describe these roads by different triangular fundamental diagrams. As before, the traffic density belongs to the interval between 0 (no traffic) and 1 (maximal dense traffic). Since we assume that the usual speed of the cars is faster at the outer roads than inside the circle, the corresponding flow function has a steeper slope outside the inner circle. We prescribe the left boundary data for the incoming roads 1 to 4. We assume, that road 1 and 3 are slightly more busy than roads 2 and 4. For simplicity we use the same boundary data for each road pair. Figure 10 gives a detailed overview of the boundary data at an average working day from 5am to 1pm. This is a test setting attempting to tackle the qualitative traffic behaviour taking the morning rush hour into account. Figure 11 shows the traffic density along the inner circle for exemplary points in time. Since the traffic at the inner roads always has the priority at junctions and outgoing roads are not blocked in our setting, no jams appear inside the roundabout. However, we observe, that at peak time the traffic density all along the inner circle is at value ρ * = 0.5, which means that the traffic moves to maximal possible flow. Also at peak time, traffic jams occur at roads leading to the inner circle. Particularly from 7am to shortly after 11am, the traffic entering the roundabout is quite dense. However, since the incoming traffic reduces drastically around 11am (see boundary condition depicted at Figure 10(a)) the jam is resolved again a while after the incoming traffic reduces.
It is easy to derive the trajectories of cars from the HJ, since we only have to compute the contour lines of function M e (x, t). In Figure 13, the trajectories of 3 cars moving along the roads 1-5-6-11 are depicted. When entering the traffic network before 6:59am the cars move freely and leave the system already about 1 minute. In contrast to that, another driver, who enters the system only 4 minutes later, already encounters dense traffic on the road and needs more than 4 minutes to move to the end of road 11. Figure 14 shows the duration of the route 1-5-6-11 depending on the starting time of the journey. While it takes only 1 minute to traverse the route during light traffic times, cars need up to 4.7 minutes between  7 and 9 am. Hence, it takes more than 4 times longer to traverse the given route during the rush-hour.     Figure 14. Travel time (in minutes) for the route depicted in Figure 13, depending on the starting time of the journey.
. The feasible region of each case is depicted in Figure 15 and the optimal solution is indicated by the black dot. This yields directly (16a) -(16c). Hence, (15) is uniquely solvable.