An efficient iterative link transmission model

ARTICLE IN PRESS

JID: TRB

[m3Gsc;February 1, 2016;21:17]

Transportation Research Part B 000 (2016) 1–21

Contents lists available at ScienceDirect

Transportation Research Part B journal homepage: www.elsevier.com/locate/trb

An efficient iterative link transmission model Willem Himpe a,∗, Ruben Corthout b, M.J. Chris Tampère a a b

Leuven Mobility Research Centre L-Mob, CIB, KU Leuven, Celestijnenlaan 300A, PO Box 2422, 3001 Heverlee, Belgium Transport & Mobility Leuven, Diestsesteenweg 57, 3010 Leuven, Belgium

a r t i c l e

i n f o

Article history: Received 30 January 2015 Revised 18 December 2015 Accepted 23 December 2015 Available online xxx Keywords: Dynamic network loading Macroscopic first order traffic simulation Link transmission model Warm start

a b s t r a c t In this paper a novel iterative algorithm is presented for the link transmission model, a fast macroscopic dynamic network loading scheme. The algorithm’s solutions are defined on a space–time discretized grid. Unlike previous numerical schemes there is no hard upper limit on the time step size for the algorithm to be numerically stable, leaving only the trade-off between accuracy and interpolation errors. This is a major benefit because mandatory small time steps in existing algorithm (required for numerical tractability) are undesirable in most strategic analyses. They lead to highly increased memory costs on larger network instances and unnecessary complex behaviour. In practice results are often aggregated for storage or analysis, which leads to the loss of computationally expensive detailed information and to the introduction of inconsistencies. The novel iterative scheme is consistent with the modelling assumptions independent of the numerical time step. A second contribution of the iterative procedure is the smart handling of repeated runs, which can be initialized (or warm started) by an earlier solution. For applications, repeatedly loading a network is often needed when evaluating traffic states under changing variables or adjusted parameter settings, or in optimization and equilibration procedures. In these cases the iterative algorithm is initialized with the solution of a previous run and iterations are performed to find a new consistent solution. Pseudo-code is provided for both a basic upwind iterative scheme and an extended algorithm that significantly accelerates convergence. The most important computational gains are achieved by ordering and reducing calculations to that part of the network which has changed (most). The properties of the algorithm are demonstrated on a theoretical network as well as on some real-world networks. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Dynamic network loading (DNL) models are designed to represent the propagation of traffic over time in road networks. They provide an accurate description of network performance including: saturation, time varying travel speed, the location of bottlenecks, queue formation/dissipation and blocking back due to upstream spillback propagation. A substantial amount of research has been carried out in the last decades. For an overview on the field of dynamic traffic flow models the reader is referred to the extensive classification of Wageningen-Kessels et al. (2014). This overview focuses on the description of traffic flow on links. Equally important for DNL is the role of intersections and their representation as nodes in the network, connecting boundaries of links and imposing capacity constraints themselves; for an overview of node models for DNL, ∗

Corresponding author. Tel.: +3216329614. E-mail address: [email protected], [email protected] (W. Himpe).

http://dx.doi.org/10.1016/j.trb.2015.12.013 0191-2615/© 2015 Elsevier Ltd. All rights reserved.

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB 2

ARTICLE IN PRESS

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

see Tampère et al. (2011). The focus here is on macroscopic traffic simulations which have the underlying assumption that traffic can be described as a continuum and compute the evolution of aggregated traffic characteristics such as density, flow or average speed. For many applications (strategic optimization, policy planning, dynamic traffic management, etc.) a simulation tool that provides an aggregated solution is sufficient and it is computationally more efficient than more detailed methods such as microscopic (car following) models. This research frames within the class of DNLs referred to as Link Transmission Models. It is a macroscopic DNL building on Newell’s simplified theory of kinetic waves (Newell, 1993a) and was first introduced by (Yperman et al., 2005). Link transmission models track traffic in a space–time discretized grid, for which they exploit the largest possible spatial discretization based on the longitudinal homogenous characteristics of a road or link – hence their name. However as discussed later in depth, temporal discretization in classical algorithms of LTM is often limited to small time intervals (∼ a few seconds). This is purely for computational reasons, as for analysis, such results are usually aggregated to larger time intervals (∼ multiple minutes). The first contribution of this paper is a method for handling sparse solution grids (with large time intervals), while still providing an accurate description of congestion patterns and the evolution of traffic states. A second contribution is an extension and optimization of this algorithm, herewith offering techniques for re-evaluating a scenario with altered input variables or parameters with a minimum of calculations. Such repeated calculations on the same network with only marginal changes (Corthout et al., 2014) are often required in practice for the computation of equilibrium routing behavior; the calibration of demand and supply; the optimization of network design, pricing instruments or control measures; the analysis of perturbations in robustness and reliability analysis. In the next section, a short overview is presented of the basic theoretical principles. On the link level, first order kinematic wave theory and the fundamental diagram are introduced to elaborate an analytical solution for a single homogenous link. Subsequently, links in a network are connected by incorporating node models and introducing a computational grid. In Section 3, a theoretical network is presented to illustrate the principles of the algorithm and motivate the core elements. First, the classical upwind scheme of the LTM is converted to an iterative approach such that large time intervals are made possible. Next, the initialisation of the iterative procedure is presented and its importance in the handling of repeated calculations, continued by the introduction of a node ordering strategy to speed up convergence. Computational efficiency is subsequently further increased by reducing redundant calculations to a minimum. This also provides significant increases in performance for repeated network loadings. The complete calculation scheme or algorithm is presented in Section 4. The algorithm is capable of simulating traffic on a general (large-scale) network for which the time varying demand and route choice is given as an exogenous input. Pseudo-code is provided for both a basic variant, limited to the iterative procedure, and an extensive variant, which includes all the speed-ups. In Section 5 the algorithm is tested on some real world networks with different morphologies and sizes. The result for different time intervals and speed-ups is investigated by comparing the congestion patterns and travel times as well as the computation times. The case studies are concluded by applying the algorithm to scenarios that require repeated calculations to illustrate the potential computation gains in practice. The paper is concluded with a short discussion on the contributions of the new scheme and its relevance. 2. Theoretical background 2.1. First order kinematic wave theory for a single homogenous link Many approaches to the macroscopic, analytical descriptions of traffic on a single homogenous road have been developed over the years. Its theory and impact on DNL and DTA literature are summarized next, to quickly converge to the elements that are important for the development of the link model for our network model. The fundaments of kinematic wave theory for traffic have been described by LWR-theory (Lighthill and Whitham, 1955; Richards, 1956) more than half a century ago. It is composed of a partial differential equation (PDE) that represents the conservation of vehicles, and a fundamental relation (Fig. 1) between instantaneous flow and local density called the fundamental diagram (FD). The LWR-PDE is described in function of the density of vehicles along a road. Classical first order finite difference approaches (Godunov, 1959) for numerically solving PDEs are widely used in practice, e.g. CTM of Daganzo (1994) or Lebacque (1996). These methods operate on a computation grid that discretizes space in cells and time in slices to provide an approximate solution to the PDE. The solutions have interesting physical properties like storage of vehicles (in terms of vehicle densities) along a link, the formation of shock waves and the relations between traffic states along characteristic speeds such as observed on real roads. However, the method is also computationally demanding requiring a link to be cut into smaller cells and updating times to be limited by the Courant Friedrich and Lewy (CFL) conditions (Courant et al., 1967). These CFL-conditions state that the time discretization should not be longer than the shortest possible travel time of traffic over a cell. With short cells and high speeds, this time step is in realistic networks typically in the order of a few seconds. In this paper, we provide a numerical scheme that circumvents these limitations. This work builds on the seminal work of Newell on simplified theory of kinematic waves (Newell, 1993a, 1993b, 1993c). He simplified the FD to a triangular shape that proves to be efficient for numerical schemes in realistic networks (Yperman et al., 2005), and more importantly, he used cumulative vehicle number (CVN) functions to represent traffic. This intermediate abstraction (Moskowitz, 1965) by CVN functions (referred to as N(x, t) at location x and time instant t) has also been of great importance for building exact solutions. It was formalized by Claudel and Bayen (2010b), Daganzo (2005b), Friesz et al. (2013) and Laval and Leclercq (2013) using theory of variations and Hamilton Jacobi (HJ) PDEs for general FD. Solutions Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB

ARTICLE IN PRESS W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

[m3Gsc;February 1, 2016;21:17] 3

Fig. 1. Triangular fundamental diagram calibrated on highway data. Diamond markers represent free flow traffic states and triangular markers represent congested traffic states.

to HJ-PDEs are described in function of CVN from which travel times and densities are easily derived (Szeto and Lo, 2005). As we explain in more detail in Section 2.3, the numerical evaluation of the CVNs involves the convex Legendre–Frechel transformation of the FD and the use of Lax-Hopf formulas to translate the initial values and boundary conditions into a traffic state at any location or time within the homogenous link (Mazaré et al., 2011). This Lax-Hopf formula can also be used to construct a grid that leads to exact and efficient solutions using dynamic programming techniques while assuming a simplified triangular FD (Daganzo, 2005a). Such grids become increasingly difficult to define in large networks (Raadsen et al., 2014) so a simpler rectangular grid is assumed, with a triangular FD defined for each link. The grid is evaluated by an iterative procedure such that only locally (at each grid point) the Lax-Hopf formula holds. This leads to faster algorithms, with substantially less grid points to be stored by the algorithm, while still computing a good approximate solution that is sufficient for practical applications. 2.2. Triangular fundamental diagram The existence of a fundamental relation between traffic states has been empirically proven and many approximations have been proposed in literature (Greenberg, 1959; Greenshields et al., 1935; Underwood, 1961). An example of observed aggregated1 traffic data (density and flow) on a highway lane and an approximate triangular shaped fundamental relation is shown in Fig. 1. Hypocritical states are represented by diamond markers and triangular markers represent hypercritical traffic states. The hypocritical and hypercritical branch of the triangular function can be found by simple regression techniques on both sets of points. The slope of both linear approximations represent the characteristic wave speeds, in this case the free flow speed (γ ) in uncongested traffic and the spillback speed of perturbations in congestion (−ω). Also the maximum storage capacity or jam density of a link (j) and the maximum throughput (c) at the critical density (kc ) are defined on the graph. The triangular shape is often used because it leads to efficient numerical solution schemes (CTM, LTM). The Lax-Hopf equations (see discussion in next section) are reduced significantly in this case because only two characteristic kinematic wave speeds (one forward and one backward) govern the prevailing traffic states (Claudel and Bayen, 2010a; Daganzo, 2005a; Han et al., 2012; Jin, 2015). In principle the iterative procedure that is proposed in the next sections can be adopted to handle any concave FD. It is chosen to omit it here because it does not contribute to our main argumentation concerning avoiding CFL-conditions and efficiently computing repeated runs. 2.3. The Newell–Luke principle and the corresponding Lax-Hopf formula The theory behind the LTM link model can be understood by solely exploiting the simplified kinematic wave theory of Newell (1993a). This is also how it was originally introduced in Yperman et al. (2005) and Yperman and Tampère (2006). In this section however, it is related additionally to the more recent stream of research that frames LTM within the more general theory of Hamilton–Jacobi partial differential equations and its Lax-Hopf solutions, and hence gives LTM a more general theoretical foundation. This is done by briefly summarizing and physically interpreting (whenever possible) the various theoretical concepts and refer to more specialized literature for the rigorous mathematical background and proofs. Let us consider the Hamiltion–Jacobi partial differential equation (HJ-PDE) solutions on a homogenous link with appropriate initial values and boundary conditions. For further details on the mathematical properties of this PDE and its solutions for the traffic flow problem, the reader is referred to Daganzo (2005b), Han et al. (2012) and Mazaré et al. (2011). To illustrate the general approach to formulating the HJ-PDE, let us consider a simple road of length l with two detectors; one at the upstream end and one at the downstream end. Traffic on this road is governed by a FD (referred to as Q(k)) with the 1 The data is collected by double loop detectors on the E313 highway near the city of Antwerp. Observations of individual cars were aggregated in intervals of 10 minutes using the formalisms described in the Handbook of Treiber and Kesting (2012).

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB 4

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

Fig. 2. Trajectories of individual cars crossing 2 detectors along a highway and the relation between the cumulative vehicle number functions at these detectors in free flow (a) and congested regimes (b).

same functional form as the previously defined triangular FD(γ , ω, c, j) . The following HJ-PDE in function of N(x, t) can be constructed by stating that all traffic states have to be on the FD while by definition the partial derivative of the CVN N(x, t) in time (t) is equal to the instantaneous flow (q) and the upstream (note the negative sign) spatial derivative is equal to the localized density (k):

  ∂ N (x, t ) ∂ N (x, t ) =0 −Q − ∂t ∂x

(1)

In HJ-PDE theory, it is customary to formulate the analytical continuous solution to Eq. (1) through introduction of an auxiliary convex transformation of the FD called Legendre–Frechel transformation, which represents the flux of vehicles that are encountered by a moving observer with speed ϑ:

R(ϑ ) = sup (Q (k ) − ϑ · k ), k∈[0, j]

∀ϑ ∈ [−ω, γ ]

(2)

The boundary values of this transformation are easily found for the triangular FD(γ , ω, c, j). If the slope ϑ is equal to γ it results in 0 and if the slope is equal to −ω the result is ω · j. Only these extreme values are important, as all traffic states will only be governed by these two kinematic waves in a HJ-PDE with a simplified triangular FD (Jin, 2015). As Han et al. (2012) show, the traffic state solution to Eq. (1) on the homogeneous link is expressed in function of the CVN boundary values C(x, t) at the upstream and downstream end of the link along these two kinematic waves as expressed by the so-called Lax-Hopf formula solving the HJ-PDE:

N (x, t ) = inf {C (t − τ , x − τ · ϑ ) + τ · R(ϑ )},

τ ≥ 0 and ϑ ∈ [−ω, γ ]

(3)

Next, it is shown how the model formulation and solution algorithm of LTM with triangular FD are based on a pragmatic solution for Eq. (3). First, let us start with a physical interpretation of the mathematical solution that is to be constructed. Consider the link to be in free flow conditions. Any vehicle that crosses the upstream detector will cross the downstream detector some time later. Looking at the most left graph in Fig. 2a, one can observe that the trajectories of individual cars, represented by continuous lines, in space and time first traverse the upstream detector to continue their way until they cross the downstream detector. Because the link operates in free flow conditions, traffic is governed by characteristic waves with slope γ . An observer moving along such a wave does not cross any vehicle trajectories while traversing the link. As explained before, this is also expressed by R(γ ) ≡ 0. As a result, the CVN the observer recorded on the upstream detector at his departure must be equal to the CVN at the downstream detector at arrival and anywhere along the wave. The relation between the cumulative counts of vehicles U at the upstream detectors and V at the downstream detector is illustrated by the right graph in Fig. 2a. The resulting downstream CVN function is a simple translation of the upstream CVN function formalized in Eq. (4).

V (t ) = U (t − l/γ )

(4)

Next, if the link is congested, it is the other way around. The upstream CVN function is now governed by the downstream cumulative exiting rate because traffic is governed by upstream moving kinematic waves with slope ω. An observer moving along a wave will now cross multiple trajectories. The observed flow rate is given by R(ω) ≡ ω · j. The total vehicle amount observed between the two detectors integrates to j · l with l the length of the link between the two detectors because the observer takes l/ω time to traverse the link. This relation is illustrated by the two graphs in Fig. 2b. The CVN function of the upstream detector is a translation of the downstream CVN function, translated in time by the spillback travel time, and in cumulative number by a fixed additional amount of vehicles as given in Eq. (5).

U (t ) = V (t − l/ω ) + j · l

(5)

Finally, let us consider the CVN function N(x, t) at an arbitrary location x along the link. The solution can be derived only by translating and combining the boundary conditions (upstream and downstream CVN). According to the Newell-Luke principle, it is given by the lower envelope of what can be sent by the upstream side (Eq. (4)) and what can be received by Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

[m3Gsc;February 1, 2016;21:17] 5

downstream side (Eq. (5)). It is based on the assumption that each driver individually tries to move as far downstream as possible. The result is equal to the Lax-Hopf formula (Eq. (3)) that is reduced to the two terms of Eq. (6) in this simplified setup and it summarizes the core principle of our link model (Han et al., 2012; Jin, 2015).

 

N (x, t ) = inf U t −

x

γ





,V t−

l−x

ω



 + j · (l − x )

(6)

The LTM link model is obtained by applying Eq. (6) to the most up- and down-stream ends of a link, i.e. for x = 0 and x = l, respectively. This is elaborated further into a conceptual and practical algorithm in Sections 3 and 4, respectively. 2.4. Transferring flow with a first order node model In DNL schemes, the node model provides the connection between the homogeneous pieces of road. An accurate description of bottlenecks and intersections is extremely important in a network because they form the main source of delays and congestion. The node model adopted here is a point-like macroscopic description of the corresponding spatial structure. With the exception of some earlier works (Adamo et al., 1999; Buisson et al., 1996; Daganzo, 1995; Holden and Risebro, 1995), these node models have only recently regained increasing attention in literature (Bliemer et al., 2014; Flötteröd and Rohde, 2011; Guido Gentile and Tiddi, 2009; Gibb, 2011; Jin and Zhang, 2003; Lebacque, 2005; Ni et al., 2006; Smits et al., 2015; Tampère et al., 2011). The node model used here is an implementation of the generalized first order macroscopic node model that is presented in Tampère et al. (2011). It fulfils the seven requirements that need to be satisfied by any node model in order to propagate flow consistently. Its solutions have properties like uniqueness and stability (Corthout et al., 2012), an important property for practical applications. In principle the node model here can simply be seen as an interface. It returns the amount of vehicles that can be transferred from the incoming links into the outgoing links during a specific time interval. In case of a simple general intersection, the only arguments that need to be provided are (Tampère et al., 2011): the maximum amount of vehicles that arrive from each incoming link during the time interval (sending flow); the maximum amount of vehicles that are allowed to flow into each outgoing link during the time interval (receiving flow); split rates or turning fractions to distribute flow over outgoing links according to some route choice model; and a redistribution rule for handling competing turn movements towards restricted outgoing links. The focus of this paper is on the introduction of a novel solution scheme of a network DNL, which is in itself not limited to any specific node model. Although in principle any consistent node model can be used, some issues related to intersection delay and capacity are worth mentioning. These phenomena are important for the resulting flows and travel times in the DNL and are governed by signal timings (red-green alteration, offsets) in signalized, and conflict handling (priorities) in priority junctions and roundabouts. The modeller has two options: to model these phenomena as (rapidly) time-varying boundary constraints or a continuum approach. In order to take advantage of the larger time steps enabled by our LTM algorithm, a continuum approach is recommendable both for capacity (Han and Gayah, 2015) and for delays (Yperman et al., 2007). However, on both issues there are still challenges left that warrant further research, like internal supply constraints (Corthout et al., 2012; Flötteröd and Rohde, 2011) and how to consistently combine transient delay formulas with LTM (Blumberg-Nitzani and Bar-Gera, 2014; Yperman et al., 2007). It is chosen to omit further details on the specific node model used here (Tampère et al., 2011) and focus on the connection between node model and link model. As mentioned in the previous sections, the link model operates in CVN space. The node model on the other hand operates on flow rates (i.e. sending and receiving flow) that are assumed to be constant within a certain time interval. To connect both models in a consistent manner, CVN function on the link ends are assumed to be piecewise linear. Although the piecewise linear shape of the CVN functions also follows from the triangular FD, we only allow CVN to change rate at predefined time steps. These time steps are not necessary coinciding with the vertices found by the analytical formulation. This conjecture leads to the statement that the link model in our algorithm will only be consistent at grid points of the space–time discretization and it is not forced to be consistent anywhere in between. 2.5. Computational grid The solution is defined on a grid of discrete points in space and time such that it can be solved by a computer program with numerical precision. The algorithm will search for matching CVN values at these grid points. A wide range of possible grids has been proposed so far (Daganzo, 2005a). Within the class of Link Transmission Models (LTM), grid points are always positioned at one of the link ends, based on the fact that links are assumed homogenous and internal states can be fully explained by the boundary conditions (Eq. (6)). Because grid points are situated at the link ends, they are grouped together at the node level such that all connecting link ends are evaluated at the same time. For this reason, the evaluation of these grouped grid points is called a node update. The original LTM (Yperman et al., 2006) was constructed on a predefined updating scheme of the nodes that minimized updates and avoided CFL violations. It required for each node to be updated at fixed time tics specific to the characteristics of the connecting links. Other algorithms such as GLTM (Gentile, 2010) require global time tic updates. This means that all grid points that make up a time slice are updated before moving to the next time slice. Again because of CFL-conditions, the time interval between two time slices is here limited to small update steps. Recently Raadsen et al. (2014) presented an event based algorithm that creates the solution grid at run time. It is completely dependent on the wave tracking of changes in flow rates, making it an efficient solver especially if departure rates and Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB 6

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

Fig. 3. Rectangular solution grid with a short time interval (a) and a long time interval (b) for a simple 2 link/3 node network.

turning fractions are stable for long time periods. Although this makes it also cumbersome for repeated evaluations, e.g. when computing the influence of small flow changes, which then will also alter the solution grid. A core contribution of our algorithm is in fact the efficient handling of such repeated evaluations and it is for this reason that we choose a predefined grid. Finally, the presented computation grid is not constrained in its time discretisation (Bellei et al., 2005), except for what is desired as an output accuracy level. For DTA procedures such output aggregation intervals are typically larger than CFL conditions (e.g. multiple minutes). 3. Methodology and motivation 3.1. Theoretical network In the next 5 sections the core elements of the algorithm are introduced and illustrated on a small two-link/three-node network (Fig. 3). This elementary network allows us to use a simplified notation and to focus on the main principles behind our numerical scheme. In this way, the reader can understand the concept of the algorithm before it is generalized. Further details, combined with a more precise notation for general networks are presented in Section 4. The numerical scheme operates on a rectangular grid; Fig. 3 illustrates the grid for the simple network. The grid points (◦) are positioned at the link ends and are evenly spread in time. For both links (a and b), the FD is identical and CVN values are stored at the upstream and downstream end. At node n2 CVN curves are identical because for this simple case this node has one incoming and one outgoing link for which all flow entering the node has to exit it as well. The time slices τ i for which consistent solutions are constructed are separated by a fixed time interval τ . In our algorithm, unlike the original LTM, it is a fixed parameter independent of time or node but no longer limited by CFL-conditions. A full solution will have CVN values for each grid point. The CVN solutions at a node are constrained by the upstream or downstream CVN functions of the connecting links along the kinematic waves (illustrated by the dashed and point-dashed line in Fig. 3). To find a constraining CVN value (•), one can interpolate linearly between two time grid points, following our earlier assumption that at the nodes, CVN functions are piecewise linear in time. 3.2. Upwind iterative procedure The algorithm is an upwind scheme because of the time dependency. This means that grid points in Fig. 3 are evaluated column by column (or more correctly time slice by time slice), moving forward in time. Within each column, CVN values X(τ i ) are obtained by repeatedly visiting each grid point, in a Gauss–Seidel type of iterative scheme (Weisstein, 2015), while all CVN values of earlier time slices (U(τ < τ i ) or V(τ < τ i )) are fixed. In general, three sequential steps are taken to update a node and to find a CVN value at its grid point. In the following list, these three steps are presented for node n2 in time slice τ i . (a) First the sending flow (SF) and receiving flow (RF) are determined based on the upstream/downstream CVN of connected links looking along the kinematic waves (with slope γ and ω). The SFa and RFb represent the maximum number of vehicles that can arrive at or depart from node n2 within the considered time interval for links a and b. For this reason, the total amount of vehicles that already crossed the node in previous time steps, has to be subtracted from the constraining CVN. Negative flow rates are not allowed and total flow also has to be limited to the capacity c of the sending or receiving link (Eqs. (7) and (8)). Both links are assumed identical in our example to reduce the indexing; more generally, the SFa and RFb need to be calculated dependent on the characteristics of each connecting link. If a kinematic wave intersects with the upstream or downstream end in between two time grid points, a linear Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

[m3Gsc;February 1, 2016;21:17] 7

interpolation is made and if it takes more than τ time for waves to traverse the link, the CFL condition holds (illustrated by Fig. 3a). Because all related CVN values have then already been fixed in earlier time slices, a single update of the node is enough (similar to a dynamic programming type of algorithm). Specifically for this example where both waves arrive in between points τi−2 and τi−1 , we get:

SFa = min (c · τ , max (0, (1 − ϕ ) · Ua (τi−2 ) + ϕ · Ua (τi−1 ) − Va (τi−1 )))

(7)

RFb = min (c · τ , max (0, (1 − ψ ) · Vb (τi−2 ) + ψ · Vb (τi−1 ) + j · l − Ub (τi−1 )))

(8)

ϕ= with :

ψ=

2 · τ − l/γ

τ

2 · τ − l/ω

τ

If it takes less than τ time for waves to traverse the link (illustrated by Fig. 3b), then the sending and receiving flows are also dependent on the solution of neighbouring grid points within the evaluated time slice itself. In this case, they are implicitly linked to each other and the result of one grid point is dependent on the solution at the other grid points and vice versa. A solution is then found by iteratively evaluating these grid points while filling in the solution of the last iteration, in iterate m + 1 the values of Xm . Because only known information is used from a previous iteration, the order of updating of the nodes within an iteration does not matter. One could however also exploit the latest values of grid points that were already updated within the current iteration (X m+1 ) similar to the Gauss–Seidel procedure for solving a system of linear equations. In this case, convergence of the iterations will depend on the order of updating nodes within each iteration. As an example of the Gauss–Seidel procedure, consider a downstream sweep from origin to destination over the nodes in the two link network. In this case the SFa of node n2 is based on values evaluated in the current iteration (Eq. (9)) while the RFb is based on values of the previous iteration (Eq. (10)). Hence, for a fully free flowing scenario, convergence would be achieved in just one iteration. Likewise, for a fully congested scenario, fast convergence would require choosing an upstream sweep over the nodes. Obviously in more generic cases, such optimal sweeping orders are not that trivial. In Section 3.4 an efficient sweeping order is developed that is able to cope with the more complex interactions that emerge in complex networks.











SFa = min c · τ , max 0, (1 − ϕ ) · Ua (τi−1 ) + ϕ · Xam+1 − Va (τi−1 )

RFb = min c · τ , max 0, (1 − ψ ) · Vb (τi−1 ) + ψ · Xbm + j · l − Ub (τi−1 )



(9) (10)

τ − l/γ τ τ − l/ω ψ= τ

ϕ=

with :

(b) Next the SF and RF are used to evaluate the node model. In case of a single incoming link and a single outgoing link, like for node n2 , the node model is reduced to Eq. (11). Recall that this is similar to the Newell–Luke principle and it states that the flow over the node is maximized. It is a transformation of the general rule that each driver individually tries to move as far downstream as possible (Ansorge, 1990). Here the node model returns a single value but in general the node model will calculate turning flow rates (TF) for each turn movement that crosses the node. In addition, the node may also impose constraints itself (e.g. because of conflicting prioritized turn movements, or traffic signal control).

T Fab = min(SFa , RFb )

(11)

(c) In the final step, the CVN values are updated. This is a trivial operation in case of the simple node model (Eq. (12)). In practical networks, multiple incoming and outgoing links will connect to a single node and the CVN values for the entire group of grid points of that node will need to be updated. If for every connected link the CFL conditions holds, all calculated CVN values are stable because they are based on CVN values that have been fixed in earlier time slices. In this case X 1 is immediately fixed to Va (τi ) = Ua (τi ) and is not considered for updating in the next iterations.

X m+1 = Va (τi−1 ) + T Fab = Ua (τi−1 ) + T Fab

(12)

After updating all grid points, a new iteration is started as long as the solution has significantly changed with respect to the previous iteration. If for all grid points of time slice τ i , Eq. (13) is fulfilled, the solution is fixed and the algorithm moves to the next time slice. The precision level prDNL is defined beforehand, typically in the range of numerical precision of the computations.

 m  X − X m+1  < prDNL

(13)

The use of iterations to overcome CFL conditions in a DNL is first introduced by (Bellei et al., 2005) which is implemented as DUE in VISUM. In their algorithm iterations are made over the entire simulation period, contrary to our approach which considers iterations for each time slice. Although no proof of convergence for our algorithm is made here, the case studies Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB 8

ARTICLE IN PRESS

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

in Section 5 and the wide application in practice of VISUM provide a clear case for further developing simulations with large time intervals. The convergence of our algorithm is guaranteed by considering that in Eqs. (9) and (10) the factor multiplied with X is always smaller than 1, because it is the interpolation fraction based on a wave crossing in the last time interval. This means that any change in CVN of a node in iteration m leads to a smaller change on the neighbouring nodes in iteration m + 1. A rigorous analysis of the convergence properties is left for future work. The interested reader is referred to the fixed-point formulations of (Gentile, 2015) and more mathematical work of Zhang et al. (2006). 3.3. Initialization and repeated calculations The iterative upwind scheme avoids small time steps, which are typically present because of CFL conditions in noniterative schemes. However, the iteration procedure needs to be initialized. If no information is available about a new time slice, the solution can be initialized for instance by setting Xa0 = Ua (τi−1 ), Xb0 = Vb (τi−1 ) and X 0 = Va (τi−1 ) = Ua (τi−1 ). This would mean that the iterative scheme starts by assuming zero flow over the entire network during the time slice considered and will update the solution until the flow is maximized over all grid points (subject to all prevailing propagation constraints of course). Alternative initializations for the iterations of a time slice are conceivable, for instance instead of zero flow assume stationary conditions, which means that initial flows in the new time slice are equal to those in the previous time slice, and let the iterations correct this initial guess. All initialization options that use no other information than results from the DNL scenario that is currently simulated are called ‘cold start’. All cold start scenarios in this paper use zero flow initialization. Most importantly, the iterative scheme can also be initialized by CVN values of a previously calculated DNL. This is referred to as ‘warm starting’ the DNL and provides interesting efficiency gains when repeated calculations of the DNL are required. The importance of efficient algorithms for repeated DNL runs was brought to attention by Corthout et al. (2014). Many applications require the DNL to be evaluated multiple times on the same network with different settings or scenarios (Cascetta, 2001; de Palma and Lindsey, 2011; Frederix et al., 2014). Often, supply (topology, link and node characteristics, control settings of traffic signals, etc.) and demand (route choice, OD matrix) parameters are very similar in the different simulations; and therefore also the resulting CVN functions show a strong resemblance with previously simulated runs. The initialization of iterations with CVN values of a resembling run leads to a lower number of iterations and shorter computation times as illustrated by the case studies of Section 5. If the computational grid is exactly the same, the use of resulting CVNs of a previous network loading as a starting point is straightforward. Simply replacing the initialization of X0 in the iterative procedure is sufficient for warm starting the procedure. The next two sections introduce further efficiency gains to the basic updating scheme of the upwind iterative scheme. For the formal notation of the basic DNL procedure on general networks, the reader is referred to Section 4 and the pseudocode. The pseudo-code in bold represents the core working of the algorithm, which was introduced up to this point. Further algorithmic refinements introduced in the next sections are in plain text and can be omitted to construct a simpler, more comprehensible version of the algorithm, which would produce results of the same accuracy but less efficiently (requiring more iterations and more computation time per iteration). 3.4. Fast sweeping The number of iterations of the algorithm can be reduced by ordering the grid points selected for updating (Zhao, 2005). As illustrated before, if all grid points can be ordered along a characteristic wave, a single iteration is enough for convergence. In practical networks, such simple ordering is not evident because of the interaction of different characteristic waves at each node. Next, a fast sweeping algorithm is proposed based on the magnitude of the CVN differences at grid points between two iteration steps. The motivation behind this heuristic approach originated from the observation that large changes tend to sprawl over a larger part of the network and overrule the smaller changes. The algorithm focuses on the large changes first to speed up the convergence. Although our numerical results do not indicate such ordering as critical for convergence, its benefits are certainly significant and more important when considering repeated calculations. Like most LTM solution schemes, the algorithm proposed here consists of a series of node updates. Each node update constructs a consistent solution for the grid points of all connected links. For this reason, the changes of all grid points influencing a node update are grouped together. It is the combined potential for change based on the neighbouring nodes that determines the ordering of the node updates. Eq. (14) provides a quantitative measure for the combined change for node n2 in the simple two link network. First, the meaning of the second and third terms is explained. The second term describes the potential influence of upstream grid point changes that move downstream. The third term captures the potential influence of changes moving upstream because of congestion at a downstream grid point. The weight functions δf f and δ con are set equal to the interpolation fractions of Eqs. (9) and (10) respectively ϕ and ψ . Consequently, they must be set to zero if the change of a grid point is unable to reach the node within that time interval (i.e. when CFL conditions are met on the link). The first term captures the change of the grid point at the node itself compared to the previous iteration. It should be included if an argument of the node model is calculated based on the prevailing flow conditions (like a turn fraction, if it is part of the simulation, or feedback control settings like traffic responsive signals, e.g. ramp metering). δn2 is set to one in this case and to zero otherwise. For nodes with multiple in- and outgoing links, all three terms should be Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB

ARTICLE IN PRESS W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

[m3Gsc;February 1, 2016;21:17] 9

summed over all influencing links.

      2 = δn2 X m+1 − X m  + δf f Xam+1 − Xam  + δcon Xbm+1 − Xbm 

(14)

Before each new iteration, the nodes are arranged in descending order using metric 14. We don’t recommend any specific sorting strategy for finding such descending order. Nevertheless, it is important to note that the sorted list is in some cases very similar between sequential iterations, which is a known worst-case scenario for certain types of sorting algorithms (like quick sort). The stopping criterion for the iterative procedure is now redefined in terms of the potential for change. If for all nodes the potential for change drops below the precision level, the iterations have converged and the algorithm moves to the next time slice. The computational properties of the fast sweeping method are further explored in Section 5. 3.5. Reducing redundant calculations The computational effort of each iteration can be further reduced by avoiding redundant calculation. Up to this point, it was assumed that every node in the network is updated within each iteration. The procedure for evaluating the CVN updates at a node is a costly operation. However, if no significant potential for change on the grouped grid points is recorded (according to Eq. (14)), updating the node will not (significantly) change the result. A parsimonious handling of these updates is proposed by intelligently shaping the set of nodes to be updated. The set of nodes that is checked in each iteration is called the active set. After each iteration, the active set is updated by selecting only those nodes that have sufficient potential to change. Nodes that are sufficiently similar and have stable neighbours are not considered in the set and are de-activated. During the iteration scheme, the active set is repeatedly changed. A node might be part of the active set multiple times alternated with periods of inactivity. The potential to change the grid points of a node is measured by the same metric (Eq. (14)) used for ordering the node updates. For example, node n2 is part of the active set if 2 is larger than a pre-set precision level. If the active set is empty, the iterative procedure has converged. Because of the Gauss–Seidel type of evaluating the variables, the calculation of the metric n is integrated into the iteration scheme. To illustrate the need for this integration, consider that the potential for updating a node is based on the change of the grid points at its neighbouring nodes. After a node update, there is no potential for updating that node again (if the first term in Eq. (14) is zero) until one of the considered grid points of the neighbouring nodes changes. For some nodes, such a change will already be taken into account because the node update is based on information of the current iteration. To handle this complicated interaction, Eq. (14) is calculated during the evaluating of the nodes in the active set. Every time a node is updated, its own potential for updating is set to zero. Since its grouped grid points have changed, it must inform its dependent nodes (also inactive ones) of this change. It does so using the same three terms of Eq. (14). The first term is added to the now zero potential of the node itself. In the simple two link network, for node n2 this leads to Eq. (15). The second term is added to the downstream node’s potential of the outgoing links (Eq. (16)). The third term is added to the upstream node’s potential of the incoming links (Eq. (17)). As a result at the end of an iteration, every node has a potential change that is based only on the set of nodes that have been updated after its own update.

  2 = δn2 X m+1 − X m 

(15)

  3 = 3 + δf f X m+1 − X m 

(16)

  1 = 1 + δcon X m+1 − X m 

(17)

The active set needs to be initialized at the start of the iterative procedure (i.e. when a new time slice is considered). First, the case of cold starting the algorithm is elaborated on. Since all vehicles start at origin nodes, those nodes are part of the active set in the first iteration of each time slice such that they can start the knock-on effect of subsequently visited nodes. In later iterations of the same time slice, origin nodes should not be considered because the CVN values are fixed to the boundary value. Vehicles that are already on the network also need to be propagated towards their destination. For this reason, those nodes that have been changed during the previous time interval are always selected in the initial active set of a new time slice. In case of repeated DNL simulations with a warm start, the CVN values at the grid points are initialized with a previous result. Under the assumption that this result is consistent with the solution of the presented DNL algorithm for a combination of nodes and time slices, only those nodes that are inconsistent are added to the initial set of the iterative procedure for the concerning time slices. Similar to the cold starting procedure, also all nodes that have changed during the previous time slices are selected. If the change to the simulations input has only a local effect in space or time, the set of active nodes identified by Eqs. (15)–(17) will remain small and calculations will be performed only there where they are necessary. This leads to a parsimonious updating scheme that allows for a very efficient calculation of marginal changes. It is also completely embedded into the general structure of the algorithm making it less complex to implement than more ad-hoc approaches that have been presented in literature (Corthoutet al., 2009, 2014). Further computational analysis of the activation procedure is presented in Section 4. The general formulation of the updating process is introduced with pseudo-code in the next section. All the added features such as the selection, updating the potential change, and sorting of the nodes for fast sweeping over the set of active nodes are written in plain text compared to the core algorithm in bold. Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB 10

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

4. Algorithm description 4.1. Framework and notation In this section, a mathematical representation of the complete algorithm for general networks is constructed. The notation is mostly adopted from Gentile (2010) and Raadsenet al. (2014). An overview of the symbols used in this section can be found in Appendix. The traffic network is defined as a directed graph G(N, A), where N is the set of nodes and A is the set of links. All nodes that are sources of flow are part of the origin nodes subset O⊆ N and all nodes that are sinks are part of the destinations nodes subset D⊆ N. Vehicles travel between origin–destination (OD) pairs (o, d) with o ∈ O and d ∈ D. An arc a ∈ A between two nodes is directed with an initial node or tail node tla ∈ N and a final node or head node hda ∈ N. Each directed arc or link is characterized by a length la (km) and a triangular FD,2 which is defined by a link specific free flow speed γ a (km/h), spillback speed ωa (km/h), capacity ca (veh/h), and jam density ja (veh/km). Most nodes are connected to multiple links. The set of incoming links (backward star) of a node n ∈ N is denoted as BSn = {a ∈ A : hda = n} and the set of outgoing links (forward star) of a node is denoted as F Sn = {b ∈ A : tlb = n}. Note that the network is constructed such that origin nodes have no incoming links BSn = ∅ if n ∈ O and destination nodes have no outgoing links F Sn = ∅ if n ∈ D. Furthermore, the outgoing links of origins and the incoming links of destination nodes are called connectors. They operate under a very specific link model (defined by an FD with ja ≡ +∞ : {a ∈ tlo∈O , a ∈ hdd∈D }), which implies that no spillback will occur for those links as they have infinite storage capacity. The advantage of using connectors is that in output statistics (like total time spent, total density), vehicles that have left the origins but that are unable to enter the physical network (due to inflow constraints), will be accounted for spontaneously and hence contribute to overall performance measures of the network. The time axis is discretised into a number of intervals of equal duration τ . The entire simulation period’s time is defined as a multiple of these intervals such that T time slices τ i with time step index i = 1 . . . T are constructed. The algorithm (see Section 3) evaluates these time slices in upwind fashion starting with the earliest time interval. The demand between each OD-pair fod (i) is assumed to be stationary for each time interval [τi−1 , τi ). Similar, also the inflow ua (i) and outflow va (i) rates of each link are assumed stationary in between two time slices as directed by the node model. The formulation of routing principles is based on (Bellei et al., 2005; Nezamuddin and Boyles, 2014) description of implicit path enumeration. It consists of time varying turn proportions φ abd (i) specified for every turn from one incoming arc a ∈ A into an outgoing arc b ∈ A towards destination d ∈ D. The turn proportions directs the incoming flow that is heading towards the same destination vad (i) into the different outgoing links such that ubd (i ) = a∈BSn φabd (i ) · vad (i ). Note that as a result, an additional dimension is needed to store the information of destination specific flows. The turn proportions are exogenously given with b∈F Sn φabd (i ) = 1 for each incoming link, destination and time slice. The use of destination-based turn rates leads to the assumption that all vehicles driving towards the same destination are handled as identical drivers once merged into a link independently of the specific origin from which they entered the network. The routing information provided by the turn proportions is used within each time slice to construct turning fractions that are consistent with the observed arrival flow:



ab (i ) =

v (i ) · φabd (i ) d∈D vad (i )

ad

d∈D

(18)

The destination-based handling of the routing and the derivations of consistent turning fractions requires that the additional information of destination specific flow is propagated throughout the network. For this reason, the solution at a grid point is defined in terms of destination-specific cumulative flow. The destination-specific cumulative inflow and outflow of a link at the end of time interval i are respectively:

Uad (i ) =

i

uad (r ) · τ

(19)

vad (r ) · τ

(20)

r=1

Vad (i ) =

i

r=1

The core of the algorithm is the evaluation of the CVN functions along a kinematic wave, which corresponds to the evaluation of Ua (τ ) = d∈D Uad (τ ) or Va (τ ) = d∈D Vad (τ ) at a specific time instance, mostly situated in between two time slices. This evaluation requires an interpolation of the CVN functions; which is a simple linear interpolation because all uda and vda are stationary in between the time slices. The total number of time intervals it takes for a forward wave to reach the end of a link (rounded up to the nearest integer) is defined as ιa = la /γa · 1/τ and the interpolation fraction in the ι ·τ −l /γ last time interval is defined as ϕa = a τ a a . The total number of time intervals it takes for a backward wave to reach the a /ωa . start of a link is defined as κa = la /|ωa | · 1/τ and the interpolation fraction of this wave is defined as ψa = κa ·τ−l τ

2 Note that a subset of three parameters is sufficient to specify the triangular FD. The 4th parameter is always chosen such that the FD is consistent with a triangular shape.

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

11

4.2. Pseudo-code formulation Next, the base algorithm is presented. First consider only the text in bold of the pseudo-code in the algorithm below. This simplified implementation contains the core iterative procedure, without fast sweeping and changing the active set of node updates. Within each time slice, all nodes are updated until for each node the potential change n is below the precision level prDNL . The origin nodes are only updated during the first iteration. Destination nodes are also updated only once because they represent sinks with infinite storage from which no spillback originates (see also the FD of connectors). For this reason destination nodes are updated, after a solution has been found for all other nodes. Every node is still updated by executing the three steps identified in the previous section. First, sending/receiving flows and the turning fractions are calculated, next the node model is run, and finally the CVN values are updated. The CVN updates are computed for each destination and stored in a temporal intermediate variable Xad such that the result can be compared to the CVN value of the previous iteration to calculate n . The weight factors for calculating the part of n that is linked to changes on neighbouring nodes are already defined as the interpolation factors ϕ a and ψ a if ιa or κ a is equal to 1 and 0 otherwise. The weight for potential changes on the node is given by the variable δ n that is equal to 1 if the node model requires repeated evaluations to find a consistent solution and 0 otherwise. In all numerical experiments of Section 5 the node model is explicitly given and δn = 0 , ∀n ∈ N. The algorithm is completed (in plain text) by adding the procedure for managing the active set and fast sweeping over the active nodes. The active set Nact contains all nodes that are checked in each iteration. It is initialized by the time dependent initial set I(i) at the start of each new time slice. As argued in the last paragraph of the previous section, the initial set of a cold start consists of origin nodes for which flow departs (Eq. (21)), while for a warm start the initial set could also include other nodes. After each iteration the active set S is filled with nodes for which n > prDNL . The algorithm stops when this set is empty, which corresponds to the stopping criteria of the base algorithm. The initial set is appended with nodes for which CVN are changed and a future update might be necessary.

 I (i ) =

o∈O:



fod (i ) > 0 ,

∀i = 1 . . . T

(21)

d∈D

Algorithm overview. Pseudo-code of the entire DNL-procedure. for i = 1 . . . T . See Algorithm part 1. . S = I (i ) . n = Inf,∀n ∈ Nact ∩ N \ O . until n ≤ prDNL ,∀n ∈ N \{D ∪ O} do . . sort Nact in descending order of n . . for each n in Nact ∩ N \ {D ∪ O} do ... SFad = (1−ϕa )·Uad (i − ιa ) + ϕa ·Uad (i − ιa + 1 ). . . . . . −Vad (i − 1 ), ∀d ∈ D, ∀a ∈ BSn . . .SFa = max(0, min(ca , d∈D SFad )),∀a ∈ BSn . . .RFb = (1 − ψb ) · d∈DVbd (i − κb ) + ψb · d∈DVbd (i − κb + 1 ). . . . . . − d∈D Ubd (i − 1 ) + jb · lb , ∀b ∈ F Sn . . .RFb = max (0, min(cb ,RFb )), ∀b ∈ BSn (i )·φabd (i ) , (i )

∀a ∈ BSn , ∀b ∈ F Sn . . .{T Fab }∀a∈BSn ,∀b∈BSn = node_model ({SFa }∀a∈BSn ,{RFb }∀b∈F Sn ,,c ) . . .ab (i ) =

SF d∈D ad

d∈D SFad

. . .See Algorithm part 2. . . end . . S = {n ∈ N: n > prDNL } . end . See Algorithm part 3. end

(go over all time steps) (update all active origin nodes) (active set initialization) (initialize difference measure) (keep updating until precision is reached) (order set of active nodes) (update all active standard nodes) (sending flow calculation) (total sending flow calculation)

(total receiving flow calculation) (turning fraction calculation) (run the node model) (update CVN values & ࢞n ) (update active set) (update all active destination nodes)

Algorithm part 1. Pseudo-code for updating all (active) origin nodes. for each o in O ∩ I (i ) do . for each a in F So do . . Xad = Uad (i−1 ) + fod (i ),∀d ∈ D . . if d∈D |Xad − Uad (i )| > prDNL . . .Uad (i ) = Xad ,∀d ∈ D . . .add o to I (i + 1 ) . . .add hda to I (i + ιa − 1 ) . . end . end end

(go over all active origins) (go over all outgoing links) (calculate upstream CVN) (track upstream moving change) (update upstream CVN values) (add origin to initial set) (add downstream node to initial set)

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB 12

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

Algorithm part 2. Pseudo-code for updating CVN values of an (active) node and distributing changes towards neighbouring nodes. active = f alse Rbd = 0,∀b ∈ F Sn ,∀d ∈ D for each a in BSn do . Xad = Vad (i−1 ) + SFadSF · b∈F Sn T Fab ,∀d ∈ D ad d∈D . if d∈D |Xad − Vad (i )| > prDNL . . active = true . . add tla to I (i + κa − 1 ) . . if κa is equal to 1 . . .tla = tla + ψa · d∈D |Xad −Vad (i )| . . end . end . for each b in F Sn . . Rbd = Rbd + SFadSF · T Fab ,∀d ∈ D d∈D

ad

. end end for each b in F Sn . Xbd = Ubd (i−1 ) + Rbd ,∀d ∈ D . if d∈D |Xbd − Ubd (i )| > prDNL . . active = true . . add hdb to I (i + ιb − 1 ) . . if ιb is equal to 1 . . .hdb = hdb + ϕb · d∈D |Xbd −Ubd (i )| . . end . end end n = 0 if active . n = n + δn · a∈BSn d∈D |Xad −Vad (i )| . Vad (i ) = Xad , ∀a ∈ BSn , ∀d ∈ D . n = n + δn · b∈F Sn d∈D |Xbd −Ubd (i )| . Ubd (i ) = Xbd , ∀b ∈ F Sn ,∀d ∈ D . add n to I (i + 1 ) end end

(CVN are updated if a significant change is found) (initialize flow towards outgoing links) (go over all incoming links) (calculate downstream CVN) (track upstream moving change) (node needs a CVN update) (add upstream node to initial set) (only if iterations are required) (update  of upstream node)

(go over all outgoing links) (update flow towards outgoing links)

(go over all outgoing links) (calculate upstream CVN) (track downstream moving change) (node needs a CVN update) (add downstream node to initial set) (only if iterations are required) (update  of downstream node)

(reset changes affecting node) (only if node needs a CVN update) (update  of node) (update downstream CVN of incoming links) (update  of node) (update upstream CVN of outgoing links) (add node to initial set)

Algorithm part 3. Pseudo-code for updating all (active) destination nodes. for each d in D ∩ S do . for each a in BSd do . .Xad = (1−ϕa ) · Uad (i−ιa ) + ϕa · Uad (i−ιa + 1 ) . . if d∈D |Xad − Vad (i )| > prDNL . . .Vad (i ) = Xad . . .add d to I (i + 1 ) . . end . end end

(go over all active destinations) (go over all incoming links) (calculate downstream CVN) (track changes) (update downstream CVN) (add destination to initial set)

5. Application and case studies 5.1. Test networks In this section, numerical examples are used to verify the working of our algorithm. Four different test networks are selected that have different scales and morphologies. One is a very simple corridor network that consists of the interaction of three highways (R0–E40–E314) in the area between Leuven and Brussels, also presented in Himpe et al. (2013). The second is a more complex network covering the ring road around Rotterdam and the main arterials through the city center. The third network is based on the region around The Hague. It is a combination of a dense city center (with a rectangular grid street plan) and three large connecting highway corridors. The last is a medium-size regional network around Leuven that consists of all major local roads and the highway system east of Brussels. The different morphologies are clearly visible in Fig. 4. The different characteristics of each of the networks are presented in Table 1. Note that for all networks a 4 h peak period is simulated. The number of links and amount of zones (and more precisely the number of destinations) has a major impact on the memory requirements of the simulation (the total size of the solution is a function of T · |A| · |D|). The route choice in all the presented cases is limited to an all-or-nothing assignment based on the free flow travel times in an empty network. As a result the splitting rates of each time slice are the same. With φabd = 1 if that turn is part of the shortest path towards the destination (based on γ a ) and 0 otherwise. Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

13

Fig. 4. Different morphology of the considered networks. Table 1 Network characteristics. Network

# of nodes

# of links

# of zones

# of valid OD pairsa

Simulation period (h)

Total trips (excl. intra zone)

(1) (2) (3) (4)

39 331 1636 2581

60 562 3722 4653

9 44 160 417

42 1890 13,603 52,441

4 4 4 4

95,544 187,328 332,873 427,334

a

R0-E40-E314, Belgium Rotterdam, The Netherlands Regional The Hague, The Netherlands Regional Leuven, Belgium

The number of valid OD pairs represents the amount of OD-relations for which a positive demand exists.

5.2. Varying the simulation time interval In this section, the effect of increasing the time interval is investigated. In Fig. 5 a space–time diagram is shown of the densities along the congested E40 highway of the corridor network in case 1. It is clear that the contour of the congestion pattern for all results is similar. For small time steps, like 1 min and less, the resulting congestion pattern shows perturbation changes with a frequency in the range of the simulation interval. This is an indication of the increased sensitivity of the dynamics inside the queues. Small fluctuations in demand or supply are the source of upstream moving waves of increased density. For longer time intervals, a smoothing effect occurs, which makes the result less sensitive to these small changes while still providing an accurate contour of the congestion. The travel time for vehicles traveling through the traffic jam is also very similar for the different time intervals. The travel time of a link is computed as the difference between the entry time of a vehicle number as given by the cumulative inflow Ua and corresponding exit time of that same vehicle as found in the cumulative outflow Va . The integration of travel

