Runge-Kutta Discontinuous Galerkin Method for Traffic Flow Model on Networks
Suncica Canic 111Department of Mathematics, University of Houston, Houston, TX, 77204. E-mail: email@example.com. The research of the first author is partially supported by NSF under grants DMS-1263572, DMS-1318763, DMS-1311709, DMS-1262385, and DMS-1109189. , Benedetto Piccoli 222Department of Mathematics and Center for Computational and Integrative Biology, Rutgers University - Camden, Camden, NJ, 08102. E-mail: firstname.lastname@example.org. The research of the second author is partially supported by NSF under grant DMS-1107444. , Jing-Mei Qiu 333Department of Mathematics, University of Houston, Houston, TX, 77204. E-mail: email@example.com. The research of the third and the fourth author is partially supported by Air Force Office of Scientific Computing YIP grant FA9550-12-0318, NSF grant DMS-1217008 and University of Houston., Tan Ren 444School of Aerospace Engineering, Beijing Institute of Technology, Beijing, 100081. E-mail: firstname.lastname@example.org
Abstract. We propose a bound-preserving Runge-Kutta (RK) discontinuous Galerkin (DG) method as an efficient, effective and compact numerical approach for numerical simulation of traffic flow problems on networks, with arbitrary high order accuracy. Road networks are modeled by graphs, composed of a finite number of roads that meet at junctions. On each road, a scalar conservation law describes the dynamics, while coupling conditions are specified at junctions to define flow separation or convergence at the points where roads meet. We incorporate such coupling conditions in the RK DG framework, and apply an arbitrary high order bound preserving limiter to the RK DG method to preserve the physical bounds on the network solutions (car density). We showcase the proposed algorithm on several benchmark test cases from the literature, as well as several new challenging examples with rich solution structures. Modeling and simulation of Cauchy problems for traffic flows on networks is notorious for lack of uniqueness or (Lipschitz) continuous dependence. The discontinuous Galerkin method proposed here deals elegantly with these problems, and is perhaps the only realistic and efficient high-order method for network problems.
Keywords: Scalar conservation laws; Traffic flow; Hyperbolic network; Discontinuous Galerkin; Bound Preserving.
In this paper we deal with vehicular traffic models on networks. More precisely, we focus on the classical Lighthill-Whitham-Richards model (see [15, 16]), which consists of a single conservation laws for the car density. The model describes the evolution of traffic load on a single road, assuming that the average velocity depends only on the density via a closure relation. The resulting density-flow function is usually called fundamental diagram in engineering literature. Such model was adapted to networks in a number of different ways [13, 14, 6] depending on the rules used to describe the dynamics at junctions. The only conservation of cars is not sufficient to isolate a unique dynamics, thus additional rules, such as traffic distribution matrices, are to be prescribed. In particular, various authors proposed a set of additional rules which isolate a unique solution for every Riemann problem at a junction, i.e. a Cauchy problem with initial density constant on each road. In that case the map providing a unique solution to Riemann problems is called Riemann solver. A fairly general theory for such models on networks is now available, see [10, 11].
It is interesting to notice that lack of continuous dependence (or uniqueness) may indeed happen for Cauchy problems even if we do have unique solutions to Riemann problems. More precisely, the phenomenon of lack of Lipschitz continuous dependence is illustrated in Section 5.4 of the book . However, some Riemann solvers do provide Lipschitz continuous dependence for Cauchy problems (and thus also uniqueness). Examples can be found in ,  (Chapter 9) and . The specific solvers considered in this paper fall in this category, except for the case of two incoming and two outgoing roads (for which Lipschitz continuous dependence is false but uniqueness and continuous dependence are still open problems.)
Due to limitations of the single conservation law to describe dynamics in case of congestion, various models consisting of two equations (conservation of car mass and balance of “momentum”) have been proposed, see e.g. [1, 7]. Numerical methods for conservation laws on networks were developed mainly based on first order schemes, see [2, 8, 12]. First order schemes on networks have the same limitations as when they are applied to problems defined on a single real line: weak solutions are not well approximated, unless the spatial mesh is very fine to resolve solution structures. For this reason, we propose to use Discontinuous Galerkin methods with arbitrary high-order accuracy, which will be adapted in this paper to graph domains. Adaptation to graph problems requires supplementing the classical DG method with coupling conditions that hold at graph’s vertices. We propose the use of Runge-Kutta DG methods with total variation bounded limiters as a straight forward way of implementing the coupling conditions, while preserving the upper and lower bounds of DG solutions with a bound preserving limiter.
Since the late 80s, DG methods have been gaining great popularity as methods of choice for solving systems of hyperbolic conservation laws, with high order accuracy for smooth solutions and good shock capturing capabilities. We refer the reader to review papers and books [5, 3] for the history, development, and applications of the methods. The high order accuracy in time evolution is realized by applying the strong stability preserving (SPP) Runge-Kutta (RK) time discretization via the method-of-line approach. Compared with the existing high order finite volume and finite difference schemes, DG methods are more flexible with general meshes and local approximations, hence more suitable for h-p adaptivity. They are very compact in the sense that the update of the solution on one element only depends on direct neighboring elements, thus allowing for easy handling of various boundary conditions with high order accuracy and great parallel efficiency. Compared with the classical continuous finite element methods, DG methods are advantageous in capturing solutions with discontinuities or sharp gradients for convection dominant problems.
We propose to use the high order RK DG method with total variation bounded limiters as a general approach for simulating hyperbolic network problems. The compactness of the DG method enables a straightforward way of implementing coupling conditions at junctions. The bound preserving property of a first order monotone scheme for our traffic flow model on networks is theoretically proved, thanks to the rule of maximizing fluxes at junctions. Such property enables the application of bound preserving limiters, while maintaining classical high order accuracy of the RK DG method. Numerical results on benchmark problems from the literature, as well as on the ones with rich solution structures that we constructed in this paper, showcase the effectiveness of the proposed approach. We emphasize that, to our best knowledge, the DG method perhaps is the only realistic and efficient high order method for network problems. Existing high order finite difference and finite volume schemes would involve a wide and one-sided stencil in reconstructing solutions at junctions; such one-sided reconstruction stencil would lead to potential accuracy and stability issues. This paper is an initial step in applying the DG method to network problems; further development of the method to nonlinear hyperbolic systems with additional challenges in resolving junction conditions and numerical stability will be subject to future investigation.
The paper is organized as follows. Section 2 is on the background of traffic flow models on networks, with a general description of coupling conditions at junctions. Section 3 presents the proposed high order Runge-Kutta discontinuous Galerkin method for network problems with bound preserving properties. Section 4 demonstrates the performance of the proposed schemes on benchmark test problems from the literature and in challenging test cases with rich solution structures. Finally, a conclusion is given in Section 5.
2 Background on traffic flow models on networks
The nonlinear traffic model based on conservation of cars is a scalar hyperbolic conservation law in the form of
where is the density of cars, with being the maximum density of cars on the road; is the flux. The main assumption of this model
is that the average velocity is a function depending only on the density ,
thus giving rise to (2.1). The usual assumptions on
is that and that is strictly concave,
thus has a unique maximum point called the critical density.
Indeed, below the traffic is said to be in free flow
and the flux is an increasing function of the density. On the other side,
above the flux is a decreasing function of the density,
For future use, we define:
Let be the map such that for every and for every
A network is described by a topological graph, i.e. a couple where is a collection of intervals representing roads, and is a collection of vertices representing the junctions. For a fixed junction , a Riemann Problem (RP) is a Cauchy Problem with initial data which are constant on each road incident at the junction. The evolution on the whole network of the solution to (2.1) is determined once one assigns a Riemann Solver at each junction, i.e. a map assigning a solution to every Riemann Problem at the junction. More precisely, given initial conditions , where ranges over incoming roads and over outgoing ones, we will assign density values so that the solution on the incoming road is given by a single wave , and on the outgoing road by the single wave .
We consider the Riemann Solver based on the following rules:
There exists traffic distribution coefficients , representing the portion of traffic from incoming road going to outgoing road . The resulting traffic distribution matrix:
is row stochastic, i.e. for every it holds:
Respecting (A), drivers behave so as to maximize the flux through the junction. In other words the sum of the flux over incoming roads is maximized.
If a yielding rule, (C), is needed.
For example, consider the case of two incoming roads and and one outgoing road . Assume that not all cars can enter the road , and let be the amount that can do it. Then, cars come from the road and cars from the road .
Now we describe the solutions generated at junctions using rules
(A), (B) and (C).
Notice that solving a Riemann Problem at a junction is equivalent to solving Initial Boundary Value Problems (IBVP) on each road. Since solutions to IBVP may not attain the boundary values, due to the nonlinearity of the equation, one has to impose admissible values on each road, which generate only waves with negative speed on incoming roads and positive on outgoing ones. Indeed, if waves would enter the junction, then the solution to the IBVP may not attain the prescribed boundary value and, for instance, even violate conservation of cars, see also . In turn, this allows to state the problem in terms of fluxes, since densities can be reconstructed due to these restrictions. Moreover, we have some bounds on maximal flows on each road, more precisely we have:
Let be the initial densities of a RP at and and be the maximum fluxes that can be obtained on incoming roads and outgoing roads, respectively. Then:
In particular, densities can be recontructed by flows at the junction.
Consider first an incoming road and indicate by the trace at the junction for positive times. Then the IBVP on road is solved by a single wave: , which must have negative speed. If then either is or belongs to In the first case, there is no wave, while in the second case the wave is a shock wave with negative speed, see Figure 2.1 (left). Therefore the maximal flux is given by . Moreover, there exists a unique value of , which is compatible with a given value of the flux in the interval .
If, instead, then and the wave is a rarefaction or a shock wave with negative speed, see Figure 2.1 (right). In this case the maximal flux is given by and, again, there exists a unique value of , which is compatible with a given value of the flux in the interval .
For an outgoing road, the analysis is analogous, see Figure 2.2 . ∎
Proposition 2.2 allows to restate rules (A), (B) and (C) as a Linear Programming problem in terms of the incoming fluxes . Indeed, rule (A) allows to determine the outgoing fluxes in terms of the incoming ones. Then rule (B) provides a linear functional in the fluxes to be maximized. The constraints are given by the formulas (2.2) and (2.3). Rule (C) allows to choose a unique solution to the Linear Programming problem in case of more incoming than outgoing roads.
In the following sections we will explicitely solve the Riemann Problems in the following cases: junctions of type (two incoming roads and one outgoing road), junctions of type (one incoming road and two outgoing roads), and junctions of type (two incoming roads and two outgoing roads). We refer the reader to  for a complete description of the general case.
2.1 The case of = 2 incoming roads and = 1 outgoing road
Let us consider the junction with two incoming roads and and one outgoing road . Given initial data we construct a solution in the following way. To maximize the through traffic (rule (B)), we set:
where is defined as in (2.2) and as in (2.3). Notice that in this
case the matrix (or rule (A)) is simply given by the column
vector , thus it gives no additional restriction.
Consider now the space and the line:
defined according to rule (C). Let be the point of intersection of the line (2.4) with the line The final fluxes must belong to the region
There are two different cases:
does not belong to
The two cases are represented in Figure 2.3.
In the first case, we set while in the second case we set where is the point of closest to the line (2.4). Once we have determined and (and ), we can find in a unique way , :
Consider a junction with incoming roads and outgoing road. For every there exists a unique admissible weak solution at the junction , satisfying rules (A), (B) and (C), such that
Moreover, there exists a unique tuple such that
and for , the solution is given by the wave , while for the outgoing road the solution is given by the wave
2.2 The case of = 1 incoming road and = 2 outgoing roads
Let us now consider the junction with the incoming road and two outgoing roads and . The distribution matrix , of rule (A), takes the form
where and indicate the percentage of cars which, from road , goes to roads and , respectively. Thanks to rule (B), the solution to a RP is:
Once we have obtained and , it is possible to find in a unique way , , reasoning as in the proof of Theorem 2.3. Then we obtain the following:
Consider a junction with incoming road and outgoing roads. For every there exists a unique admissible weak solution at the junction , respecting rules (A) and (B), such that
Moreover, there exists a unique tuple such that
and for the incoming road the solution is given by the wave , while for the solution is given by the wave
2.3 The case of = 2 incoming roads and = 2 outgoing roads
Let us now consider the junction with two incoming roads and and two outgoing roads and . The distribution matrix , of rule (A), takes the form
where . We assume that , otherwise we may have more than one solutions to the Linear Programming problem, see  for details.
First notice that constraints from outgoing roads fluxes can be expressed as:
Define to be the point of intersection of the two lines:
To express the solution we need to distinguish some cases:
Case a). If and then the solution is given by:
Case b). If and then the solution is given by:
Case c). Assume and . If (thus ) then the constraint given by outgoing road is more stringent than that of outgoing road , thus the solution is given by:
Otherwise, i.e. if , then the solution is given by:
Case d). Assume and . If (thus ) then the constraint given by outgoing road is more stringent than that of outgoing road , thus the solution is given by:
Otherwise, i.e. if , then the solution is given by:
3 Numerical methods: RKDG
Below, we will first describe the RKDG method to discretize the 1D nonlinear traffic flow equations; then we will extend the algorithm to 1D network problems incorporating coupling conditions at the junctions. Finally, we will apply a high order limiter to preserve the upper and lower bounds of the high order solutions.
3.1 RKDG for 1D hyperbolic equations.
DG spatial discretization. Consider the following spatial discretization: let for be a partition of with being the cell center and being the cell size. The semi-discrete DG scheme for the equation (2.1) can be designed as finding numerical solutions in a finite dimensional space consisting of piecewise polynomials of degree ,
so that for any test functions ,
Here and below, the superscripts denotes the right/left limit of the function at a point. The function is a single valued function at the cell interface, which is defined via an approximate Riemann solver depending on the left and right limits of the DG solutions. One example is the global Lax-Friedrich flux, for which
For implementation, the approximate solution on mesh can be expressed as
where is the set of basis functions of . For example, the Legendre polynomials are a local orthogonal basis of with
where is the spatial operator, which denotes the R.H.S. of eq. (3.1), and is the numerical time step.
TVB limiters. When the solutions contain shocks, the TVB limiters proposed by Cockburn and Shu  will be used to eliminate spurious oscillations and enforce stability. Let be the cell averages of the numerical solution on cell , and let
The TVB limiter is used to adjust as follows,
where , and , and the modified minmod function is defined by
The limited and are then used to recover the new point values,
With the modified , as numerical solutions at the cell boundaries, as well as the cell average , we can construct a unique polynomial as the modified numerical solution.
Boundary conditions. Depending on the flow directions at the boundaries, one can prescribe either the inflow or outflow boundary conditions for the open boundaries. The RKDG method is well-known for its compactness, for which the inflow information or outflow information from extrapolation can be directly used in evaluating the flux specified in equation (3.2) at the boundary. Treatment of the coupling (boundary) conditions at network’s vertices is described next.
3.2 RKDG for hyperbolic networks.
The coupling conditions within a network, consisting of many incoming and outgoing roads with different junction points, are described in a general setting in Section 2. Below, we consider specific ways of constructing numerical coupling conditions, i.e., the numerical fluxes, at a single junction point of the following types: (a) one incoming and one outgoing road, (b) two incoming roads and one outgoing road, (c) one incoming and two outgoing roads, (d) two incoming and two outgoing roads. Similar coupling conditions can be derived for more complicated cases based on the rules (A), (B), (C) specified in Section 2.
Below, we assume that is chosen following equation (2.2) with being the right limit of the numerical solution on the cell of an incoming road at the junction point; is chosen following equation (2.3) with being the left limit of the numerical solution on the first cell of an outgoing road at the junction point.
One incoming road and one outgoing road. Let us consider the junction with one incoming road and one outgoing road . The coupling conditions are
The numerical fluxes of the junction point at incoming and outgoing roads are set to be
One incoming road and two outgoing roads. Let us consider the junction with one incoming road and two outgoing roads and . With the notation introduced for in Section 2.2, the coupling conditions are
The numerical fluxes at the junction point at the incoming and outgoing roads are
Two incoming roads and one outgoing road. Let us consider the junction with two incoming roads , and one outgoing road . With the notation introduced for in Section 2.1, the coupling conditions are the following. For the case of , we have
For the case of , the coupling conditions are
The numerical fluxes at the junction point between the incoming and outgoing roads are
Two incoming roads and two outgoing roads. Let us consider the junction with two incoming roads , and two outgoing roads . The coupling conditions follow those described in Section 2.3 for , with . The numerical fluxes at the junction point on incoming and outgoing roads are
One distinct property of the DG method compared with the other high order methods, e.g. finite volume or finite difference methods, for solving hyperbolic equations is the compactness of the scheme. Specifically, the high order fluxes depend only on the direct neighboring cells. Such property shows great advantage in prescribing boundary conditions at the junction points of hyperbolic networks.
3.3 Bound preserving numerical solutions
In the traffic flow model, it is known that . However, such a property does not hold for high order numerical solutions in general. To preserve the theoretical bounds on the RKDG solution, in this subsection we propose to apply the limiter proposed in . The application of this limiter is based on the fact that a first order monotone scheme with piecewise constant numerical solutions for network problems satisfies the following property.
Let denote numerical approximations to solutions at cell at time , and let be the spatial mesh size for a uniform mesh. (Similar results hold for nonuniform meshes.)
(Monotone scheme) A first order monotone scheme for the 1-D hyperbolic equation (2.1) can be written in the following form,
Here is a monotone flux, that is, it is a non-decreasing function with respect to the first argument, and non-increasing function with respect to the second argument. Similarly for .
It can be easily shown that in the monotone scheme, is a non-decreasing function with respect to , and , if a proper CFL condition is satisfied. Hence,
provided that the initial condition , leading to the bound preserving property of the numerical solution. The proposition below is a generalization of such a result to hyperbolic network problems. Without loss of generality, we consider the Godunov flux as our numerical flux, with
Consider a first order monotone scheme with Godnov flux as a numerical flux for the 1-D hyperbolic equation (2.1) holding on each road in a network, satisfying the coupling conditions at the junctions respecting the rules (A), (B), and ( C), with equations (2.2) and (2.3) specified in Section 2. Then, the numerical solution satisfies the bounds (3.17).
It is sufficient to prove the statement for the boundary elements adjacent to junctions. We consider the left-most element on a road (), which is an outgoing road at a junction. From equation (2.3), together with the rules (A), (B), and (C), we have . Hence
where the last inequality is due to the monotonicity of the scheme. In order to prove , we discuss two cases:
When , we have , hence
where the last inequality is due to the monotonicity of the scheme.
Similarly, when , we have , hence
Similar procedures can be done to prove the property for the right-most element on a road. ∎
In , a maximum principle preserving limiter is introduced for the RKDG scheme to preserve the maximum principle of the numerical solutions for hyperbolic PDEs, with the assumption that a first order monotone scheme satisfies the same property. The procedure of the maximum principle preserving limiter can be viewed as controlling the maximum and minimum of the numerical solution (polynomials on discretized cells) by a linear rescaling around cell averages. Such a procedure can be applied to control the bounds of the high order RKDG solutions for hyperbolic network problems. In particular, we would like to modify the numerical solution to , approximating a function on a cell , such that it satisfies
Accuracy: for smooth function , , on ;
Mass conservation property: ;
Bounds-preserving: on .
In order to achieve the above mentioned properties, one can apply the following limiter
where and are the maximum and the minimum of at Legendre Gauss-Lobatto quadrature points for the cell . It can be easily checked that with the application of such a limiter, the conservation and bound preserving properties of the numerical solution are satisfied. Furthermore, it was proved  that such a limiting process maintains the original order accuracy of the approximation.
Since the first order monotone scheme preserves the bounds for hyperbolic network problems, following similar procedures as in  one can show that the cell averages of the high order scheme are also well bounded, i.e. , , under the additional CFL constraint:
where ’s are the quadrature weights in the Legendre Gauss-Lobatto quadrature rule on a standard interval . Hence, the above limiter can be applied to the proposed RKDG scheme for hyperbolic networks. We also remark that if varies among different roads within a network, one can apply the similar limiter with the appropriate upper bounds on density.
4 Numerical examples
In this section, we reproduce simulation results from , using the high order RKDG method, discussed above, to compare the performance of the proposed high order scheme with the first order scheme used in . We also present several new examples with more complicated solution structures to showcase the advantages of high order schemes. In our numerical examples, for the third-order TVD Runge-Kutta method (3.4), we take CFL=1.0, 0.33, 0.20, 0.14 for , , and solution spaces corresponding to DG schemes with first to fourth spatial orders respectively. The time step for , and solution spaces, while for the solution space. The cell size is and the reference solutions are obtained by first order RKDG (finite volume method) with cell size in all examples, except the accuracy test.
4.1 Accuracy test
The first test is to solve the traffic flow equation (2.1) with the following flux function
with the initial condition
The computational domain is with periodic boundary condition. We compute the solutions up to time . Newton’s method is used to get the reference solution. We use a smaller CFL number, i.e. CFL=0.05 for the and cases to ensure that the spatial error dominates, so that the spatial order of accuracy can be observed for the scheme with the BP limiter. The results without and with BP limiter are shown in Tables 4.1 and 4.2. One can observe the -order convergence rate for solution spaces for the scheme with or without the BP limiter. The results without BP limiter show that the regular RKDG scheme produces numerical solutions that overshoot and undershoot the bounds of the exact solution. With the BP limiter, one can see that the scheme produces results that respect the bounds of the physical solutions.
We also tested the scheme with the following initial condition,
Figure 4.1 shows the RKDG solutions with solution space with and without BP limiter, superimposed over the exact solution.The numerical solution obtained using the scheme without the BP limiter displays some oscillations with the solution outside the physical bounds, while the results obtained using the scheme with the BP limiter fall inside the bounds, and approximate well the exact solution.
The simplest traffic flow model on networks is represented by the bottleneck problem. The conservation of cars is always expressed by (2.1), supplemented with initial and boundary conditions. The bottleneck problem models a road with different widths, hence different flux functions along different parts of the road. Denote the separation point between the two parts of the road by . We consider a road parametrized on with . We may consider the road as composed of two different roads. Let be the traffic density to the left of on (wider part) and be the traffic density to the right of on (narrower part). The wider part can be viewed as incoming road and the narrower part can be viewed as outgoing road. The flux function in the wider part is given by eq. (4.1), while the flux function in the narrower part is given by
The maximum of the fluxes is unique: