An easy-to-use algorithm for simulating traffic flow on networks: theoretical study

In this paper we study a model for traffic flow on networks based on a hyperbolic system of conservation laws with discontinuous flux. Each equation describes the density evolution of vehicles having a common path along the network. In this formulation the junctions apparently disappear since each path is considered as a single uninterrupted road. We consider a Godunov-based approximation scheme for the system which is very easy to implement. Besides basic properties like the conservation of cars and positive bounded solutions, the scheme exhibits nice properties, being able to select automatically a reasonable solution at junctions without requiring external procedures (e.g., maximization of the flux via a linear programming method). Moreover, the scheme can be interpreted as a discretization of the traffic models with buffer, although any buffer is introduced here. Finally, we show how the scheme can be recast in the framework of the classical theory of traffic flow on networks, where a conservation law has to be solved on each arc of the network. This is achieved by solving the Riemann problem for a modified equation, and showing that its solution corresponds to the one computed by the numerical scheme.


Introduction and model setting
Starting from the introduction of the LWR model [24,26], 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,21]. 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 [17], 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 first-order version of the model proposed in [20,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, ∂ ∂t µ p (x (p) , t) + ∂ ∂x (p) µ p (x (p) , t) ω(x (p) , t) f ω(x (p) , t) = 0, x (p) ∈ P p , t > 0, 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 so-called multi-population or multi-class models, see, e.g., [2,29]. 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 [22]). In that case, vehicles are divided in different populations on the basis on their source-destination pair. Given the total density ω of all vehicles, the density µ p of the p-th population is given by µ p (x, t) = π p x, t, O(p), D(p) ω(x, t), 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 [23] 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 [25], 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 [23]. We also point out the paper [19], where a numerical method for (systems of) scalar conservation laws with discontinuous flux is proposed, and the paper [28], where the convergence of a Godunov-based scheme for scalar conservation laws with discontinuous flux is investigated.
Goal. In the "twin" paper [5], the authors propose a simple Godunov-based numerical scheme for the system (1),(3). In the same paper the authors performed a numerical comparison between that scheme and the algorithm described in [4,14] based again on the LWR model and the Godunov scheme, showing that the two numerical solutions do not coincide in general.
In this paper we want to investigate theoretically the properties of the solution of the proposed scheme. This is a challenging task because it is well known that the system is ill-posed and infinite solutions are admissible. The question arises which solution is automatically chosen by the numerical scheme if no additional constraints are given (typically the maximization of the flux at junctions). We recast the problem in the framework of the classical theory of traffic flow on networks [14], and we show a modified model (together with a Riemann problem for the associated equation) which is correctly solved by the numerical scheme without ambiguity. To define the modified equation, a new flux h, which depends on f , is used near the junctions.
Paper organization. In Section 1 we recall some basic facts about the LWR model on networks, the Godunov approximation and the notions of Riemann and half-Riemann problem. The main references are [4,8,14,25].
In Section 2 we introduce a Godunov-based numerical scheme to solve the multi-path model (1),(3) and a local version of the model suitable for large networks.
In Section 3 we investigate in detail the case of a merge, i.e. a junction with two incoming roads and one outgoing road. In particular we recast the problem in the framework of the classical theory of traffic flow on networks [14] and we exhibit a modified problem which is correctly solved by the proposed numerical scheme. We also investigate two generalizations of the model: the case of different flux functions and the case with N incoming roads.
In Section 4 we investigate the case of a diverge, i.e. a junction with one incoming road and two outgoing roads.
In Section 5 we sketch some conclusions and in the Appendix we collect the technical proofs.

Background
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,25].

LWR on networks
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, and we denote by {I i := [a i , b i ]}, i = 1, . . . , n the incoming roads and by {I j := [a j , b j ]}, j = n + 1, . . . , n + m the 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 (γ 1 , . . . , γ n ) ∈ {qs, s ∈ R + }, 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,

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 f -based Godunov numerical flux G f is defined as if ρ ℓ > ρ r and ρ r > σ.
In the following we drop the subscript f from G whenever the underlying flux is deducible without ambiguity. 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

Riemann and half-Riemann problem
Let us first introduce the following notation, which will be used through the paper.