Fig. 5. Congestion pattern under various time intervals; graph (a): 1 min, (b): 5 min, (c): 10 min.

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB 14

ARTICLE IN PRESS

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

Fig. 6. Travel time pattern of simulations with various time intervals for vehicles departing during the peak hour along the congested corridor of case 1.

time over a route leads unavoidably to interpolation errors as departure times in the next link do not correspond to time grid points and hence need to be interpolated. Alternatively, route travel times can be calculated based on the cumulative inflow/outflow of a route, if those CVN correspond to vehicles only traveling over the analysed route. In Fig. 6, the total travel time for vehicles traveling along the congested corridor is plotted for different departure times along the simulation period. The CVN values used for the analysis are extracted from a designated destination-layer for which traffic only flows along the specified route. The different travel time lapses show a reasonable fit. Compared to the results of small time intervals, the simulations with larger time intervals are biased towards lower total travel times. This is a consequence of the increased smoothing of flow that traverses multiple links in one time interval, which also results in shorter and lower saturation rates of potential bottlenecks. As can be seen, the total travel time based on route CVN is also sensitive to interpolation errors. Especially in free flow states, the errors might lead to unrealistic travel times (shorter than free flow travel time). These errors arise from the assumption that the link model is only enforced at the predefined grid points. Consider for example the first interval of a simulation with a large time interval. Assume that for each link a it holds that la /γ a < τ and iterations are required to find a consistent solution, but the time interval τ is shorter than total free flow travel time along the route. Now at the start of the simulation all CVN values Ua (τ 0 ) and Va (τ 0 ) are set equal to 0 and during the first interval some flow is send over the considered route. After convergence, all CVN values Ua (τ 1 ) and Va (τ 1 ) along that route will be nonzero. Such result implies that an infinitesimal number of vehicles departing right after τ 0 will reach the end of the route instantaneous. The errors disappear in stable conditions (constant flow rates over multiple time intervals). A quick fix would be to enforce a lower limit on travel times equal to the free flow travel time. However, interpolation errors will always persist in a rectangular space–time grid and the trade-off between accuracy and computation speed should be made case by case. 5.3. Convergence properties of the iteration scheme Now, the difference between the core algorithm and the different improvements, such as fast sweeping and active set managing, is analysed. For this, the solutions of test networks 1 and 2 are analysed with a 5 min time interval. The base algorithm updates all nodes in every iteration, so each iteration requires approximately the same computation time. In Fig. 7, this can be verified by checking the horizontal spacing between each point of the red line. Each triangularly shaped point represents the sum of the potential changes over all nodes at the end of an iteration. This sum needs to be small to reach convergence. The number of iterations before convergence increases for larger networks and increasing levels of congestion because more nodes are interacting. The horizontal spacing also increases for larger networks because more nodes are evaluated in every iteration (note that horizontal axis scale is different in the two sub-plots). Next, by ordering the nodes within each iteration and thus increasing efficiency by fast sweeping through the node updates, the convergence of the algorithm is accelerated. The blue line (diamond shaped points) in Fig. 7 is clearly positioned below the red line indicating a significant computation gain (take into account that the vertical axis has a logarithmic scale). The horizontal difference between two iterations is again relatively stable because all nodes are updated in each iteration. The computation gain of the fast sweeping procedure is more important if traffic states are related by kinematic waves that move in similar directions. This is the case in corridors with homogenous traffic states. Finally, the full algorithm with updating of the active set of nodes and fast sweeping is represented by the green line (with rectangular points). It is clear that by only updating a subset of nodes the efficiency is significantly increased especially near convergence. In the last iterations an ever smaller subset of nodes is still being updated because most nodes are already consistent with their neighbours, so horizontal distance (i.e. time per iteration) decreases rapidly, improving total Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB

ARTICLE IN PRESS W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

[m3Gsc;February 1, 2016;21:17] 15

Fig. 7. Total sum of potential differences of all nodes over iterations (time) for base algorithm and its extensions in case 1 and case 2 for a time interval in free flow and with some congestion.

Fig. 8. Total computation time for the entire simulation period (a) and average computation time for the calculation of a time slice (b) in each of the 4 test networks.

convergence speed. Clearly the results presented here are dependent on the specified network and the size of the time interval. Such influences are further investigated in the next section. Note that for the given precision bound, the average absolute difference between all three solutions is very similar and well within the range of the precision bound, so although the total number of node updates has decreased and the order of the nodes is different, the solution is still (numerically) the same. 5.4. Overview of performance Computation time and memory size are typically the largest resource bottlenecks for many large scale DNLs. The iterative approach allows for much larger time intervals to be used in the simulation. This is very effective for reducing the calculation time and coping with the available memory size as can be seen by looking at Fig. 8a. Calculations are performed in C++ on a basic laptop equipped with an Intel® Core i3 m380 CPU and 2GB of RAM memory. Note that both axes of the graph have a logarithmic scale. On the largest test networks, the analysis of the smaller time intervals was not possible because of Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB 16

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21 Table 2 Performance: total number of node updates. Network

R0-E40-E314, Belgium Rotterdam, The Netherlands Regional The Hague, The Netherlands Regional Leuven, Belgium

LTMa

35,567 2,047,129 16,108,924 15,687,402

Iterative LTM

τ = 6 s

τ = 60 s

τ = 600 s

50,283 577,296 – –

25,860 161,123 2,465,917 –

17,025 150,192 630,436 1,412,277

a These numbers are the theoretical minimum number of node updates of the original LTM. Node updating is different for each node and based on the most restrictive CFL condition of the connecting links. As a result, the total number of node updates is known beforehand and no actual simulation runs are needed, compared to the iterative scheme.

Fig. 9. Calculation time of warm vs cold start for increasing and decreasing trip matrix in the congested corridor of case 1.

the limited memory of our machine. A trend (dashed) line was added for a shorter simulation time period to approximate the computation time for these time intervals. The average computation time per time interval increases for larger time intervals (Fig. 8b). Yet, as fewer time intervals are required to cover the entire simulation time period, total computation time (Fig. 8a) decreases with increasing time interval. For the smaller networks a sigmoid shape is visible in Fig. 8b. For time intervals below the CFL bound, a single iteration is sufficient to reach a consistent solution. The average computation time for a time slice is in this case the same for any shorter time interval, because all nodes are updated only once. For large time intervals, a second plateau is observed. This indicates that also an absolute upper limit exists on the maximum number of iterations, needed to reach convergence. Alternatively, the performance can be expressed by the total number of node updates. This indicator is platform independent and can be used to compare our method with alternative discretization schemes. The computation time of the three steps required for a node update represents the core operations of a LTM-based scheme. It is important to emphasise that our scheme propagates destination specific information and that in such multi-commodity systems (Nezamuddin and Boyles, 2014; Yperman and Tampère, 2006) computation effort is larger for each node update compared to a single commodity system (Gentile 2010; Raadsen et al., 2014). In Table 2 we show the total number of node updates for each of the 4 cases, simulated with the original LTM and our iterative approach. Note that this measure includes iterative node updates of the same node within a time interval for the iterative LTM. From these results it is clear that the iterative scheme with large time intervals performs much less node iterations than the original LTM. The Regional Leuven network is constructed such that the CFL conditions for all links where maximized which is a clear benefit for the original LTM compared to the more fine-grained network structure of the region of The Hague. The Rotterdam case shows that although the number of node updates is similar for τ = 60 s and τ = 600 s there is a computational gain in Fig. 8 when working with a larger time interval, the additional node updates stem from the increased number of iterations for each time slice. This is a clear example of how implementation and data structures can influence the computation speed. In this case computing more iterations for a time slice has the benefit (of less overhead for memory access) compared to making less iterations and moving more from time slice to time slice (and data block to data block). Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB

ARTICLE IN PRESS W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

[m3Gsc;February 1, 2016;21:17] 17

Fig. 10. Left column top (a) to bottom (e): analysed route of test network 1, space–time plot of densities for the base scenario (0.43 s), space–time plot of differences in density for adding vehicles before the peak period 0.04 s (10.8×), during the peak period 0.11 s (3.9×), and after the peak period 0.03 s (14.3×). Right column top (f) to bottom (j): analysed route of test network 2, space–time plot of densities for the base scenario (9.54 s), space–time plot of differences in density for adding vehicles before the peak period 0.37 s (25.8×), during the peak period 1.1 s (8.7×), and after the peak period 0.26 s (36.7×).

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB 18

ARTICLE IN PRESS

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

