An easy-to-use algorithm for simulating traffic flow on networks: numerical experiments

In this paper we propose a Godunov-based discretization of a hyperbolic system of conservation laws with discontinuous flux, modeling vehicular flow on a network. Each equation describes the density evolution of vehicles having a common path along the network. We show that the algorithm selects automatically an admissible solution at junctions, hence ad hoc external procedures (e.g., maximization of the flux via a linear programming method) usually employed in classical approaches are no needed. Since users have not to deal explicitly with vehicle dynamics at junction, the numerical code can be implemented in minutes. We perform a detailed numerical comparison with a Godunov-based scheme coming from the classical theory of traffic flow on networks which maximizes the flux at junctions.

1. Introduction. Starting from the introduction of the LWR model [23,25], a huge literature about macroscopic fluid-dynamic models for traffic flow was developed. More recently, models, theory and numerical approximations for traffic flow on networks became a hot topic [8,11,14,20]. The interest in forecasting traffic flow on large networks became even stronger in the very last years, due to the increasing number of GPS devices (smartphones, satellite navigators, black boxes) which provide real-time traffic data. Private companies like GOOGLE, WAZE MOBILE, NOKIA [16], INRIX, OCTOTELEMATICS [9], ZEROPIU, YANDEX, started collecting data and, in some cases, broadcasting traffic forecast.
The model. In this paper we study from the numerical point of view a firstorder version of the model proposed in [19,Sect.4]. Let us consider a network, i.e. a directed graph with N R arcs (roads) and N J nodes (junctions). Vehicles moving on the network are divided on the basis of their desired path. Let us assume that the number of possible paths on the graph is N P and denote those paths by P 1 , . . . , P p , . . . , P NP . We stress that paths can share some arcs of the networks. A point x (p) of the network is characterized by both the path p it belongs to and the distance x from the origin of that path. We denote by µ p (x (q) , t) the density of the vehicles following the p-th path at point x (q) at time t > 0, and we assume that µ p (x (q) , t) ∈ [0, ρ max ] for some maximal density ρ max . Note that we have, by definition, µ p (x (q) , t) = 0 if x (q) / ∈ P p . We also define i.e. ω(x (p) , t) is the sum of all densities living at x (p) at time t. The function ω is discontinuous and takes into account the topology of the network. Note that, for any point x (q) , the densities µ p (x (q) , t), p = 1, . . . , N P , are admissible if their sum ω(x (q) , t) ≤ ρ max . Let us denote by v(ω) the velocity of vehicles (given as a function of the density) and by f (ω) = ωv(ω) the flux of vehicles. The LWR-based model is constituted by the following system of N P conservation laws with space-dependent and discontinuous flux ∂ ∂t µ p (x (p) , t) + ∂ ∂x (p) µ p (x (p) , t) v ω(x (p) , t) = 0, x (p) ∈ P p , t > 0, (2) or, equivalently, for p = 1, . . . , N P . If ω = 0 we have, a fortiori, µ p = 0, then it is convenient to set µ p ω = 0 in (3) to avoid singularities. In the following we assume that the flux f ∈ C 0 ([0, ρ max ]) ∩ C 1 ((0, ρ max )) and Equations of the system (2) (or (3)) are coupled by means of the velocity v, which depends on the total density ω. On the other hand, not all the equations of the system are coupled with each other because paths do not have necessarily arcs in common.
It is plain that the number of equations in the system grows rapidly when the number of nodes of the graph increases, making unfeasible the computation of a numerical solution. To keep the computational load within reasonable limits, we also propose another version of the model which splits the vehicles on the basis on their path only at junctions. In this way, along the arcs a single equation for the total density ω is considered. The price to pay is that the global behavior of drivers is lost.
Relevant literature. The multi-path model described above differs from the socalled multi-population or multi-class models, see, e.g., [2,27]. In those cases, the models consist of one equation for a single road (extension to networks is also possible) with different velocity functions v i , one for each class of vehicles. Typically, the populations have different maximal velocities, in order to take into account different types of vehicles or drivers' behaviors.
An apparently similar model is the one presented in [13] (see also [21]). In that case, vehicles are divided in different populations on the basis on their sourcedestination pair. Given the total density ω of all vehicles, the density µ p of the p-th population is given by where π p (x, t, O, D) specifies the percentage of the total density that started from source O, it is moving towards the final destination D, and it is at x at time t. Moreover, O(p), D(p) are the origin and the destination associated to the p-th path, respectively. Then, a standard PDE for ω is coupled with a system of PDEs for the coefficients π p 's.
Several papers investigate from the theoretical point of view the (systems of) scalar conservation laws. The interested reader can find in the book [22] an introduction to the field, in [1,6] some references for the case of discontinuous flux, and in the book [3] the analysis of the systems of conservation laws. Systems of scalar conservation laws with discontinuous flux are instead less studied. An attempt related to traffic flow can be found in [24], where a model very similar to the one considered here is investigated. From the numerical point of view, a good basic reference is again the book [22]. We also point out the paper [18], where a numerical method for (systems of) scalar conservation laws with discontinuous flux is proposed, and the paper [26], where the convergence of a Godunov-based scheme for scalar conservation laws with discontinuous flux is investigated.
Goal. The goal of this paper is proposing an algorithm to solve numerically the system (3). The numerical method is based on the Godunov scheme. The algorithm does not require use of ad hoc procedures (e.g., linear programming) to solve the dynamics at junctions, as it is done in classical approaches [4,11,14,18]. This leads to a very simple algorithm which can be implemented in minutes. We perform a detailed comparison with the algorithm described in [4,14] based on the LWR model and the Godunov scheme. We show that the two algorithms do not coincide in general since the proposed one does not always maximize the flux at junction. Nevertheless, the solution is admissible in the framework of the classical theory and in some cases it is expected to provide a more reasonable behavior of the car flow. A theoretical investigation of the algorithm proposed here can be found in [5].
Paper organization. The paper is organized as follows: In Section 2 we recall some basic facts of the classical theory of traffic flow on networks and Godunov discretization, following [4,8,14]. In Section 3 we present the algorithm for small and large networks. In Sections 4, 5 and 6 we report the numerical results and comparisons for the junctions with 2 incoming roads and 1 outgoing road, 1 incoming road and 2 outgoing roads, and 2 incoming roads and 2 outgoing roads, respectively. In Section 7 we report the result of an experiment performed on a real network in Rome, Italy. Finally, in Section 8 we sketch some conclusions.
2. Classical theory. In this section we recall some basic facts of the classical theory of traffic flow on networks, including a numerical approximation based on the Godunov scheme. Our main references are [4,8,14].
2.1. Basic definitions, assumptions, and results. Let I = [a, b] be a generic arc of the graph, i.e. a road. At any time t, the evolution of the vehicle density on the network is computed by a two-step procedure: First, a classical conservation law is solved at any internal point of the arcs; Second, the densities at endpoints a, b which correspond to a junction are computed. The latter step has not in general a unique admissible solution, so that additional constraints must be added. Beside the conservation of vehicles at junctions, we assume here that drivers behave in order to maximize the flux at junctions and that incoming roads are regulated by priorities (right of way). The second step is performed by a linear programming method.
On each arc, the density ρ(x, t) of all vehicles (no distinction among vehicles is considered here) is given by the entropic solution of We consider now a generic junction J, with I 1 , . . . , I n incoming roads and I n+1 , . . ., I n+m outgoing roads. We assume that the choice of the outgoing road is prescribed by a matrix of preferences, where 0 ≤ α j,i ≤ 1 for every i ∈ {1, . . . , n} and j ∈ {n + 1, . . . , n + m}, and n+m j=n+1 α j,i = 1, i = 1, . . . , n. The i-th column of A describes how the traffic from I i distributes in percentage to the outgoing roads.
A basic requirement for a vector (ρ 1 , . . . , ρ n+m ) to be an admissible solution to the problem at J is that which translates the fact that the vehicles are conserved at junction. Note that (6) can be seen as a generalization of the Rankine-Hugoniot condition at junctions. We define and The quantities γ i max (ρ(b i , t)) and γ j max (ρ(a j , t)) represent, respectively, the maximum incoming and the maximum outgoing flux that can be obtained on each road. Then we define The sets Ω i , Ω j contain all the possible fluxes for the solution at junctions and then the set Ω contains all the possible admissible fluxes at the end of the incoming roads, taking into account the matrix of preferences A. Since we want to maximize the flux at junction, we find the solution(s) of the maximization problem with linear constraints Note that the problem (12) has not in general a unique solution. For example, if n = 2 and m = 1 (two incoming and one outgoing roads) we have A = (1 1), and the constraint reads as γ 1 + γ 2 ≤ γ 3 max . If γ 1 max + γ 2 max > γ 3 max , we have infinite solutions. To fix this, the additional constraint can be introduced, where q = (q 1 , . . . , q n ) are the so-called priorities (right of way) coefficients such that q i ≥ 0 ∀i and i q i = 1. Equation (13) translates the fact that some incoming roads have priority with respect to the other roads in assigning their flux to the outgoing roads. The constraint (13) guarantees uniqueness of the solution of problem (12) (unless a further projection onto Ω is needed, see [7,14] for details). Finally, define (γ * 1 , . . . , γ * n ) as the unique solution of problem (12) (with the additional constraint (13) if necessary), and, consequently, 2.2. Numerical approximation by the Godunov scheme. Let us describe briefly the Godunov scheme for solving a conservation law of the form Let us denote by (x k , t n ) := (k∆x, n∆t), k ∈ Z, n ∈ N, the generic grid node and by ρ n k the density ρ at (x k , t n ). Once the initial conditionρ has been projected on the grid, in the internal nodes of the interval [a, b] the density at time t n+1 is given by where the numerical flux g is defined as At the boundary nodes we proceed as follows: • If the node is not linked to any junction, then we assign the desired boundary condition (Dirichlet or Neumann). • If the node is a right endpoint and corresponds to the i-th incoming road of a junction, we use the maximal flux (14) and set • If the node is a left endpoint and corresponds to the j-th outgoing road of a junction, we use the maximal flux (15) and set 3. Numerical approximation of the multi-path model on small and large networks. In this section we present a numerical approximation for the system (1),(3) based on the Godunov scheme (17). Then, we discuss how the multi-path model can be modified to deal easily with large networks.
3.1. Numerical approximation. Let us denote by µ n,p k (q) the approximate density is the k-th node along the path P q . Then, analogously to (1), we define From now on, to avoid cumbersome notations, we write k p instead of k (p) . Computing properly ω n k p at every node is the sole part of the algorithm which requires some effort. Indeed, the rest of the algorithm is constituted by the imposition of the boundary conditions at the beginning and the end of the N P paths and by the computation of the discrete solution at any internal nodes k p by means of the following conservative scheme: for n > 0 and p = 1, . . . , N P . Note the intrinsic asymmetry of (22). The coefficients in front of the fluxes involve only the nodes k p and k p − 1, and not k p + 1.
We stress again that no special management of the junctions is needed. The simplicity of the scheme is the main strength of this approach, making it possible to simulate traffic flow on networks in minutes.

