FAST ALGORITHMS FOR THE APPROXIMATION OF A TRAFFIC FLOW MODEL ON NETWORKS DIIMA

New computation algorithms for a fluid-dynamic mathematical model of flows on networks are proposed, described and tested. First we improve the classical Godunov scheme (G) for a special flux function, thus obtaining a more efficient method, the Fast Godunov scheme (FG) which reduces the number of evaluations for the numerical flux. Then a new method, namely the Fast Shock Fitting method (FSF), based on good theorical properties of the solution of the problem is introduced. Numerical results and efficiency tests are presented in order to show the behaviour of FSF in comparison with G, FG and a conservative scheme of second order.

1. Introduction.In this paper we introduce new computation algorithms to deal with a fluid-dynamic mathematical model for flows on networks.In particular, we consider traffic and telecommunication networks.
More precisely, we deal with the models proposed in [6,8], which are macroscopic fluid-dynamic models, respectively, for traffic flow on road networks and for flows of information on a telecommunication network.In the 1950s James Lighthill and Gerald Whitham in [14], and independently Richards in [18], proposed to apply fluid-dynamics concepts to traffic.This nonlinear model is based on the single conservation law expressing the conservation of cars: where ρ = ρ(t, x) is the density, with ρ ∈ [0, ρ max ], (t, x) ∈ R 2 and ρ max is the maximum density; f (ρ) is the flux, which can be written as f (ρ) = ρv and v is the average velocity of cars.Typically v is assumed to be a smooth decreasing function of the density ρ and f is assumed concave.Fluid-dynamic models can describe macroscopic phenomena as shocks formation and propagation.Since they can develop discontinuities in finite time even starting from smooth initial data, the study of the analytical and numerical aspects is fundamental.Recently, in [5,6,8], some models based on the conservation formulation (1.1) have been proposed for flows on networks and it was proved existence of solutions to Cauchy problems.Road (transmission lines) networks are a finite set of roads (lines) linked by junctions and, in such models, some rules (see [8]) are introduced to uniquely solve Riemann problems at junctions.The numerical approximation of the traffic models described in [5,6] has been studied in [4], where simulation algorithms based on some numerical schemes such as Godunov scheme and 3-velocities Kinetic schemes of first and second order have been proposed, with the introduction of suitable boundary conditions at the junctions.
Motivated by real applications we address the following issues demanding for fast algorithms: • Real urban networks are usually formed by thousands of arcs, representing roads, walking paths etc., and nodes, representing junctions of various types.• In case of accidents or other reasons for partial interruption of part of the network, it is important to have a real-time simulation tool in order to study the evolution of heavy traffic load and planning a strategy for traffic redirection.• The future aim is to have an important ability of communication with drivers (via information panels, GPS system etc.)This open the possibility of optimization via traffic regulation even in real time.Clearly any optimization algorithm needs some computational effort, whose basic step is the solution of an evolution equation for the flow on the network.
First, starting from the numerical approximation methods described in [4], we improve the performance of the Godunov scheme using a simplified flux function.Indeed, by considering only the possible values that the numerical flux correspondent to the analytical flux (2.2) can assume, we introduce the Fast Godunov scheme (FG).The latter can be implemented substituting the usual computations for numerical flux with some "if then" instructions based on a-priori given flux values.This clearly strongly speeds up computations.Then, supported by theoretical considerations (cfr.Section 2), we define a new numerical method, which is based on the existence of a single separating shock.Indeed, such shock separates two phases, one of light traffic and one of heavy traffic.More precisely, the uniqueness of the separating shock is proved for solutions starting from an empty network, but it does not hold for an arbitrary initial data.Focusing again on the analytical flux (2.2), an exact procedure to approximate solutions is introduced: the Fast Shock Fitting (FSF) method, which is a fast algorithm that constructs solutions by localizing shocks.We exploit the existence of a unique separating shock, that is in fact a generalised characteristic, by tracking it with an exact fitting procedure.If we start from a not identically zero density, which was anyway obtained loading the network from an empty one, it is possible first to apply an approximate fitting procedure and then use our algorithm.Once the separating shock is tracked, thanks to the special flux function considered, the rest of the solution is easily computed by a shift procedure.
Notice that FG is applicable also for general initial configurations or in the case of conservation laws with a source term, while FSF cannot be used for these problems.
To complete the analysis, we also consider the 3-velocities Kinetic scheme of second order (K3V) for the same flux function.Then we compare through some numerical experiments the performances of the schemes regarding the accuracy and the time of execution.We observe that FG allows to save more than 50% of CPU time with respect to the classical Godunov.On the other hand, FSF outperforms FG of an increasing percentage, which is around 70% for large networks and small space discretization steps.
The paper is organized as follows.Section 2 is devoted to the description of the model.In Section 3 we describe the numerical schemes to approximate the problem, namely FG, K3V and FSF.In Section 4 we present some numerical experiments and analyse the features of the new approximation methods in terms of accuracy and CPU time.
2. Analytical framework.We consider the equation (1.1) on a traffic or a telecommunication network.From now on we will focus only on road networks, since the case of telecommunication networks, where roads are replaced by transmission lines and the role of junctions is played by routers, presents the same features.The parametrization of arcs composing a network is made through a set of intervals . ., N , with the endpoints possibly infinite.The datum is given by a finite set of functions ρ i defined on [0, +∞[×I i .Each road can be incoming at most for one junction and outgoing for at most one junction, hence the complete model is given by a pair (I, J ), with I = {I i : i = 1, . . ., N } the collection of arcs and J the set of junctions.Let us consider the conservation equation: If the velocity function is chosen as follows: then the flux f = vρ is If we fix ρ max = 1, v = 1 and σ = 1 2 , then the graph of the velocity is depicted in Fig. 1, while flux function assumes the shape described by Fig. 2. Notice that the velocities of Lax waves not crossing σ is v or −v.
Remark 1.Let ρ = (ρ 1 , . . ., ρ n+m ) be a weak solution at the junction such that each x → ρ i (t, x) has bounded variation.We can deduce that ρ satisfies the conservation of fluxes at the junction J, for almost every t > 0, namely A Riemann problem is a Cauchy problem with Heaviside type initial data.In case of concave or convex fluxes, the Riemann solutions are of two types: continuous waves called rarefactions and travelling discontinuities called shocks.For a junction a Riemann problem is a Cauchy problem with an initial data constant on each road (the natural discontinuity being the junction).We look for self-similar solutions with one wave on each road.Definition 1.A Riemann solver for the junction J is a map RS : associates to Riemann data ρ 0 = (ρ 10 , . . ., ρ n+m0 ) at J a vector ρ = (ρ 1 , . . ., ρn+m ) so that the solution on an incoming transmission line I i , i = 1, . . ., n is given by the wave (ρ i , ρi ) and on an outgoing one I j , j = n+1, . . ., n+m, is given by the wave (ρ j , ρ j ).We require the consistency condition (CC): RS(RS(ρ 0 )) = RS(ρ 0 ).
Since condition (2.4) does not prescribe the uniqueness of solutions, it is necessary to impose further conditions.Here we refer to two possible strategies that can be followed to solve the problem at the junction.In particular, we deal with the Riemann solver at a junction introduced in [6] and represented by the criteria: RS1: The RS is described by two rules: (A): there are some fixed coefficients, the prescribed preferences of drivers, that express the distribution of traffic from incoming to outgoing roads; (B): respecting (A), drivers choices are made in order to maximize the flux.
Another rule to solve the Riemann problem at a junction was introduced in [8]: RS2: The flow through the junction is maximized both over incoming and outgoing arcs.Focusing on RS1, to deal with rule (A) we need to introduce the traffic distribution matrix, whose elements 0 < α ji < 1 describe how traffic coming from incoming roads i distributes through outgoing roads j.This rule implies condition (2.4).Rule (B) gives a criterion to select a unique solution, by maximizing the flux at junctions.Once the solutions to Riemann problems are given, one can use a wave front tracking algorithm to construct a sequence of approximate solutions.
Definition 2. Fix an approximate wave front tracking solution ρ and a road + for concave flux.Definition 3. Fix an approximate wave front tracking solution ρ and a junction J.We say that an incoming road I i has a good datum at J at time t > 0 if and a bad datum otherwise.We say that an outgoing road I i has a good datum at and a bad datum otherwise.
For the proof of the next Lemma, we need to introduce some more notation: Obviously τ is well defined and satisfies: Lemma 1.If road I i of a junction J has a good datum, then it remains good after interactions with J of waves coming from other roads.Then, no big shock can be produced in this way.If road I i has a bad datum, then after interactions with J of waves coming from other roads, either the datum of I i is unchanged or a big shock is produced (and the new datum is good).
Proof.Let us consider the case in which I i is an incoming road, depicted in Fig. 3.In order to enter from the junction J, a wave must have a negative speed.Therefore, if , thus the resulting wave can be a rarefaction or a shock.If the initial data , hence, if a wave is created, it is a shock.In the case of outgoing roads, represented in Fig. 4, the wave must have a positive speed.Thus, if ρ i (t, a i +) ∈ [0, σ] (good datum), then the new state ρi ∈ [0, σ], thus the resulting wave can be a shock or a rarefaction.If   Now we focus on solutions generated starting by an empty network, i.e. ρ(0, x) = 0, to give some theoretical results on which the FSF scheme is based.We introduce below the definition of separating shock, which is in fact a generalised characteristic, which separates each road in two zones of light and heavy traffic.Roughly speaking, a generalised characteristic is a curve which solves the classical equation for characteristic as long as the solution is smooth, while it follows a shock when some discontinuity appears.For precise definition and theoretical results we refer to [7].
Definition 5. Consider a road network (I, J ) and an approximate wave front tracking solution ρ, such that ρ(0, x) = 0.For every road I i , we define S i (t) to be the generalised characteristic such that S i (0) = b i , with the following conventions.If S i (t) = a i (respectively b i ), then S i ≡ a i (respectively b i ) up to the first time in which a generalised characteristic starts from a i (respectively b i ).Moreover, we assume the characteristic velocity to be equal to zero at σ.We call S i the separating shock of I i .Definition 6.We say that there is a latent separation shock on road I i at a i (on the left) if S i (t) = a i .Analogously, there is a latent separation shock on road Lemma 2. For every solution ρ with ρ(0, x) = 0 and for every I i the following statements hold: a) For every t ≥ 0 such that S i (t) ∈]a i , b i [, there exists a (possibly null) shock

¥@ ¥@
A ¥A B ¥B Proof.We prove both claims at the same time by recursion.Let t be the first time such that S i (t) < b i for t in a right neighbourhood of t.Then both claims are easily verified up to time t.Then, from Lemma 1, we have that for t in a neighbourhood of t, S i (t) is a big The main idea for FSF is the following.Since we have at most one separation shock per arc, the procedure of the algorithm for determining exactly the solutions consists of the two main steps below: • exact localization of the separating shock; • shift by velocities v or −v.
3. Approximation schemes.In order to find approximate solutions, we use two new methods, namely a modified version of the classical Godunov scheme (FG) and the Fast Shock Fitting (FSF) scheme, which solves exactly the problem.These schemes will be compared to an efficient scheme, the 3-Velocities Kinetic scheme of second order (K3V), already presented and discussed in [4].In particular, FG brings a simplification in the computations of the numerical flux respect to the classical scheme.In fact, by analysing the expression of the flux function (2.2) we are able to restrict the range of possible values which can be assumed by the flux itself (see paragraph 3.1).FSF represents a completely original method, which is based on a localization technique of separating shocks and allows us to compute the exact solution, see paragraph 3.2.Concerning the discrete kinetic scheme, we recall that is a quite recent scheme for conservation laws [15,1], adapted to traffic flow problem in [4].The kinetic scheme we consider are known for the Cauchy problem.They were first introduced in the framework of the Boltzmann approach of hydrodynamic problems, see [9,16,17].A kinetic interpretation of flux splitting schemes is given in the paper by A. Harten, P.D. Lax, and B. van Leer [13].For general conservation laws, S. Jin and Z. Xin introduced a relaxation approximation and constructed related numerical schemes, which are equivalent to kinetic schemes with discrete velocities, for the Cauchy problem [12].A quite complete investigation on second order relaxation and discrete kinetic schemes for general systems of conservation laws in several space variables and with boundary conditions is developed in [1] and [3].
The interactions at junctions are solved by the use of a Linear Programming algorithm that compute the maximized fluxes for all the schemes.We define a numerical grid in R N × (0, T ) using the following notations: • ∆x is the space grid size; • ∆t is the time grid size; • (t n , x m ) = (n∆t, m∆x) for n ∈ N and m ∈ Z are the grid points.For a function v defined on the grid we write v n m = v(t n , x m ) for m, n varying on a subset of Z and N respectively.We also use the notation u n m for u(t n , x m ) when u is a continuous function on the (t, x) plane.

Fast Godunov scheme (FG).
A good numerical method to solve the equations along roads is represented by the Godunov scheme, which is based on exact solutions to Riemann problems, [10,11].The idea is the following: first the initial datum is approximated 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 and proceeds by induction.Let us now detail the scheme.We take a piecewise constant approximation of the initial datum: and the scheme defines v n m recursively starting from v 0 m .Solutions to Riemann problems from x m−1/2 are taken and then projected on a piecewise constant function by Since v is an exact solution of (1.1), we can use the Gauss-Green formula in (1.1) to compute v n+1 .Under the CFL condition ∆t sup the waves, generated by Riemann solutions, do not influence the solution in x = x m , for t ∈ (t n , t n+1 ).Hence the flux in x = x m for t ∈ (t n , t n+1 ) is given by , where W R x t ; v − , v + is the self-similar solution between v − and v + .As the flux is time invariant and continuous, we can put it out of the integral and, setting g G (u, v) = F (W R (0; u, v)) under the condition (3.3), the scheme can be written as: The expression of the numerical flux for Godunov method is in general given by and for the flux function under consideration: In the particular case of v = 1 and σ = 1 2 , the numerical flux determines the regions depicted in Fig. 6.
Now we display all the admissible configurations in the next table 7. Note that there are only twelve possible cases, schematized in Fig. 7.The numerical flux, for v = 1 and σ = 1 2 , reads: Finally, the scheme, under the CFL condition ∆t = ∆x, (3.6) Boundary conditions.Suppose to assign a condition at the incoming boundary x = 0: u(0, t) = ρ 1 (t) t > 0 and study equation only for x > 0. Now we are considering the initial-boundary value problem It is not easy to find a function u that satisfies (3.10) in a classical sense, because, in general, the boundary data cannot be assumed.One seeks a condition which is to be effective only in the inflow part of the boundary.Following [2] the rigorous way of assigning the boundary condition is: We practically proceed by inserting a ghost cell and defining where takes the place of v n −1 .An outgoing boundary can be treated analogously.Let x j < L = x N .Then the discretization reads: where takes the place of v n N +1 , that is a ghost cell value.
3.2.Fast Shock Fitting (FSF) method.The main framework of this method is determined by theoretical results of Section 2. In particular, we know that for solutions starting from an empty network, taking the flux as in (2.2), there exists one separating shock (possibly latent) on each road (see Lemma 2).Such property allows us to construct an exact procedure for the computation of the approximate solutions.This scheme is a finite volumes method, since we consider the values of densities on the cells C j = [x j−1 , x j ), with x j = x0 + j∆x and j = 1, . . ., N , instead of considering the densities on the nodes of the grid.FSF is composed by two main steps: 1) separating shock tracking; 2) shift by constant velocities.
Let us now describe the method starting from step 2, since is the simplest.Consider Riemann problems on each road, with constant states ρ − , ρ + .If ρ − , ρ + < σ we have a wave moving with speed v as showed by the case a) in Fig. 8.If ρ − , ρ + > σ, instead, the speed is −v, as represented in Fig. 8  in the following way.The points on the left respect to σ are shifted on the right of one index: and, analogously, the points on the right are shifted on the left of one index: with j the index of the space cells, n the index of the time nodes.For the numerical scheme, the shift reads as u n+1 j = u n j+1 for j > j * (forward) or u n+1 j = u n j−1 for j ≤ j * (backward).Now we describe step 1 of the algorithm.Under a suitable CFL condition which avoids the interaction of two neighbouring cells, namely ∆tv = ∆x, we capture the separating shock as follows.The separating shock is described by the 4-tuple (x, ρ − , ρ + , j), where x is the starting point, ρ − , ρ + are the densities, respectively, on the left and on the right respect to x and j indicates the cell containing x.We indicate by (x * , ρ * − , ρ * + , j * ) the shock parameters after a time step.For simplicity, we consider a scaling in the parametrization of each road I i .In particular, we set . The big shock is described by four parameters: (x, ρ − , ρ + , j).