5.6. Repeated calculations The analysis of the algorithm is concluded with the illustration of repeated DNL calculations to illustrate the advantage of warm starting. As presented in Section 3, in this case the iterative scheme is initialized by the result of a previously calculated DNL. Two cases will be used to present the efficiency gains possible with the iterative procedure. Note that the full algorithm (with active set updating and fast sweeping) is used for these scenarios. In all remaining case studies, a 5 min interval is used. For the first example, consider the network loading of test case 1 based on the trip matrix fod . The solution of a new DNL with a new trip matrix fod +  fod differs only marginally from the base network loading. In this specific case, demand is decreased or increased for each time slice of each non-zero cell of the original base trip matrix. If the new trip matrix is close to the old one, the computational gains are obviously the largest. In Fig. 9, this can be observed as all warm start scenarios take less computation time than the respective cold-started counter parts, with the largest gains for small deviations. There is a clear difference in computation gain for decreasing and increasing trip matrix. This is a result of the increased between-link interactions (spillback) in congested networks. A similar effect is observed in the general trend of the computation time of the cold-started network loadings, which is gradually increasing as the demand is increased and more congestion is present. The accuracy of the warm starting solution is sufficiently similar to that of a cold start: the average absolute differences between cold and warm start solutions are in the same order of magnitude as the precision bound. The upper bound on the error reaches a 10-fold of the average error. This is an important indicator when using warm start for numerical finite difference approximations (e.g. OD estimation, control/toll optimization). The finite difference should then be chosen such that observable differences are well above the numerical approximation ‘noise’, which according to this experiment is approximately one order of magnitude above the precision bound parameter. The next example shows that the largest computational gains are to be expected if the effect of a change is local. Now, only a single OD cell is modified such that flow along a single path is increased. Observe the congestion pattern (graphs b and g in Fig. 10) during the peak hour in test networks 1 and 2 along the paths highlighted in graphs a and b of Fig. 10. Next, 5 vehicles are added to the simulation before, during, and after the peak hour. The differences in density between the base scenario and the scenario with increased path flow are shown in the bottom graphs of Fig. 10. The computational gains, comparing to the base scenarios, are indicated between brackets above each graph. Note that these are now much larger than in the previous example. The largest gains are observed for adding vehicles near the end of the simulation period (graphs c and j). Remember that solutions are expressed in terms of CVN and adding vehicles at the beginning of the simulation period will trigger a cumulative change for the entire period, requiring more computations. Also congestion effects are influencing the result as they trigger more complex updating patterns for small changes (also note the different scales of the difference plots). Comparing the networks of test cases 1 and 2 it is clear that the computational gains are larger in bigger networks. The speed-ups observed on the largest two networks are even in the range of 100 up to 1000 times less computation time compared to the base scenario. 6. Conclusion and discussion An efficient iterative dynamic network loading algorithm is presented that is useful for the analysis of real world networks. Computation resources are a major concern for practical case studies that apply dynamic traffic models. Our implementation reduces the calculation effort of the network loading to a minimum by avoiding all unnecessary calculations. A first contribution solves the inefficient sampling of space and time by eliminating the need for small discretization steps for numerical stability. The only remaining limits are those of topological concerns (with homogenous stretches of road as the smallest building blocks of a network) and the time resolution required (for an accurate description of congestion) by the intended application. The second contribution expands the usefulness of dynamic network loadings for repeated analysis of a network. Evaluating traffic states under changing variables or adjusted parameter settings is typically required in many applications. The concept of a warm start, with the initialization of the algorithm with previous results, provides substantial computation gains. Both novelties make the iterative link transmission model a valuable candidate for providing practical DNL solutions. With the algorithmic improvements presented in this paper, we may be close to the lower limit of computational budget required for a DNL consistent with first order link propagation and intersection models. We hypothesize that significant further efficiency gain may only be possible at the cost of simplifying the theoretical assumptions behind those models, like for instance: link models without spillback (vertical queues), simplified intersection models that obviate the need for iterative solutions or approximations of the multi-commodity nature of the DNL by fixed time-dependent turning fractions. However, we feel that for many applications, the modelling assumptions used in our DNL might on the contrary require some further refinement rather than simplification. We briefly discuss two of them that we feel most crucial, both related to modelling of heavily loaded, yet non-saturated conditions. With the assumption in this paper of triangular fundamental diagram, no delay would be incurred by traffic on links unless there is congestion (saturation). The same holds for intersections: only when demand exceeds intersection supply (on average during the possibly long time increments of multiple minutes), will our DNL model predict delay because of saturation. The expected delay due to deterministic and stochastic short-lived queuing phenomena at nearly saturated junctions (Gerlough and Huber, 1975; Sections 8 and 9 on signalized Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

ARTICLE IN PRESS

JID: TRB

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

19

and unsignalized intersections) is hence neglected entirely. Our striving for long time steps even conflicts fundamentally with any traditional approach to this problem like explicitly simulating green–red phases or gap acceptance behaviour. Yet such non-saturation delays are important in DNL as by disregarding them, neither routes over motorway corridors nor urban arterials will exhibit realistic travel times and sensitivity of travel time for (non-saturated) traffic volumes. Remedies for both issues have been proposed in the past for LTM-based DNL: GLTM for concave fundamental diagrams (Gentile, 2010), and explicitly adding intersection delays (Yperman and Tampère, 2006). Integrating these two innovations consistently and efficiently in the newly developed iterative LTM algorithm of this paper, is therefore high on our priority list. Acknowledgments We would like to thank the three anonymous reviewers for their constructive comments. This study was supported by the Flemish Agency for Innovation and Technology (IWT/SBO 101684). Appendix Symbols N A G(N, A) O D tla hda la

γa ωa

ca ja BSn FSn

τi τ

T fod (i) φ abd (i) ua (i) va (i) uad (i) vad (i) Ua (i) Va (i) Uad (i) Vad (i)

ιa κa ϕa ψa n

Xad

δn

I(i) Nact

set of nodes set of links network graph set of origins set of destinations tail node of link a head node of link a length of link a free flow speed of link a spillback speed of link a capacity of link a jam density of link a set of incoming links of node n set of outgoing links of node n time slice with index i size of the time interval between consecutive time slices total number of time intervals within the simulation period demand between origin o and destination d for time interval [τi−1 , τi ] turn proportion form link a into link b towards destination d in time interval [τi−1 , τi ] inflow rate of link a in time interval [τi−1 , τi ] outflow rate of link a in time interval [τi−1 , τi ] inflow rate of link a in time interval [τi−1 , τi ] towards destination d outflow rate of link a in time interval [τi−1 , τi ] towards destination d accumulated inflow of link a up to time slice τ i accumulated outflow of link a up to time slice τ i accumulated inflow of link a up to time slice τ i towards destination d accumulated outflow of link a up to time slice τ i towards destination d number of time periods for the forward wave to traverse link a rounded up to the nearest integer number of time periods for the backward wave to traverse link a rounded up to the nearest integer interpolation fraction of the forward wave when it arrives in between two time slices interpolation fraction of the backward wave when it arrives in between two time slices potential for a change at node n temporal variable for storing cumulative values at intermediate steps in the algorithm weight of a change at node n in the potential change n of that same node initial set of active nodes for time slice τ i , listed for a potential update in the first iteration set of active nodes, listed for a potential update in the next iteration

References Adamo, V., Astarita, V., Florian, M., Mahut, M., Wu, J.H., 1999. Analytical modeling of intersections in traffic flow models with queue spill-back. In: Proceedings of the 15th Triennial Conference (ORCE), IFORS’99. Beijing, PR China. Ansorge, Rainer, 1990. What does the entropy condition mean in traffic flow theory? Transportation Research Part B 24 (2), 133–143. doi:10.1016/ 0191-2615(90)90024-S. Bellei, G., Gentile, G., Papola, N., 2005. A within-day dynamic traffic assignment model for urban road networks. Transportation Research Part B 39 (1), 1–29. Bliemer, Michiel C.J., Raadsen, MarkP.H., Smits, Erik-Sander, Zhou, Bojian, Bell, MichaelG.H., 2014. Quasi-dynamic traffic assignment with residual point queues incorporating a first order node model. Transportation Research Part B 68, 363–384. doi:10.1016/j.trb.2014.07.001.

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB 20

ARTICLE IN PRESS

[m3Gsc;February 1, 2016;21:17]

W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