3.2.
A local model for large networks. If the network under consideration is small, the number of possible paths is limited. Then, the number of equations in the system (1),(3) fits a manageable size. This is true even if we deal with large networks where many paths are impracticable or negligible. Conversely, if the network is large and has a large number of paths, the computation of ω n k p becomes a hard task. In this case, we adopt a hybrid point of view, creating an algorithm which merges the features of the multi-path algorithm described above and the classical algorithm described in Section 2. In the internal nodes of the arcs a single equation for the total density ρ (or ω) is solved. Then, just before each junction, vehicles are split on the basis of their direction by means of the matrix of preferences A, see (5), and the scheme (22) is applied. Just after the junction, densities associated to the same arc are summed again. Considering that typical real junctions have at most three incoming and three outgoing roads, we have at most nine different paths at each junction.
As already noted in [10,13], the original (global) multi-path approach of Section 3.1 differs from the hybrid (local) method, because splitting the density just before the junctions causes the lost of the global behavior of the drivers. To better understand the difference, let us consider the network in Fig. 1. Let us consider only two paths, If, at some time t, vehicles along [a 2 , b 2 ] are stopped by, say, an accident, the two algorithms will result in two different outcomes in [a 5 , b 5 ]. The global algorithm will empty the arc [a 5 , b 5 ], while the local algorithm will fill [a 5 , b 5 ] with a percentage of PSfrag replacements We also note that the classical algorithm described in Section 2 is local, while the one presented in [13] is global.