3.2.1.
Evolution: separating shock tracking.We distinguish three cases: 1. Shock inside I i .2. Latent shock on the left.3. Latent shock on the right.
We set: intending that the index of the road i is fixed.

Case 1: Shock inside
We call x ∆/2 the position of the big shock at time ∆t 2 if no interaction occurs: We consider two cases: . Possible evolution of the big shock inside the cell.
In this case there is an interaction with the wave from j∆x and we want to compute the interaction point, namely (µ, y), where µ ∈ (0, ∆t/2), y ∈ ((j − 1)∆x, (j − 1/2)∆x).We get Then the position of the big shock at time ∆t, if no other interaction occurs, would be: We have two possibilities: Case IA.If j = 1 then we set x * = 0 and j * = 0, otherwise a second interaction verifies.After this interaction, the big shock position x * is: and x * is surely contained in the cell C j−1 due to its speed (< 1).The coordinates of the interaction point are (ν, z), with ν ∈ (0, ∆t), z ∈ ((j − 3/2)∆x, (j − 1)∆x).
After the second interaction, the big shock position x * is: The coordinates of the interaction point are (ν, z), with ν ∈ (0, ∆t), z ∈ ((j − 1/2)∆x, j∆x).We get: and we set: Case IIB.If j = N we set x * = N ∆x and j * = N + 1 , otherwise we have the second interaction.After this interaction the big shock position is: The coordinates of the interaction point are (ν, z), with ν ∈ (0, ∆t), z ∈ (j∆x, (j + 1)∆x).We obtain: and we let: Case 2: latent shock on the left boundary Lemma 3. If we have a latent shock at a i (on the left), then S i enters the first cell under the condition ρ 0 < τ (ρ 1 ).
In the case of a shock at the left boundary (as an example see Fig. 11), the big shock position at the interaction is: where y = λµ = ∆x − µ and µ = ∆x λ + 1 .
Now we have only the case x ∆ > 0, since the case x ∆ < 0 is already treated in IA.
In order to determine the position after the interaction we set: Proof.We have ρ N < σ.In order to exit from the last cell, the shock wave must have a negative speed, thus ρ N +1 > τ (ρ N ).
If we have a shock at the right boundary (as an example see Fig. 12), shock position at the interaction can be determined in the following way: where Now we have only the case x ∆ < 1, since the case x ∆ > 1 is already treated in IIB.We determine the position after interactions: Boundary conditions.If there is no latent shock on the left, we set: 1 < σ and we need to have a positive speed of the wave.If there is no latent shock on the right, we put: N > σ and we want the wave to have a negative speed.At the outgoing boundary it can also be applied a Neumann condition.Finally, the structure of the FSF method at each iteration in time can be schematized as in Fig. 13.Notice that for a network with junctions, the program necessarily involves the solution of Linear Programming (LP) problems for the maximization of the flux.Therefore, a LP solver must be taken into account, but the time required for the execution of any of such algorithms can be larger than the time taken by the rest of the scheme.This depends on which technique is chosen and on the fact that the solver can be internal or external to the simulation program.In order to better evaluate the efficiency of the proposed algorithms, we neglect the cost of the LP solver and focus on that of the numerical schemes.
Example 4.1.In the numerical tests we consider the following initial data: • constant null initial data ρ(0, x) = 0. (4.2) We take the boundary datum below: • reflected Brownian-motion boundary datum 4.1.Orders and errors.We introduce the formal numerical order γ of a numerical method for a single road in the following way: The L 1 -error on each road is where w M m (h) denotes the numerical solution obtained with the space step discretization equal to h, calculated in x m at the final time t M = T .The quantity w M m (h) denotes the numerical solution obtained with the space step discretization equal to h, computed in x m at the final time t M = T .In the following table T1 we used the exact solution to compute L 1 errors of the approximation schemes.Recalling that numerical order is computed by the approximate formula (4.5), whenever it is undefined we put the symbol "-" into the tables above.4.2.Accuracy: correct reconstruction of solutions.In the next Figures 14  and 15, we compare the solutions obtained by the Fast Godunov scheme (FG), by the 3-velocities Kinetic scheme of second order (K3V) and by the Fast Shock Fitting method (FSF) for different space (and time) steps.From the analysis of the previous tables it is easy to see that the Fast Godunov scheme (FG) performs better than the classical scheme.Moreover, we have that FSF performs better than FG and K3V as regards as the CPU time.Notice that in tables T6 and T7 the CPU time has an increment of about 0.50 seconds for FG and FSF respect to CPU times reported in tables T4 and T5, due to the fact that boundary data (4.4) is loaded at each time step.All the simulations have been performed by a personal computer, processor AMD Athlon XP 2400 Mhz, RAM 512 Mb.