The Riemann problem
Let us consider the Riemann problem for some constant initial data ρ ℓ and ρ r in [0, ρ max ]. The unique weak entropy solution is given by: • If ρ ℓ < ρ r (shock wave): We have a shock wave with positive speed where the function ψ(ξ) is defined by the solution of f ′ (ψ(ξ)) = ξ. We have a rarefaction wave with

Half Riemann problem
Following [25], we call half Riemann problem the simple case of an initial-boundary value problem in the quarter of plane {x ≤ 0} or {x ≥ 0} when the initial condition is a constant. The problem is then to find the acceptable boundary condition at x = 0.
The study of the left-half problem is equivalent to searching an artificial right state in the problem (22) which lead to waves with negative speed, in order to know which are the states attainable along the line x = 0. Definition 1.2. For any ρ ℓ ∈ [0, ρ max ], we define N (ρ ℓ ) to be the set of pointsρ ∈ [0, ρ max ]\{ρ ℓ } such that the solution to the Riemann problem   contains only waves with negative speed.
It is trivial to verify the following result.
The study of the right-half problem is equivalent to searching an artificial left state in the problem (22) which lead to waves with positive speed, in order to know which are the states attainable along the line x = 0. Definition 1.4. For any ρ r ∈ [0, ρ max ], we define P (ρ r ) to be the set of pointsρ ∈ [0, ρ max ]\{ρ r } such that the solution to the Riemann problem   contains only waves with positive speed.
Analogously to the previous case, we get the following result.

The (half ) Riemann problem with different fluxes
To our purposes, we need a trivial generalization of the Riemann problem to different fluxes. Let us consider the following problem where h ℓ , h r : [0, ρ max ] → R are two C 1 fluxes. This problem is different from what is commonly called "conservation law with discontinuous flux", since in our case the discontinuity does not depend explicitly on x.
Instead, it depends on the solution ρ itself and it is located at the boundary of two adjacent advection problems of a constant state. Let us clarify this point by computing the entropy solution to (29).
, characteristic lines starting from the line {x < 0} meet those starting from the line {x > 0}, creating a shock x = ξ(t), see Fig. 2a. In order to have the mass conserved across the shock, a Rankine-Hugoniot-like condition must be verified. This condition is easy found by standard arguments [16, Chapt.77] asξ Then we have a shock wave with positive speed ifξ > 0, and with negative speed ifξ < 0.
a rarefaction is created, see Fig. 2b. As in the classical Riemann problem, the "boundaries" of the rarefaction are given by the lines Note that the sign of the waves can be computed even if the solution inside the rarefaction is not explicitly given. Following Definitions 1.2 and 1.4 we introduce the following two definitions: contains only waves with negative speed.

Numerical approximation of the multi-path model on small and large networks
In this section we present the numerical approximation we use to discretize the system (1),(3). Same approximation is also considered in the "twin" paper [5]. After that, we discuss how the multi-path model can be modified to deal easily with large networks.

Numerical approximation
Let us denote by µ n,p k (q) the approximate density µ p (x (q) k , t n ), where k (q) 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 to (3) at any internal nodes k p as for n ≥ 0 and p = 1, . . . , N P . Note the intrinsic asymmetry of this conservative scheme: 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.
Considering that the system (3) is ill-posed and the scheme (34) itself does not impose any additional constraints to the solution, the question arises which behavior for the density is described by the scheme, especially at junctions. We try to answer to this question is Section 3.

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 allows a large number of paths, the computation of ω n k p become 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 1. 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 (34) 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.
The original (global) multi-path approach of Section 2.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 (cfr. [10,13,27] on this point). To better understand the difference, let us consider the network in Fig. 3. 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 the density in [a 3 , b 3 ]. Clearly, only the global algorithm gives the correct solution.
We also note that the classical algorithm described in Section 1 is local, while the one presented in [13] is global.
J -1 Figure 4. A network with 3 arcs and 1 junction, representing a merge. Path P 1 (left) and path P 2 (right). and the scheme (34) becomes 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 3.1. Let the initial densities around the junction be admissible, namely If the following CFL-like condition holds Proof. Let us first prove that (37) implies Let us define M := sup ρ∈(0,ρmax) |f ′ (ρ)|. By (37) we have Noting that f (ρ max ) = 0, and using the Lagrange theorem and (39), we have To simplify the notations, let us introduce the auxiliary variable z n := µ n,1 J + µ n,2 J . The worst case happens when µ n,1 J−1 = µ n,2 J−1 = σ (incoming roads try to transfer the maximal flux to cell J) and µ n,1 J+1 + µ n,2 J+1 = ρ max (no flux from cell J to cell J + 1). In this case the equation for z is We proceed by induction: Assume that z n ≤ ρ max and prove that z n+1 ≤ ρ max . We have The conclusion follows easily by (38), in particular by the fact that We have The conclusion follows easily by (38), in particular by the fact that