4.
Two incoming roads and one outgoing road. In this section we consider a network with three arcs and one junction, with two incoming roads and one outgoing road. On this network two paths P 1 and P 2 are defined, see Fig. 2. We denote PSfrag replacements J -1 Figure 2. A network with 3 arcs and 1 junction. Path P 1 (left) and path P 2 (right).
by J the node just after the junction, see Fig. 2. Note that the nodes before the junction, namely J − 1, J − 2, etc., can refer to one path or the other one, depending on the context. We have Let us first prove that the density ω n k p is admissible for any n and k p , p = 1, 2, i.e. it does not exceed ρ max . This is a crucial property which the scheme must satisfy. It is plain that the risk of exceeding the maximal density is only at node J. In the other nodes the result comes from the usual properties of the Godunov scheme.
Theorem 4.1. Let the initial densities around the junction be admissible, namely If the following CFL-like condition holds then (µ n,1 J + µ n,2 J ) ≤ ρ max ∀n.
We proceed by induction: Assume that z n ≤ ρ max and prove that z n+1 ≤ ρ max . We have • CASE 1: z n ≤ σ We have The conclusion follows easily by (25), in particular by the fact that • CASE 2: z n > σ -CASE 2.1: z n = ρ max We have f (z n ) = f (ρ max ) = 0 and then z n+1 = z n = ρ max . -CASE 2.2: z n < ρ max We have The conclusion follows easily by (25), in particular by the fact that The two cases conclude the proof.
In the following we present the results of some numerical tests in which we compare the proposed algorithm (21)- (22) with the classical one (17)- (20). A theoretical comparison can be found in [5]. For all tests we choose ρ max = 1 and f (ρ) = ρ(1−ρ). We also assume that each arc has the same length, equal to 1, then the densities ρ 1 , ρ 2 , ρ 3 computed by the classical algorithm are defined in [0, 1] × [0, t f ], while the densities µ 1 , µ 2 computed by the proposed algorithm are defined in [0, 2] × [0, t f ]. We discretize each arc by means of 20 cells, then the junction is between the node J − 1 = 20 and the node J = 21. We choose t f = 5 and we discretize the time interval with 241 nodes. At time t = 0 the network is empty. The boundary conditions at any time t are: (2) , t) = 0.
(28) By this choice, two queues are formed behind the junction, since the outgoing road is not able to receive the car flux coming from the two incoming roads. Hereafter, we shall denote by (ρ i , ρ j ) the density on the path defined by two consecutive arcs I i , I j . In Fig. 3 we report the solution at the final time computed by the classical algorithm with equidistributed priority coefficients q 1 = q 2 = 1/2, and that computed by the proposed algorithm. For path P 1 , the solution (ρ 1 , ρ 3 ) of the classical algorithm must be compared with the solution (µ 1 , µ 1 + µ 2 ) of the proposed algorithm, while for path P 2 , the solution (ρ 2 , ρ 3 ) of the classical algorithm must be compared with the solution (µ 2 , µ 1 + µ 2 ) of the proposed algorithm. It can be seen that the two solutions are very similar, except for the fact that the junction is shifted by one cell. Indeed, the new algorithm makes the cell J = 21 play the role of the last one of the incoming roads even if it is defined as the first cell of the outgoing road. We also note that the proposed algorithm propagates the queue back in space slightly slower than the classical algorithm, probably because of the different behavior at the junction (see next test).
In Fig. 4 we show the evolution in time of the density at cells k = J − 1, J, J + 1. We note that, before the junction (k = J − 1), the two algorithms do not coincide in a transient, in particular when the queue is forming and the characteristic curves change direction (from left→right to right→left). The proposed algorithm leads to a smaller density than the classical algorithm. A cell after the junction (k = J + 1), instead, the situation is inverted: starting from the formation of the queue, the density of the proposed algorithm is slightly higher. Nevertheless, these differences are in some sense balanced, since the total number of vehicles which leave the network is exactly the same for the two algorithms, i.e.
Just after the junction (k = J), the two algorithms show a remarkably different behavior, with two "alien" values µ 1 J=21 and µ 2 J=21 which are not immediately justified in the framework of the classical theory, see [5]. These values cannot be considered by their own because they leave after the junction (where the two densities coexist). Rather, the sum of the two alien values is the value by means of which the two incoming roads perceive the presence of the junction ahead. The cell J may be seen as a region with an oversize capacity which gathers the flows  coming from the incoming roads. We may therefore say that the cell J has been augmented by a "buffer". However, our approach differs from the traffic models with buffer [12,15,17]. In fact, in those models the load of the buffer at any time is described by a dedicated function which evolves according to an additional ordinary differential equation.
Finally, in Fig. 5 we report the solution at the final time t f = 8 computed by the classical algorithm with equidistributed priority coefficients q 1 = q 2 = 1/2 and by the proposed algorithm in the case of a time-dependent left boundary conditions In this test the transient period is kept alive by the time-dependent boundary conditions. Nevertheless, we observe the same behavior as in the previous test, meaning that the different solution around the junction do not affect the solution elsewhere.    We observe the same behavior as in the previous test.