5.
Conclusions.This paper introduces two new algorithms, Fast Shock Fitting (FSF) and Fast Godunov (FG), to compute solutions for traffic flow on a network for a piecewise constant flux function as in (2.2).Our tests show that FSF, whether applicable, represents an efficient and good approximation scheme for this problem.
FSF and 3-velocities Kinetic scheme of second order (K3V) substantially provide the same accuracy, thus, whether applicable, FSF is preferable regarding the speed of computation.FSF, indeed, results to be faster than FG and K3V and at the same time more accurate.The speed of computation represents an important aspect, also in consideration of the fact that road networks of cities can be composed by a considerably large number of roads.On the other hand, FG results faster than the classical Godunov scheme.

Figure 1 .Figure 2 .
Figure 1.Velocity as a function of the density.
shock and a) and b) hold.Now, if S i (t) ∈]a i , b i [ and a), b) hold, then a), b) hold true for every interaction with a wave from the right or on the left.If S i (t) = a i , then S i may leave the left endpoint only if a shock (ρ, σ), with ρ < σ, is generated from a i .Thus a) and b) still hold.If S i (t) = b i , then S i may leave the right endpoint only if a shock (σ, ρ), with ρ > σ, is generated from b i .Hence a) and b) still hold.
Wave with negative speed.