Analogies and differences with respect to the classical algorithm
Here we present two brief examples which show the analogies and the differences between the classical scheme (17)- (20) and the proposed scheme (34).
Example 3.2. Assume that, at some time step n, we have This choice guarantees that the outgoing road can receive the flux of cars coming from both the two incoming roads, then no queue is formed. The solution of the maximization problem (14)-(15) is Then, after one time step, the two incoming roads transfer a density equal to ∆t ∆x (f (u ℓ ) + f (v ℓ )) to the outgoing road, see (19)- (20). Now consider (36) at J − 1. The outgoing density along the first path is ∆t ∆x Analogously, the outgoing density along the second path is ∆t ∆x f (v ℓ ). Considering the cell J, we see that the incoming density along the first path is ∆t ∆x and the incoming density along the second path is ∆t ∆x f (v ℓ ). Then, the two incoming roads transfer a density equal to ∆t ∆x (f (u ℓ ) + f (v ℓ )) to the outgoing road and the two algorithms coincide. Example 3.3. Assume that, at some time step n, we have In this case the outgoing road cannot receive the flux of cars coming from both the two incoming roads, then a queue is formed. The solution of the maximization problem (14)- (15), with priority coefficients q 1 and q 2 , is (we do not consider here the case such that a further projection onto Ω is needed, see Section 1.1). Then, after one time step, the two incoming road transfer a density equal to ∆t ∆x f (σ) to the outgoing road, see (19)- (20). As in the Example 3.2, consider (36) at J − 1. The outgoing density along the first path is ∆t ∆x f (σ) and the outgoing density along the second path is ∆t ∆x f (σ). Considering the cell J, we see that the incoming density along the first path is ∆t ∆x f (σ) and the incoming density along the second path is ∆t ∆x f (σ). In conclusion, the total density transferred at junction is ∆t ∆x 2f (σ) and then the two algorithms do not coincide.