5.
One incoming road and two outgoing roads. In this section we consider a network with three arcs and one junction, with one incoming road and two outgoing roads. On the network two paths P 1 and P 2 are defined, see Fig. 6. As before, we PSfrag replacements Figure 6. A network with 3 arcs and 1 junction. Path P 1 (left) and path P 2 (right). assume that each arc has the same length, equal to 1. We denote by J the node just before the junction, see Fig. 6. Note that the nodes after the junction, namely J + 1, J + 2, etc., can refer to one path or the other one, depending on the context. We have and the scheme (22) becomes (23).
As in the previous case, we discretize each arc by means of 20 cells. The junction is between the node J = 20 and the node J + 1 = 21. At time t = 0 the network is empty. The boundary conditions for the classical algorithm are ρ 1 (0, t) = 0.5, ρ 2 (1, t) = 0, ρ 3 (1, t) = 0.9 (31) and we set the preference coefficients (5) as α 21 = 0.8 and α 31 = 0.2. By this choice, a queue is formed behind the junction, since the road 3 is almost full. According to (31), we set In Fig. 7 we report the solution at the final time t f = 11 computed by the classical  6. Two incoming roads and two outgoing roads. In this section we consider a network with four arcs and one junction, with two incoming roads and two outgoing roads. On the network four paths P 1 , P 2 , P 3 , P 4 are defined. Boundary conditions and preference coefficients are: In Fig. 8 we report the solution at the final time t f = 6 computed by the classical  algorithm and by the proposed algorithm. In this case the solutions do not coincide. Also, the total number of vehicles that leave the network is not the same (cfr. (29)), then we infer that the flux at junction must be different. As a consequence, the proposed algorithm does not maximize the flux at junction. Remarkably, the proposed algorithm gives the same density in the incoming roads, independently from the choice of the boundary conditions. The question arises which solution is chosen by the proposed algorithm. Numerical evidence shows that the multi-path algorithm coincides with the classical one in case of no formation of queues. If, instead, the constraint A(γ 1 , γ 2 ) T ∈ (Ω 3 × Ω 4 ) in (11) is active, the multi-path algorithm seems to find the solution of the problem (12), but with the additional (and not required for uniqueness of the solution) constraint (13), assuming q 1 = q 2 , see Fig. 9. . The solution of (12) at junction found by maximizing the flux (black circle), and the one which seems to be found by the multi-path algorithm (yellow circle). In case of no queues the two solutions overlap (left), in case of queue formation the two solutions are different (right).