Blumberg-Nitzani, Michal, Bar-Gera, Hillel, 2014. The effect of signalised intersections on dynamic traffic assignment solution stability. Transportmetrica A: Transport Science 10 (7), 622–646. doi:10.1080/18128602.2012.751682. Buisson, C., Lebacque, J.P., Lesort, J.B., 1996. Strada, a discretized macroscopic model of vehicular traffic flow in complex networks based on the godunov scheme. In: CESA’96 IMACS multiconference: computational engineering in systems applications, pp. 976–981. http://trid.trb.org/view.aspx?id=495534. Cascetta, Ennio, 2001. Transportation Systems Engineering: Theory and Methods. Springer Science & Business Media. Claudel, C.G., Bayen, A.M., 2010. Lax-Hopf based incorporation of internal boundary conditions into Hamilton–Jacobi equation. Part II: computational methods. IEEE Transactions on Automatic Control 55 (5), 1158–1174. doi:10.1109/TAC.2010.2045439. Claudel, C.G., Bayen, A.M., 2010. Lax-Hopf based incorporation of internal boundary conditions into Hamilton–Jacobi equation. Part I: theory. IEEE Transactions on Automatic Control 55 (5), 1142–1157. doi:10.1109/TAC.2010.2041976. Corthout, Ruben, Flötteröd, Gunnar, Viti, Francesco, Tampère, ChrisM.J., 2012. Non-unique flows in macroscopic first-order intersection models. Transportation Research Part B 46 (3), 343–359. doi:10.1016/j.trb.2011.10.011. Corthout, Ruben, Himpe, Willem, Viti, Francesco, Frederix, Rodric, Tampère, ChrisM.J., 2014. Improving the efficiency of repeated dynamic network loading through marginal simulation. Transportation Research Part C 41, 90–109. doi:10.1016/j.trc.2013.12.009. Corthout, Ruben, Tampère, Chris, Immers, Lambertus, 2009. Marginal incident computation. Transportation Research Record: Journal of the Transportation Research Board 2099, 22–29. doi:10.3141/2099-03. Courant, R., Friedrichs, K., Lewy, H., 1967. On the partial difference equations of mathematical physics. IBM Journal of Research and Development 11 (2), 215–234. doi:10.1147/rd.112.0215. Daganzo, Carlos F., 1994. The cell transmission model: a dynamic representation of highway traffic consistent with the hydrodynamic theory. Transportation Research Part B 28 (4), 269–287. Daganzo, Carlos F., 1995. The cell transmission model, part ii: network traffic. Transportation Research Part B 29 (2), 79–93. doi:10.1016/0191-2615(94) 00022-R. Daganzo, Carlos F., 2005a. A variational formulation of kinematic waves: solution methods. Transportation Research Part B 39 (10), 934–950. Daganzo, Carlos F., 2005b. A variational formulation of kinematic waves: basic theory and complex boundary conditions. Transportation Research Part B 39 (2), 187–196. doi:10.1016/j.trb.2004.04.003. de Palma, André, Lindsey, Robin, 2011. Traffic congestion pricing methodologies and technologies. Transportation Research Part C 19 (6), 1377–1399. doi:10. 1016/j.trc.2011.02.010. Flötteröd, Gunnar, Rohde, Jannis, 2011. Operational macroscopic modeling of complex urban road intersections. Transportation Research Part B 45 (6), 903–922. doi:10.1016/j.trb.2011.04.001. Frederix, Rodric, Viti, Francesco, Himpe, WillemW.E., Tampère, ChrisM.J., 2014. Dynamic origin–destination matrix estimation on large-scale congested networks using a hierarchical decomposition scheme. Journal of Intelligent Transportation Systems 18 (1), 51–66. doi:10.1080/15472450.2013.773249. Friesz, Terry L., Han, Ke, Neto, PedroA., Meimand, Amir, Yao, Tao, 2013. Dynamic user equilibrium based on a hydrodynamic model. Transportation Research Part B 47, 102–126. doi:10.1016/j.trb.2012.10.001. Gentile, G., Immers, L.G.H., Tampere, C.M.J., Viti, F., 2010. The general link transmission model for dynamic network loading and a comparison with the DUE algorithm. In: Transport Economics, Management and Policy Series. Edward Elgar Publishing, MA, USA. http://guidogentile.sistemaits.com/wp-content/ uploads/2010/04/dyn-3-52-GLTMandDUE-DTA2008.pdf. Gentile, Guido, 2015. Using the general link transmission model in a dynamic traffic assignment to simulate congestion on urban networks. Transportation Research Procedia 5, 66—81. doi:10.1016/j.trpro.2015.01.011, SIDT Scientific Seminar 2013. Gentile, Guido, Daniele, Tiddi, 2009. Synchronization of traffic signals through a heuristic-modified genetic algorithm with GLTM. In: Proceedings of the XIII Meeting of the Euro Working Group on Transportation. Padua, Italy: Padova University Press, pp. 1–9. https://iris.uniroma1.it/handle/11573/57762. Gerlough, Daniel L., Huber, MatthewJ., 1975. Traffic Flow Theory: A Monograph. Transportation Research Board, National Research Council. Gibb, John, 2011. Model of traffic flow capacity constraint through nodes for dynamic network loading with queue spillback. Transportation Research Record 2263, 113–122. doi:10.3141/2263-13. Godunov, S.K., 1959. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. Mat. Sb. (N.S.) 47 (89), 271–306. Greenberg, Harold, 1959. An analysis of traffic flow. Operations Research 7 (1), 79–85. doi:10.1287/opre.7.1.79. Greenshields, B.D., Bibbins, J.R., Channing, W.S., Miller, H.H., 1935. A study of traffic capacity. In: Highway Research Board Proceedings: Proceedings of the Annual Meeting. http://www.safetylit.org/citations/index.php?fuseaction=citations.viewdetails&citationIds[]=citjournalarticle_248667_38. Han, Ke, Gayah, Vikash V., 2015. Continuum signalized junction model for dynamic traffic networks: offset, spillback, and multiple signal phases. Transportation Research Part B: Methodological 77 (July), 213–239. doi:10.1016/j.trb.2015.03.005. Han, Ke, Piccoli, Benedetto, Szeto, W. Y., 2015. Continuous-time link-based kinematic wave model: formulation, solution existence, and well-posedness. Transportmetrica B: Transport Dynamics, 1–36. doi:10.1080/21680566.2015.1064793. Himpe, W.W.E., Corthout, R, Tampere, C.M.J., 2013. An implicit solution scheme for the link transmission model. In: 2013 16th International IEEE Conference on Intelligent Transportation Systems - (ITSC), pp. 572–577. doi:10.1109/ITSC.2013.6728292. Han, Ke, Benedetto, Piccoli, W.Y., Szeto., 2012. Continuous-time link based kinematic wave model: formulation, solution existence and well-posedness. arXiv:1208.5141 [math]. http://arxiv.org/abs/1208.5141. Holden, H., Risebro, N., 1995. A mathematical model of traffic flow on a network of unidirectional roads. SIAM Journal on Mathematical Analysis 26 (4), 999–1017. doi:10.1137/S0036141093243289. Jin, Wen-Long, 2015. Continuous formulations and analytical properties of the link transmission model. Transportation Research Part B 74, 88–103. doi:10. 1016/j.trb.2014.12.006. Jin, W.L., Zhang, H.M., 2003. On the distribution schemes for determining flows through a merge. Transportation Research Part B 37 (6), 521–540. doi:10. 1016/S0191-2615(02)00026-7. Laval, Jorge A., Leclercq, Ludovic, 2013. The Hamilton–Jacobi partial differential equation and the three representations of traffic flow. Transportation Research Part B 52, 17–30. doi:10.1016/j.trb.2013.02.008. Lebacque, Jean-Patrick, 2005. First-order macroscopic traffic flow models: intersection modeling, network modeling. In: Mahmassani, H.S. (Ed.), Proceedings of the 16th International Symposium on the Transportation and Traffic Theory. College Park, Maryland, USA, Elsevier, Oxford, pp. 365–386. http://trid. trb.org/view.aspx?id=758861. Lebacque, J.P., 1996. The godunov scheme and what it means for first order traffic flow models. In: Lesort, J.B. (Ed.), 13th ISTTT Symposium. Elsevier, New York, pp. 647–678. http://trid.trb.org/view.aspx?id=481282. Lighthill, Michael J., Whitham, GeraldBeresford, 1955. On kinematic waves. II. A theory of traffic flow on long crowded roads. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 229 (1178), 317–345. Mazaré, Pierre-Emmanuel, Dehwah, AhmadH., Claudel, ChristianG., Bayen, AlexandreM., 2011. Analytical and grid-free solutions to the Lighthill–Whitham– Richards traffic flow model. Transportation Research Part B 45 (10), 1727–1748. doi:10.1016/j.trb.2011.07.004. Moskowitz, K., 1965. Discussion of ‘freeway level of service as influenced by volume and capacity characteristics’. Highway Research Record 99, 43–44. Newell, Gordon F., 1993a. A simplified theory of kinematic waves in highway traffic, part I: general theory. Transportation Research Part B 27 (4), 281–287. Newell, Gordon F., 1993b. A simplified theory of kinematic waves in highway traffic, part III: multi-destination flows. Transportation Research Part B 27 (4), 305–313. Newell, Gordon F., 1993c. A simplified theory of kinematic waves in highway traffic, part II: queueing at freeway bottlenecks. Transportation Research Part B 27 (4), 289–303.

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013

JID: TRB

ARTICLE IN PRESS W. Himpe et al. / Transportation Research Part B 000 (2016) 1–21

[m3Gsc;February 1, 2016;21:17] 21

Nezamuddin, N., Boyles, Stephen D., 2014. A continuous DUE algorithm using the link transmission model. Networks and Spatial Economics 15 (3), 465–483. doi:10.1007/s11067-014-9234-x. Ni, Daiheng, Leonard, JohnD., Williams, BillyM., 2006. The network kinematic waves model: a simplified approach to network traffic. Journal of Intelligent Transportation Systems 10 (1), 1–14. doi:10.1080/15472450500455070. Raadsen, Mark P.H., Bliemer, Michiel C.J., Bell, Michael G.H., 2016. An efficient and exact event-based algorithm for solving simplified first order dynamic network loading problems in continuous time. Transportation Research Part B: Methodological. doi:10.1016/j.trb.2015.08.004, Accessed January 25. Richards, Paul I., 1956. Shock waves on the highway. Operations Research 4 (1), 42–51. Smits, Erik-Sander, Bliemer, MichielC.J., Pel, AdamJ., van Arem, Bart, 2015. A family of macroscopic node models. Transportation Research Part B 74, 20–39. doi:10.1016/j.trb.2015.01.002. Szeto, W.Y., Lo, H.K., 2005. Properties of dynamic traffic assignment with physical queues. Journal of the Eastern Asia Society for Transportation Studies 6, 2108–2123. Tampère, Chris M.J., Corthout, Ruben, Cattrysse, Dirk, Immers, LambertusH., 2011. A generic class of first order node models for dynamic macroscopic simulation of traffic flows. Transportation Research Part B 45 (1), 289–309. doi:10.1016/j.trb.2010.06.004. Treiber, Martin, Kesting, Arne, 2012. Traffic Flow Dynamics: Data, Models and Simulation. Springer, http://books.google.com/books?hl=en&lr= &id=sWiUjVAGTVgC&oi=fnd&pg=PR5&dq=%22focus+is+on+new+applications+ranging+from+novel+driver-assistance+systems,%22+%22modeling. +The+second+and+main+part+is+devoted+to+the%22+%22travel+times+(the+basis+of+dynamic+navigation),+and+how+to%22+&ots=CPerPei1io&sig= t3h9tTkLKc-sEXB1IciO3a4ZGqA. Underwood, R.T., 1961. Speed, volume and density relationships. Quality and Theory of Traffic Flow. Bureau of Highway Traffic, Yale University, New Haven, pp. 141–187. Wageningen-Kessels, Femke van, Lint, Hans van, Vuik, Kees, Hoogendoorn, Serge, 2014. Genealogy of traffic flow models. EURO Journal on Transportation and Logistics 4 (4), 445–473. doi:10.1007/s13676-014-0045-5. Weisstein, Eric W. 2015. Gauss–Seidel method (text accessed June 19). http://mathworld.wolfram.com/Gauss-SeidelMethod.html. Yperman, Isaak, Logghe, Steven, Immers, Ben, 2005. The link transmission model: an efficient implementation of the kinematic wave theory in traffic networks. Advanced OR and AI Methods in Transportation. Publishing House of Poznan University of Technology, Poznan, Poland, http://www.researchgate.net/publication/237532918_THE_LINK_TRANSMISSION_MODEL_AN_EFFICIENT_IMPLEMENTATION_OF_THE_KINEMATIC_ WAVE_THEORY_IN_TRAFFIC_NETWORKS/file/5046352a32c8519dc1.pdf. Yperman, Isaak, Steven, Logghe, Chris, M.J., Tampere, Ben, Immers., 2006. Multicommodity link transmission model for dynamic network loading. http: //trid.trb.org/view.aspx?id=776681. Yperman, Isaak, Chris M.J., Tampere, Ben, Immers., 2007. A kinematic wave dynamic network loading model including intersection delays. In. http://trid.trb. org/view.aspx?id=801218. Yperman, I., Tampère, Chris, 2006. Multi-commodity dynamic network loading with kinematic waves and intersection delays. In: Proceedings of the DTA 2006. Leeds, UK. https://lirias.kuleuven.be/handle/123456789/196093. Zhang, Yong-Tao, Zhao, Hong-Kai, Chen, Shanqin, 2006. Fixed-point iterative sweeping methods for static Hamilton–Jacobi equations. Methods and Applications of Analysis 13 (3), 299–320. Zhao, Hongkai, 2005. A fast sweeping method for eikonal equations. Mathematics of Computation 74 (250), 603–627. doi:10.1090/S0025-5718-04-01678-3.

Please cite this article as: W. Himpe et al., An efficient iterative link transmission model, Transportation Research Part B (2016), http://dx.doi.org/10.1016/j.trb.2015.12.013