The modified equation
In Section 3.1 we have shown that the scheme (34) is not always compatible with the classical theory, the latter imposing the vector of the incoming fluxes (γ 1 , . . . , γ n ) to be in Ω, see (11). Then the question arises which analytical problem the proposed scheme solves in the framework of the classical theory. In analogy with the "modified equation" in the sense of LeVeque [23,Sect. 11.1], in this section we exhibit a new problem which is correctly solved by the new approach. We start with a useful technical Lemma.
Remark 3.5. We point out that the quantities G(ρ ℓ , σ) and G(σ, ρ r ) appearing in (42)-(43) match exactly the definitions of the maximum incoming flux (7) and of the maximum outgoing flux (8), respectively, i.e. for all incoming roads and, for all outgoing roads, They also correspond to the demand and supply functions used in [18,22]. Consider again the system (35)-(36). After the junction, it is convenient dealing with the sum ω n of the two densities µ n,1 and µ n,2 rather than the two densities separately. Then, across the junction we have the following equations: where γ * i , γ * j are defined in (14)-(15) (with the obvious correspondence between the path index p and the road index i, j). Focusing on node J, the total left flux is given by and the total right flux is given by Interestingly, from (48) and the properties of G (Lemma 3.4 and Remark 3.5) we get that the flux is actually maximized on every single path, and not on the whole junction. We will come back on this point in Section 5.
We are now ready to introduce a first version of the modified problem, which is valid for the numerical scheme (35)-(36).
(2) The flux function f r is an admissible flux for the bottleneck problem with one incoming road with flux h(·) = f (·) + f (σ) and one outgoing road with flux f (·). More precisely, by applying the relations (44)-(45), we have . Proof. Let us prove the two cases separately.
(1) (Left flux) We have to prove that Let us now introduce the modified equation. To this end, it is convenient to slightly reformulate the problem assuming roads to be half lines and changing notations. We denote by u( with In the next section we show that, for a particular choice of C, the solution to the Riemann problem associated to equations (50)-(51) equals the solution of the numerical scheme around the junction, namely the scheme described in Section 2.1 is consistent with equations (50)-(51).
Finally note that the flux h defined in Proposition 3.6 is just an upper bound of the flux h defined in (51).
Remark 3.7. (Analogies and differences with respect to traffic models with buffer). As shown in Fig. 5, the road at cell J may be seen as an area with an oversize capacity which gathers the flows coming from the incoming roads. We may therefore say that the cell J acts as a "buffer", i.e. it is a region of size ∆x used to temporarily store cars densities while there are being moved from the incoming roads to the outgoing road. In [12,15,18] traffic models with buffer are proposed. In such models, the buffer is 0-dimensional, and its load is described by the number of cars r(t) lying at time t in it. The load r varies accordingly to the difference between the inflow and the outflow at the buffer and evolves according to an ODE. It is also assumed that the maximum number of cars which can enter or exit the buffer per unit of time is a constant parameter (µ in the original notation). In the same spirit of [12,18], in our modeling approach the function ω n J may be assumed to describe the evolution of cars densities in the buffer, which is here represented by an interval of size ∆x. Using Lemma 3.4, equation (47) can be seen as the first-order Euler approximation of an ODE for a semi-discrete function ω J (t) Note that the function ω J (t) is related to the number of cars r(t) by the relation r(t) = ω J (t)∆x.

The Riemann problem at junction and asymptotic solution
In this section we study the Riemann problem associated to the system (50)-(51). Initial data are given by the constants (53) Generally, this definition does not select a unique solution to the problem (50)-(52). To fix this, we select particular fluxes at x = 0, ∆x somehow "suggested" by the numerical scheme, and we study the corresponding evolution of the density. Doing this, we verify that the Riemann problem is not degenerate (waves do not rebound forever) and we compute the constant asymptotic solution for t → +∞. In Section 3.4 we will show that the numerical scheme proposed in Section 2.1 computes that asymptotic solution.
From now on, we assume that We also consider in the plane (f (u ℓ ), f (v ℓ )) the four regions (A),(B),(B ′ ), and (C) depicted in Fig. 6 and characterized by the relations For symmetry reasons, we will restrict our attention only to regions (A), (B), and (C). Note that the case (A) corresponds to no formation of queues, the case (B) corresponds to the formation of a queue along P 1 , and the case (C) corresponds to the formation of two queues along P 1 and P 2 . We also assume that at the initial time (or, equivalently, C(x, 0) ≡ f (σ)), meaning that the maximal capacity of the junction [0, ∆x] doubles the maximal capacity of the other roads. In addition, we assume that Assumption (58) translates the fact that at the initial time the junction ensures maximal flow, equal to 2f (σ). Therefore, the junction area itself is not responsible for a possible congestion, while, if any, this will be due to the choice of w r . Assumptions (55) and (58) could be relaxed and are taken mainly to ease the proofs of the next theorems 3.9 and 3.10, where the continuous and the discrete problem are analyzed and compared. Indeed, we still cover all the interesting cases (free and congested state, queue formation) while limiting the number of cases to be studied, otherwise astronomical.
where G · (·, ·) is the Godunov flux defined in (18) and h has the form (51) for some particular function C (to be defined). For every initial data u ℓ , v ℓ , z c , w r satisfying (55) and (58) there exists one solution (in the sense of Definition 3.8) of the problem. Moreover, by assuming the initial data to be in the regions (A), (B), (C) defined in (56), the following constant quadruplets (ū,v,z,w) are, respectively, asymptotic stationary solutions of the problem.
The proof, together with the exact expression of C, is given in the Appendix.

The numerical scheme at junction and asymptotic approximate solution
Let us focus on the behaviour of the scheme (35)-(36) at the junction, considering only the two nodes J − 1 and the node J, see (46)-(47). We set the following boundary conditions: with w r := w 1 r + w 2 r . To be coherent with the theoretical investigation of Section 3.3, we assume (55) still holds (note that here z c does not appear). The next result concerns the stationary solutions of the scheme.
is determined by the following conditions: Remark 3.12. Borderline cases are relatively easy to study. The case w r = σ is essentially analogous to the case w r > σ and all the results described above are still valid. Instead, a more complicated situation arises by assuming that the boundary conditions u ℓ , v ℓ , and w r = w 1 r + w 2 r are chosen in such a way that the point (f (u ℓ ), f (v ℓ )) lies at the boundaries of all four regions depicted in Fig. 6. In that case, infinite stationary points are present. Indeed, if (55) holds true and we necessarily have u ℓ = v ℓ and all the points (μ 1 J−1 ,μ 2 J−1 ,μ 1 J ,μ 2 J ) such that are stationary, but not in general stable.
Remark 3.13. The main result of the paper comes from the comparison between the two theorems 3.9 and 3.10.
Clearly we cannot compare exactly the theoretical and numerical solution because the first one is composed by a quadruplet (ū,v,z,w) while the second one by a triplet (μ 1 J−1 ,μ 2 J−1 ,ω J =μ 1 J +μ 2 J ). The reason for this is mainly technical, i.e. we chosen to have only three unknowns in the numerical setting to make the proof doable. Anyway, in the case (A) we have the correspondence (ū,v,z =w) = (μ 1 J−1 ,μ 2 J−1 ,ω J ), and in the case (B),(C) we have the correspondence (ū,v,z) = (μ 1 J−1 ,μ 2 J−1 ,ω J ). Remark 3.14. As already noted in [5], when queues are formed along both incoming roads, the density split in two equal values, regardless the ratio between the initial values u ℓ and v ℓ . In the classical theory, this effect can be achieved by dropping the maximization of the flux and keeping the priorities coefficients defined in (13), setting them to q 1 = q 2 = 1 2 .

Numerical tests
In this section we present three numerical tests, in order to confirm experimentally the results described in Sections 3.3 and 3.4. A more complete numerical study of the multi-path model presented here can be found in [5]. We assume that each arc has the same length, equal to 1, then each path has length equal to 2. We also assume that the flux has the classical form f (ρ) := ρ(1 − ρ), so that ρ max = 1 and σ = 0.5. We divide each arc in 25 nodes (then the junction is found at J = 26 along each path) and we impose Dirichlet boundary conditions at the beginning and the end of each path, i.e. at points 0 (1) , 0 (2) , 2 (1) , and 2 (2) , see Fig. 4. Note that the points 2 (1) and 2 (2) correspond to the same physical point. Mutatis mutandis, these boundary conditions have the same role of u ℓ , v ℓ , w 1 r , and w 2 r in (61). We assume roads are empty at the initial time and we look for the stationary solution obtained for t → ∞. Test 1. We consider the case (A) in (56) (no formation of queue along the incoming roads). For any t > 0 we set µ 1 (0 (1) , t) = 0.1, µ 2 (0 (2) , t) = 0.15, µ 1 (2 (1) , t) = 0.3, µ 2 (2 (2) , t) = 0.3.
In Fig. 8 we report the stationary solution for µ 1 along P 1 , µ 2 along P 2 and the total density ω along P 1 .
Here  Figure 7. Test 1. Stationary solution for µ 1 along P 1 , µ 2 along P 2 and the total density ω along P 1 .  Figure 8. Test 2. Stationary solution for µ 1 along P 1 , µ 2 along P 2 and the total density ω along P 1 . corresponds toū =z in Theorem 3.9 and toμ 1 J−1 =ω J in Theorem 3.10. Looking at the single densities, at the node J we find two isolated values equal to corresponding respectively toμ 1 J andμ 2 J in Theorem 3. 10. Surprisingly, if we look at the total density ω we see that the node J plays the role of the last node of the incoming roads even if it defined as the first node of the outgoing road. In other words, the scheme shifts the junction to one cell to the right. This happens also in the case (C) but does not happen in the case (A). Test 3. We consider the case (C) in (56) (formation of two queues along the incoming roads). For any t > 0 we set In Fig. 9 we report the stationary solution for µ 1 along P 1 , µ 2 along P 2 and the total density ω along P 1 . Here,  Figure 9. Test 3. Stationary solution for µ 1 along P 1 , µ 2 along P 2 and the total density ω along P 1 .
the right boundary conditions affect the solution, but only by means of their sum 0.8 = 0.3 + 0.5, while the left boundary conditions have no effect. The densities split equally along the two paths, regardless the fact that the boundary conditions are not equal. We refer the reader to [5] for a comparison with the classical theory and the relationship with priority coefficients q. Before the junction, the densities are equal toλ ≈ 0.9123 > σ with f (λ) = f (0.8) 2 . The valueλ corresponds toū =v in Theorem 3.9 and toμ 1 J−1 =μ 2 J−1 in Theorem 3.10. At the node J we find an isolated value equal toλ 2 , corresponding toμ 1 J =μ 2 J in Theorem 3.10. We refer again to [5] for a numerical study of the density evolution at this node.

Extension to roads with different fluxes
We consider here the case such that the three arcs have different flux functions, i.e. we have three triples (f 1 (ρ), ρ 1 max , σ 1 ), (f 2 (ρ), ρ 2 max , σ 2 ), (f 3 (ρ), ρ 3 max , σ 3 ), all of them satisfying assumptions (4). In order to deal with the new problem, the scheme (34) must be generalized. Indeed, vehicles following path P p now perceive the change of flux at junction. As in the classical theory, we simply introduce, in each path, a trivial junction with one incoming and one outgoing road. More precisely, the first path passes from f 1 to f 3 while the second path passes from f 2 to f 3 .

Extension to merges with N incoming roads
We consider here the case of a junction with N incoming roads and one outgoing road. All the results presented above extend easily. In particular the CFL-like condition (37) becomes N ∆t ∆x sup

Conclusions
In this paper we investigated the properties of a first-order macroscopic model for traffic flow on networks where vehicles are divided on the basis of their desired path. The resulting Godunov-based numerical scheme can be implemented in minutes since it requires no additional procedures to manage the solution at junctions. Numerical and computational details have been left to the "twin" paper [5].
Despite the simplicity of the model and of the numerical discretization, the method shows many interesting properties, some of which do not require ad hoc ingredients necessary in other models. In the following we summarize the main features of our approach.
• The model selects automatically a solution at junctions that maximizes the flow from an user point of view: Formulas (48), coupled with Lemma 3.4 and Remark 3.5, state indeed that on each single path the density is maximized (user optimum). Therefore, the scheme does not compute in general the maximal flow that could possibly be transferred over the node (global optimum), as it happens in more standard approaches. See [27,Sect. 3.1] for a discussion on this point. • In the challenging case of a merge, we proved that the proposed numerical scheme is consistent with the Riemann problem (50)-(52). As a consequence, the numerical approximation automatically assigns to the junction a finite spatial dimension and its density evolution is managed by a special flux function h, as defined in (51). Specifically, the function h balances the incoming flows and the outgoing flow according to the fact that drivers want to maximize their own flow. It is able to widen temporarily the spaced-junction capacity, allowing the passage of vehicles and making the junction act as a "buffer", see Remark 3.7. The exact expression of the space-time-dependent flux h can be found in the proof of Theorem 3.9 (see Appendix A.1). • The multi-path model, together with the proposed discretization, fulfills all the 7 requirements for junction models proposed in [27,Sect. 3.1]. In particular, it is generally applicable irrespectively of the number of incoming and outgoing roads, the traffic never flows backwards and all flows are non-negative, vehicles are conserved, and turning fractions are preserved. • Finally we have shown that the standard CFL condition (16) is not in general sufficient to ensure that the solution of the numerical scheme (34) is admissible at junctions. When more than one queue is formed, it is required a stronger condition which depends on the number of incoming roads at the junction. This condition is given, in the case of a merge with N incoming roads, by (67). and we look for the fluxesĥ(·),h(·) and the constantsũ,ṽ,ẑ,z,ŵ (see Fig. 11) such that u(x=0−, t=0+) =ũ, v(0−, 0+) =ṽ, z(0+, 0+) =ẑ, z(∆x−, 0+) =z, w(∆x+, 0+) =ŵ with (ŵ = w r orŵ ∈ P (w r )) and f (ŵ) =h(z), where the sets N , P , N , P are defined in Section 1.3. Figure 11. Solution of the Riemann problem in the plane (x, t) just after the initial time.
We claim that the solution is: It is important to note that such aẑ actually exists since Let us now verify that the fluxes/constants in (78) are actually solution of (72)-(76). (72) We haveũ = u ℓ and (73)ṽ is analogous.
Summarizing, we have that the waves L u 1 = (u ℓ ,ũ), L v 1 = (v ℓ ,ṽ), and L 4 = (ŵ, w r ) in Fig. 11 are not even created, while the wave L 2 is a strictly positive shock and the wave L 3 is a strictly negative shock. When L 2 meets L 3 , a new wave L 5 := (ĥ,ẑ), (h,z) arises, see Fig. 12. Let us prove that L 5 is a strictly positive shock. We have and then, recalling Sect. 1.3.3, L 5 is a shock. Moreover, and then, by (30) andz >ẑ, L 5 is strictly positive.
• CASE (B) Let us start again from the initial time t = 0 and find the fluxesĥ(·),h(·) and the constantsũ,ṽ,ẑ,z,ŵ such that (72)-(76) hold true. We claim that the solution is (83) and then such az exists.
(74) We prove that (ĥ,ẑ) ∈ P(h c , z c ). To this end, we show that the wave L 2 is a strictly positive shock, see Fig. 13. We have and then, recalling Sect. 1.3.3, L 2 is a shock. Moreover, and then, by (30) and z c >ẑ, L 2 is strictly positive. Finally, the equalityĥ(ẑ) = f (ũ) + f (ṽ) is trivially verified. (75) We prove that (h,z) ∈ N (h c , z c ). To this end, we show that L 3 is a strictly negative shock, see Fig. 13. We have and then, recalling Sect. 1.3.3, L 3 is a shock. Moreover, and then, by (30) andz > z c , L 3 is strictly negative. Finally, we havẽ (76) We haveŵ = w r and f (ŵ) =h(z). Summarizing, the situation is similar to the case (A), see Fig. 13: L u 1 , L v 1 , and L 4 are not created. L 2 is a strictly positive shock and L 3 is a strictly negative shock. When L 2 meets L 3 , a new wave L 5 = (ĥ,ẑ), (h,z) arises. Let us prove that L 5 is a strictly negative shock.
Summarizing, we have that the wave L 6 is not created. Let us study the waves L u 7 and L v 7 . L u 7 is a strictly negative shock sinceũ L v 7 instead is not created since (ṽ,ṽ 2 ) = (v ℓ , v ℓ ). The proof is concluded by choosinḡ u =ũ 2 ,v =ṽ 2 ,z =z,w =ŵ.
We can also assume that because if this does not hold true, we come back to CASE A.1. We have (55) + (105) + (106) ⇒ v ℓ < y and f (y) < f (v ℓ ) and then y > v ♯ ℓ > σ.
As in CASE A.4, we can also assume that and we get analogously x > u ♯ ℓ > σ and y > v ♯ ℓ > σ.
The quantity G(z, w r ) can be in principle equal to f (z), f (w r ), or f (σ). But (55) ⇒ G(z, w r ) = f (σ). By contradiction, we can also prove that G(z, w r ) = f (z). Indeed, if G(z, w r ) = f (z) we would have z < w ♯ r < σ (120) and f (σ) which contradicts (117). As a consequence, we must have G(z, w r ) = f (w r ), and then, by (98) and (117), Moreover, we have that z > σ, since, if z ≤ σ, we get a contradiction as before in (121). Finally, we have and, by (122) and (124), This is indeed a stationary point for the system.
Proof. Consider again the discrete dynamical system (DS) and the auxiliary variables (95).
• CASE (A) In a neighbourhood of the stationary point (62) the system has the form          x n+1 = F 1 (x n , y n , z n ) := x n − ∆t ∆x f (x n ) − f (u ℓ ) y n+1 = F 2 (x n , y n , z n ) := y n − ∆t ∆x f (y n ) − f (v ℓ ) z n+1 = F 3 (x n , y n , z n ) := z n − ∆t ∆x f (z n ) − f (x n ) − f (y n )