7.
A real application. In order to test the proposed algorithm on a real case, we considered a network located in Rome (Italy), constituted by 6 two-lane roads and 7 junctions. The whole network is 328.2 km. We used the local version of the model described in Section 3.2 and we set ∆x = 100 m, ∆t = 2.5 s. We also placed four traffic lights coordinated in pairs. The final time for the simulation is t f = 3 4 h, see Fig. 10 for a screenshot. The code is written in C++ (serial) and run on an Intel i3 2.27 GHz processor. The CPU time for the entire simulation was 0.5 s. This result suggests that the proposed technique can be actually used to forecast traffic flow in large networks, keeping to a minimum the implementing effort.
8. Conclusions and future work. The outcomes of the numerical tests allow us to sketch some preliminary considerations. Avoiding any additional procedures at junctions (i.e. the most annoying part of existing algorithms), the multi-path algorithm is able to compute a solution which is admissible in the framework of the classical theory but does not provide in general the maximization of the flux at junctions. It is useful to stress here that the maximization of flux is only one of the possible conditions to ensure uniqueness of the solution. Moreover, in some cases this condition may prescribe strongly unbalanced fluxes from the incoming roads. The multi-path algorithm, instead, overrides the flux maximization favouring the flux balance.
We also stress that, once the topology of the network is suitably implemented in the function ω (total density), adding or removing a single road or path requires minimal effort. Figure 10. A screenshot of the simulator. The density is represented as "walls" along the roads.