Computing Dynamic User Equilibria on Large-Scale Networks with Software Implementation

Dynamic user equilibrium (DUE) is the most widely studied form of dynamic traffic assignment (DTA), in which road travelers engage in a non-cooperative Nash-like game with departure time and route choices. DUE models describe and predict the time-varying traffic flows on a network consistent with traffic flow theory and travel behavior. This paper documents theoretical and numerical advances in synthesizing traffic flow theory and DUE modeling, by presenting a holistic computational theory of DUE, which is numerically implemented in a MATLAB package. In particular, the dynamic network loading (DNL) sub-problem is formulated as a system of differential algebraic equations based on the Lighthill-Whitham-Richards fluid dynamic model, which captures the formation, propagation and dissipation of physical queues as well as vehicle spillback on networks. Then, the fixed-point algorithm is employed to solve the DUE problems with simultaneous route and departure time choices on several large-scale networks. We make openly available the MATLAB package, which can be used to solve DUE problems on user-defined networks, aiming to not only facilitate benchmarking a wide range of DUE algorithms and solutions, but also offer researchers a platform to further develop their own models and applications. The MATLAB package and computational examples are available at https://github.com/DrKeHan/DTA.


Introduction
This paper is concerned with a class of models known as dynamic user equilibrium (DUE). DUE problems have been studied within the broader context of dynamic traffic assignment (DTA), which is viewed as the modeling of time-varying flows on traffic networks consistent with established travel demand and traffic flow theory.
DTA models, from the early 1990s onward, have been greatly influenced by Wardrop's principles (Wardrop 1952).
• Wardrop's first principle, also known as the user optimal principle, views travelers as Nash agents competing on a network for road capacity. Specifically, the travelers selfishly seek to minimize their own travel costs by adjusting route choices. A user equilibrium is envisaged where the travel costs of all travelers in the same origin-destination (O-D) pair are equal, and no traveler can lower his/her cost by unilaterally switching to a different route. • Wardrop's second principle, known as the system optimal principle, assumes that travelers behave cooperatively in making their travel decisions such that the total travel cost on the entire network is minimized. In this case, the travel costs experienced by travelers in the same O-D pair are not necessarily identical.
Since the seminal work of Merchant and Nemhauser (1978a, b) the DTA literature has been focusing on the dynamic extension of Wardrop's principles, which gives rise to the notions of dynamic user equilibrium (DUE) and dynamic system optimal (DSO) models. The DUE model stipulates that experienced travel cost, including travel time and early/late arrival penalties, is identical for those route and departure time choices selected by travelers between a given O-D pair. The DSO model seeks system-wide minimization of travel costs incurred by all the travelers subject to the constraints of fixed travel demand and network flow dynamics.
For a comprehensive review of DTA models, the reader is referred to Boyce et al. (2001); Peeta and Ziliaskopoulos (2001); Lo (2005, 2006); Jeihani (2007) and Bliemer et al. (2017). A discussion of DTA from the perspective of intelligent transportation system can be found in Ran and Boyce (1996a). Chiu et al. (2011) present a primer on simulation-based DTA modeling. Garavello et al. (2016) focus on the traffic flow modeling aspect of DTA, namely the hydrodynamic models for vehicular traffic and their network extensions. Wang et al. (2018) review relevant DTA literature concerning environmental sustainability.
Dynamic user equilibrium (DUE), which is one type of DTA, remains a modern perspective on traffic network modeling that enjoys wide scholarly support. It is conventionally studied as a open-loop, non-atomic Nash-like game (Friesz et al. 1993).
The notion of open loop refers to the assumption that the travelers' route choices do not change in response to dynamic network conditions after they leave the origin. The non-atomic nature refers to the prevailing technique of flow-based modeling, instead of treating the traffic as individual vehicles; this is in contrast to agent-based modeling (Balmer et al. 2004;Cetin et al. 2003;Shang et al. 2017). DUE captures two aspects of travel behavior quite well: departure time choice and route choice. Within the DUE model, travel cost for the same trip purpose is identical for all utilized path and departure time choices. The relevant notion of travel cost is a weighted sum of travel time and arrival penalty.
On the other hand, the DNL sub-problem captures the relationship between dynamic traffic flows and travel delays, by articulating link dynamics, junction interactions, flow propagation, link delay, and path delay. It has been the focus of a significant number of studies in traffic modeling (Ran and Boyce 1996b;Xu et al. 1999;Tong and Wong 2000;Lo and Szeto 2002;Wie et al. 2002;Szeto and Lo 2004;Yperman et al. 2005;Perakis and Roels 2006;Nie and Zhang 2010;Ban et al. 2011;Ukkusuri et al. 2012;Han et al. 2013a). The DNL model aims at describing and predicting the spatial and temporal evolution of traffic flows on a network that is consistent with established route and departure time choices of travelers. This is done by introducing appropriate dynamics to flow propagation, flow conservation, link delay, and path delay on a network level. Any DNL must be consistent with the established path flows, and is usually performed under the first-in-first-out (FIFO) principle (Szeto and Lo 2006). In general, DNL models have the following components: 1. some form of link and/or path dynamics; 2. an analytical relationship between flow/speed/density and link traversal time 3. flow propagation constraints; 4. a model of junction dynamics and delays; 5. a model of path traversal time; and 6. appropriate initial conditions. DNL gives rise to the path delay operator, which is analogous to the pay-off function in classical Nash games, and plays a pivotal role in DTA and DUE problems. The properties of the delay operator are critical to the existence and computation of DUE models. However, it is widely recognized that the DNL model or the delay operator is not available in closed form; instead, they have to be numerically evaluated via a computational procedure. As a result, the mathematical properties of the delay operator remain largely unknown. This has significantly impacted the computation of DUE problems due to the lack of provable convergence theories, which all require certain form of generalized monotonicity of the delay operators. Table 1 shows some relevant computational algorithms for DUE and their convergence conditions with respect to the continuity and monotonicity of the delay operators. The reader is referred to Han et al. (2015a) for definitions of different types of generalized monotonicity. This paper documents theoretical and numerical advances in synthesizing traffic flow theory and traffic assignment models, by presenting a computational theory of DUE, which includes algorithms and software implementation. While there have been numerous studies on the modeling and computation of DUEs, including those reviewed in this paper, little agreement exists regarding an appropriate mathematical formulation of DUE or DNL models, as well as the extent to which certain models should/can be applied. This is partially due to the lack of open-source solvers and a set of benchmarking test problems for DUE models. In addition, large-scale computational examples of DUE were rarely reported and, when they were, little detail was provided that allows the results to be validated, reproduced and compared.
This paper aims to bridge the aforementioned gaps by presenting a computable theory for the simultaneous route-and-departure-time (SRDT) choice DUE model along with open-source software packages. In particular, the DNL model is based on the Lighthill-Whitham-Richards fluid dynamic model (Lighthill and Whitham 1955;Richards 1956), and is formulated as a system of differential algebraic equations (DAEs) by invoking the variational theory for kinematic wave models. This technique allows the DAE system to be solved with straightforward time-stepping without the need to solve any partial differential equations. Moreover, this paper presents the fixed-point algorithm for solving the DUE problem, which was derived from the differential variational inequality formalism . Both the DNL procedure and the fixed-point algorithm are implemented in MATLAB, and In addition, we make openly available the MATLAB package, which can be used to solve DUEs on user-defined networks. The package is documented in this paper (Appendix), and the codes and examples are available at https://github. com/DrKeHan/DTA. It is our intention that the open-source package will not only help DTA modelers with benchmarking a wide range of algorithms and solutions, but also offer researchers a platform to further develop their own models and applications.
The rest of the paper is organized as follows. Section 2 introduces some key notions and mathematical formulations of DUE. Section 3 details the dynamic network loading procedure and the DAE system formulation. The fixed-point algorithm for computing DUE is presented in Section 4. The computational results obtained from the proposed DUE solver are presented in Section 5. Section 6 offers some concluding remarks. Finally, the MATLAB software package is documented in the Appendix.

Formulations of Dynamic User Equilibrium
We introduce a few notations and terminologies for the ease of presentation below.
We define the effective delay operator as follows: : Here, the term 'effective delay' (Friesz et al. 1993) is a generalized notion of travel cost that may include not only a linear combination of travel time and arrival penalty, but also other forms of costs such as road pricing. The effective delay operator is essential to the DUE model as it encapsulates the physics of the traffic network by capturing the dynamics of traffic flows at the link, junction, path, and network levels.
We will discuss this operator in greater detail in Section 3. The travel demand satisfaction constraint is expressed as Therefore, the set of feasible path departure vector can be expressed as The following definition of dynamic user equilibrium is first proposed by Friesz et al. (1993) in a measure-theoretic context.

Definition 2.1 (SRDT DUE) A vector of departures h * ∈ is a dynamic user equilibrium with simultaneous route and departure time (SRDT) choice if
where 'a.e.', standing for 'for almost every', is a technical term employed by measure-theoretic arguments to indicate that Eq. (2.4) only needs to hold in [t 0 , t f ] with the exception of any subset that has zero Lebesgue measure.

Variational Inequality Formulation of DUE
Using measure-theoretic arguments, Friesz et al. (1993) establish that a SRDT DUE is equivalent to the following variational inequality under suitable regularity conditions: The VI (2.5) may be written in a more generic form by invoking the inner product · in the Hilbert space L 2 [t 0 , t f ] |P| : This leads to the following VI representation of DUE:

Nonlinear Complementarity Formulation of Due
The variational inequality formulation (2.5) is equivalent to the following nonlinear complementarity problem Ukkusuri et al. 2012): where a ⊥ b means that a · b = 0. While the equivalence between the VI and complementarity formulations is obvious in a discrete-time setting, the continuous-time version requires a non-trivial measure-theoretic argument. This is yet to be rigorously proven.

Differential Variational Inequality Formulation of DUE
It is noted in (Friesz et al. 2011) that the VI formulation (2.5) of DUE is equivalent to a differential variational inequality (DVI). This is most easily seen by noting that the demand satisfaction constraints may be re-stated as which is recognized as a two-point boundary value problem. The DUE may be consequently expressed as a DVI: The equivalence of DUE and DVI (2.9) can be shown with elementary optimal control theory applied to a linear-quadratic problem as demonstrated in Friesz et al. (2011) and Friesz and Han (2018). The DVI formulation is significant because it allows the still emerging theory of differential variational inequalities (Friesz 2010;Pang and Stewart 2008) to be employed for the analysis and computation of solutions of the DUE problem when simultaneous departure time and route choices are within the purview of users (Friesz and Meimand 2014;Friesz and Mookherjee 2006;Han et al. 2015a). However, the emerging literature on abstract differential variational inequalities has not been well exploited for either modeling or computing simultaneous route-and-departure-time choice equilibria; this gap was recently bridged by Friesz and Han (2018).

Fixed-point Formulation of DUE
Experience with variational inequalities suggests that there exists a fixed-point re-statement of the DUE problem in a proper functional space. The fixed-point formulation of continuous-time DUE is first articulated in Friesz et al. (2011) using optimal control theory. We define P 0 [·] to be the minimum-norm projection opera- Then the following fixed-point problem is equivalent to the DUE problem: where α > 0 is a fixed constant. The equivalence result may be proven by applying the minimum principle and associated optimality conditions to the minimum-norm problem intrinsic to the projection operator P 0 [·]. See Friesz (2010) for the details suppressed here.

Dynamic Network Loading
An integral component of the DUE formulation is the effective delay operator , which is constructed using the dynamic network loading (DNL) procedure. This section details one type of DNL models based on the fluid dynamic approximation of traffic flow on networks, known as the Lighthill-Whitham-Richards (LWR) model (Lighthill and Whitham 1955;Richards 1956). The model, including its various forms (Newell 1993a;Daganzo 1994Daganzo , 1995Yperman et al. 2005;Osorio et al. 2011), is widely studied in the DTA literature. This section presents a complete DNL procedure based on the LWR model and its variational representation, which lead to a differential algebraic equation (DAE) system.

The Lighthill-Whitham-Richards Link Model
The LWR model is capable of describing the physics of kinematic waves (e.g. shock waves, rarefaction waves), and allows network extensions that capture the formation and propagation of vehicle queues as well as spillback. It describes the spatial and temporal evolution of vehicle density ρ(t, x) on a road link using the following partial differential equation: where the link of interest is represented as a spatial interval [a, b]. The fundamental diagram f (·) is a continuous and concave function of density ρ, and satisfies f (0) = f (ρ jam ) = 0 where ρ jam denotes the jam (maximum) density. Furthermore, there exists a unique critical density value ρ c at which f (·) attains its maximum f (ρ c ) = C were C denotes the flow capacity of the link. A few widely adopted forms of f (·) include the Greenshields (Greenshields 1935), the trapezoidal (Daganzo 1994(Daganzo , 1995, and the triangular (Newell 1993a, b, c) fundamental diagrams. In the remainder of the paper (and also in the MATLAB package), we consider the following triangular fundamental diagram: where v > 0 and w > 0 denote the forward and backward kinematic wave speeds, respectively. While (3.1) captures the within-link dynamics, the inter-link propagation of congestion requires a careful treatment of junction dynamics, which is underpinned by the notions of link demand and supply.

Link Demand and Supply
We consider a road junction with m incoming links and n outgoing links. The dynamic on each of the m + n links is governed by the LWR model (3.1). Understandably these m + n equations are coupled via their relevant boundary conditions. In particular, the following flow conservation constraint must hold: where, without causing any confusion, we always use subscript i or j to indicate the association with link i or j . Equation (3.3) simply means that the total flow through the junction is conserved. However, this condition alone does not guarantee a unique flow profile on these m + n links, and additional conditions need to be imposed (Garavello et al. 2016). To this end, we define the link demand and supply (Lebacque and Khoshyaran 1999), where the demand (supply) is viewed as a function of the density near the exit (entrance) of the link: Intuitively, the demand (supply) indicates the maximum flow that can exit (enter) the link. That is, . . , n}. Similar to Eq. (3.3), Eq. (3.6) ensures the physical validity of the flows through the junction. Nevertheless, additional conditions are still needed to isolate a unique flow profile at this junction; these conditions are often derived based on driving behavior or traffic management measures, such as flow distribution (Daganzo 1995), right of way (Daganzo 1995;Coclite et al. 2005;Jin 2010), and traffic signal control (Han et al. 2014;Han and Gayah 2015). A review of these different models is provided in Garavello et al. (2016).

The Variational Representation of Link Dynamics
The variational solution representation of Hamilton-Jacobi equations has been widely applied to the LWR-type traffic modeling (Newell 1993a, b, c;Daganzo 2005Daganzo , 2006Claudel and Bayen 2010;Laval and Leclercq 2013;Costeseque and Lebacque 2014;Han et al. 2017). Our paper follows this approach. In particular, we consider the Moskowitz function (Moskowitz 1965), N(t, x), which measures the cumulative number of vehicles that have passed location x along a link by time t. The following identities hold: where 'almost everywhere' indicates possible spatial discontinuities in the density profile, i.e. shockwaves. It is straightforward to show that N(t, x) satisfies the following Hamilton-Jacobi equation: (3.7) Next, we denote by f in (t) and f out (t) the link inflow and outflow, respectively. The cumulative link entering and exiting vehicle counts are defined as where the superscripts 'up' and 'dn' represent the upstream and downstream boundaries of the link, respectively. Han et al. (2016b) derive explicit formulae for the link demand and supply based on a variational formulation known as the Lax-Hopf formula (Aubin et al. 2008;Claudel and Bayen 2010), as follows: where L = b − a denotes the link length. Note that Eqs. (3.8) and (3.9) express the link demand and supply, which are inputs of the junction model, in terms of N up (·), N dn (·), f in (·) and f out (·). This means that one no longer needs to compute the dynamics within the link, but to focus instead on the flows or cumulative counts at the two boundaries of the link. This tremendously simplifies the link dynamics, and gives rise to the link-based formulation (Jin 2015;Han et al. 2016b) or, in its discrete form, the link transmission model (Yperman et al. 2005).

Junction Dynamics that Incorporate Route Information
Essential to the network extension of the LWR model are the junction dynamics. Unlike many existing junction models such as those reviewed in Section 3.2, in a path-based DNL formulation, one must incorporate established routing information into the junction model. Such information is manifested in an endogenous flow distribution matrix, which specifies the proportion of exit flow from a certain incoming link that advances into a given outgoing link. This can be done by explicitly tracking the route composition in every unit of flow along the link. We begin by defining the link entry time function τ (t) where t denotes the corresponding exit time. Such a function can be obtained by evaluating the horizontal difference between the cumulative curves N up (·) and N dn (·) (Friesz et al. 2013). Next, for link i and path p such that i ∈ p, we define μ p i (t, x) to be the percentage of flow on link i that belongs to path p at any point (t, x). The first-in-first-out principle yields the following identity: (3.10) Now we consider a junction J with incoming links labeled as i ∈ {1, . . . , m} and outgoing links labeled as j ∈ {1, . . . , n}. The distribution matrix A J (t) can be expressed as Here, the scalar α ij (t) represents the proportion, among traffic exiting link i, that advances to link j . There exist a number of choices for the junction models, all of which need to satisfy the flow conservation constraint (3.3), the demand-supply constraints (3.6), and depend on the flow distribution matrix A J (t). Any such model can be conceptually expressed as: where D i (t), S j (t) and A J (t) are treated as input arguments of the junction model. The output, shown as the left hand side of (3.11), include the outflows (inflows) of the incoming (outgoing) links. The mapping for simple merge or diverge junction is illustrated by Daganzo (1995). In the conservation law literature is sometimes referred to as the Riemann Solver, and we refer the reader to Lebacque and Khoshyaran (1999), Coclite et al. (2005), Jin (2010), Han and Gayah (2015), Han et al. (2016b), and Garavello et al. (2016) for more discussion.

Dynamics at the Origin Nodes
A model at the origin (source) nodes is needed since the path flows h p (·), defined by Eq. (2.3), are not bounded from above. In this case, a queuing model is needed at the origin node to accommodate departure flow exceeding the supply of relevant downstream link. We employ a simple point-queue type dynamic for the origin node o. Denote by q o (t) the size of the point queue, and let link j be the link connected to the origin node. We have that where P o denotes the set of paths originating from o. The first term on the right hand side of Eq. (3.12) represents flow into the point queue, while the second term represents flow leaving the queue, where the demand at the origin is defined as and M is a sufficiently large number, e.g. larger than the flow capacity of link j . We note that the proposed queuing dynamics differ from the classical point-queue model (Vickrey 1969;Ban et al. 2011;Han et al. 2013a) in that our model assumes varying queue exit capacity, as constrained by the downstream supply S j (t), instead of a fixed flow capacity.

Calculating Path Travel Times
The DNL procedure calculates the path travel times for a given set of path departure rates. The path travel time consists of link travel times plus possible queuing time at the origin. We define the link exit time function λ(t) by measuring the horizontal difference between the cumulative entering and exiting counts: We note that Eq. (3.13) is also applied to calculate the queue exit time function at the origin following the dynamics presented in Section 3.5. For a path expressed as p = {1, 2, . . . , K}, the path exit time for a given departure time t is calculated as where f • g(t) .
= g(f (t)) denotes the composition of two functions. λ o (·) is the exit time function for the potential queuing at the origin o.

The Differential Algebraic Equation System Formulation of DNL
As a summary of the individual sections presented so far, we formulate the complete system of differential algebraic equations (DAEs). We begin with the following list of key notations.
supply of link i μ p i (t, x) percentage of flow on link i that belongs to path p q o (t) point queue at the origin node o ∈ S τ i (t) entry time of link i corresponding to exit time t λ i (t) exit time of link i corresponding to entry time t The DAE system reads: Equations (3.15)-(3.23) form the DAE system for the DNL procedure. Compared to the partial differential algebraic equation (PDAE) system presented in Han et al. (2016a), the DAE system does not involve any spatial derivative as one would expect from the LWR-type equations, by virtue of the variational formulation. The proposed DAE system may be time-discretized and solved in a forward fashion. Figure 1 explains the time-stepping logic, which is implemented in the Matlab package alongside this publication.

The Fixed-Point Algorithm for Computing DUE
The computation of DUE is facilitated by the equivalent mathematical formulations such as variational inequality, differential variational inequality, fixed-point problem, and nonlinear complementarity problem. Some solution algorithms and their convergence conditions are mentioned in Table 1. In this section, we present the following algorithm based on the fixed-point formulation (2.11): where α > 0 is a constant representing the step size, h k+1 and h k respectively represent the path departure rate vector at the (k + 1)-th and k-th iterations; (h k ) denotes the effective path delay vector. Recalling the definition of 0 from Eq. (2.10), the right hand side of Eq. (4.1), which involves projection, amounts to a linear-quadratic optimal control problem whose dual variables can be explicitly found following Pontryagin's minimum principle. The reader is referred to Friesz et al. (2011) for detailed discussion. Fig. 1 Time-stepping logic of the discretized DAE system (3.15)-(3.23). Here, t = 0, 1, 2 · · · N , represents discrete time steps with step size denoted t We outline the main steps of the fixed-point algorithm below.

Remark 4.1
In the fixed-point algorithm, the critical step is to find the dual variable v ij in (4.2). Note that this amounts to finding x such that G(x) = 0, where is a continuous function with a single argument x. Therefore, v ij can be found via standard root-finding algorithms such as Bracketing or Bisection methods.

Computational Examples of DNL and DUE
We present several computational examples of the simultaneous route-and-departuretime choice dynamic user equilibria on four networks of varying sizes and shapes, as illustrated in Table 2 and Fig. 2. In particular, the Nguyen network was initially studied by Nguyen (1984), and the other three networks are based on real-world cities in the US, although different levels of network aggregation and simplification   Fig. 2 The four test networks have been applied. Detailed network parameters, including coordinates of nodes and link attributes, are sourced and adapted from Transportation Networks for Research Core Team (2018). Given that the DUE and DNL formulations in this paper are pathbased, enumeration of paths was applied to generate path sets using the Frank-Wolfe algorithm.
The time horizon of analysis is [0, 5] (hrs), and time is in hr. The effective path delays (in hr) are defined for all four networks as follows: where we employ asymmetric quadratic arrival penalties with TA ij being the target arrival time of O-D pair (i, j ). We apply the fixed-point algorithm (Section 4) with the embedded DNL procedure (Section 3). The fixed-point algorithm is chosen among many other alternatives in the literature, as our extensive experience with DUE computation suggests that this method tends to exhibit satisfactory empirical convergence within limited number of iterations, despite that its theoretical convergence requires strong monotonicity of the delay operator. The DNL sub-model is solved as a DAE system (3.15)-(3.23), following the time stepping logic in Fig. 1. All the computations reported in this section were performed using the MATLAB (R2017b) package on a standard desktop with Intel i5 processor and 8 GB of RAM.

Performance of the Fixed-Point Algorithm
The termination criterion for the fixed-point algorithm is set as follows: where h k denotes the path departure vector in the k-th iteration. The threshold is set to be 10 −4 for the Nguyen and Sioux Falls networks, and 10 −3 for the Anaheim and Chicago Sketch networks. These different thresholds were chosen to accommodate the varying convergence performances of the algorithm on different networks (also see Fig. 3).   Table 3 summarizes the computational performance of the fixed-point algorithm and DNL procedure on different networks based on the termination criterion (5.1). It is shown that the number of fixed-point iterations is not significantly impacted by the varying sizes of the test networks, which suggests the scalability of the algorithms. Figure 3 shows the relative gaps, i.e. left hand side of Eq. (5.1), for a total of 100 fixed-point iterations on the four networks. It can be seen that for relatively small networks (Nguyen and Sioux Falls), the convergence can be achieved relatively quickly and to a satisfactory degree; and the corresponding curves are monotonically decreasing and smooth. For Anaheim and Chicago Sketch networks, the decreasing trend of the gap can sometimes stall and experience fluctuations locally.

DUE Solutions
In this section we examine the DUE solutions obtained upon convergence of the fixed-point algorithm. We begin by selecting four arbitrary paths per network to illustrate the properties of the solutions. Figure 4 shows the path departure rates as well as the corresponding effective path delays. We observe that the departure rates are nonzero only when the corresponding effective delays are equal and minimum, which conforms to the notion of DUE. Note that the bottoms of the effective delay curves should theoretically be flat, indicating equal travel costs. This is not exactly the case in the figures since we can only obtain approximate DUE solutions in a numerical sense, given the finite number of fixed-point iterations performed to reach those solutions. To rigorously assess the quality of the DUE solutions, we define the gap function between each O-D pair (i, j ) ∈ W as Here, GAP ij represents the range of travel costs experienced by travelers in the O-D pair (i, j ). In an exact DUE solution, the gap should be zero for all O-D pairs. Figure 5 summarizes all the O-D gaps of the DUE solutions on the four test networks. It can be seen that the majority of the O-D gaps are within 0.2 (hrs) across all four networks. Even for large-scale networks (Anaheim and Chicago Sketch), the 75th percentiles are within 0.16, and the whiskers extend to 0.3. A comparison between Anaheim and Chicago Sketch also reveals that the latter yields better solution quality in terms of O-D gaps, despite the obviously larger size of the problem. This suggests that the solution quality is not necessarily compromised by the size of the network.

Conclusion
This paper presents a computational theory for dynamic user equilibrium (DUE) on large-scale networks. We begin by presenting a complete and generic dynamic network loading (DNL) procedure based on the network extension of the LWR model and the variational theory, which allows us to formulate the DNL problem as a system of differential algebraic equations (DAEs). The DNL model is capable of capturing the formation, propagation and dissipation of physical queues as well as vehicle spillback. The DAE system can be discretized and solved in a time-forward fashion. In addition, the fixed-point algorithm for solving DUE problems is presented. Both the DAE system and the fixed-point algorithm are implemented in MAT-LAB, and the programs are developed in such a way that they can be applied to solve DUE and DNL problems on any user-defined networks. The MATLAB package is documented in the Appendix of this paper, which provides detailed instructions for using the solvers.
The MATLAB program is applied to solve DUE problems on several test networks of varying sizes. The largest one is the Chicago Sketch network with 86,179 O-D pairs and 250,000 paths. To the authors' knowledge this is by far the largest instance of SRDT DUE solution reported in the literature. Hopefully, our efforts in making these codes and data openly accessible could facilitate the testing and benchmarking of dynamic traffic assignment algorithms, and promote synergies between model development and applications.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A1. Dynamic Network Loading Solver
The MATLAB solver for the dynamic network loading problem is based on the timediscretized DAE system following the procedure in Fig. 1. The individual scripts (.m files) and input data files are illustrated in Fig. 6.
The DNL solver contains three main script files:

Pre-processor: 'CREATE NETWORK PROPERTIES.m'
The purpose of this script is to prepare all the variables required for the main dynamic network loading script and save them to file for reuse to save computational time. Follow the steps below to run this script. 1 The DNL code employs junction models that incorporate the continuous signal model proposed by Han et al. (2014) and Han and Gayah (2015). For any given junction (node), each of its incoming links, including the source if the node happens to be an origin, will be assigned a priority parameter between 0 and 1, such that the sum of relevant priority parameters equals 1. Once the priority of the source is specified, the priorities of the remaining incoming links are set proportional to the links' capacities. These parameters may be changed dynamically to accommodate different signal control scenarios or strategies. 2.3. The output file will be saved to the working directory with ' out.mat' suffix, containing data required for graphical display.

Graphical display: 'OPEN NETWORK GUI.m'
3.1. Run script (can be run after loading any saved output files into the workspace) Remark A.1 While the main DNL program is written as a script file here, it can be easily modified into a function file. This is the case in our DUE solver presented later. Fig. 7 The Braess network (example) Fig. 8 Network data for the Braess network (example). Note that, for simplicity, the MATLAB program assumes that v = 3w where v and w are respectively the speeds of the forward and backward kinematic waves; see Eq. 3.2. Therefore, the link attributes listed in this figure will identify a unique triangular fundamental diagram Example network data, path data, and path departure vector are illustrated in Figs. 8, 9 and 10, respectively.
To run the dynamic network loading program for the Braess network, we follow the diagram in Fig. 6  On the lower-left corner of the user interface, one can specify the path number and press the "Path travel time" button. The interface then shows the path on the network and the corresponding path travel time function; see Fig. 12 for an example.

A3. Dynamic User Equilibrium Solver
The dynamic user equilibrium solver is based on the fixed-point algorithm presented in Section 4. The individual scripts or functions (.m files) and input data files are illustrated in Fig. 13.
The MATLAB implementation of the DUE solution procedure consists of three main steps:

Pre-processing the O-D information: 'PROCESS OD INFO.m'
The purpose of this script is to process the network data (e.g. 'Braess8 pp.mat') in order to extract information about all the O-D pairs and their association with the paths. Follow the steps below to run this script.
1.1. Prepare the pre-processed network data (e.g. 'Braess8 pp.mat') and place it in the same working directory as 'PROCESS OD INFO.m'. 1.3. The output file, named 'OD info.mat', will be saved to the working directory for use by the DUE solver.

Specify input file 'Network planning parameters.mat'
This step mainly specifies the O-D demand and target arrival times. The O-D demand data 'OD demand' is a vector of positive numbers, whose dimension is equal to the number of O-D pairs in the network. The target arrival time data 'T A' is a vector of target arrival times, whose dimension is also equal to the number of O-D pairs. This means that each O-D pair is associated with one target arrival time.

Main program: 'MAIN.m'
3.1. Set user inputs (at the top of file): a) Time step (in second) of the discretized problem. Note that ideally the time step should be no less than the minimum link free-flow time in the network to ensure numerical stability. However, the code is conditioned to accommodate violation of this rule, in which case a warning message will be displayed. b) Threshold that serves as the convergence criterion of the fixedpoint algorithm; see Eq. 5.1. c) Tolerance (indifference band) if bounded rationality (BR) is to be considered (Han et al. 2015b); the default value is set to be 0, which means no BR is considered. Empirically, and understandably, a positive tolerance tends to facilitate convergence of the fixed-point algorithm. d) The maximum number of fixed-point iterations, upon which the algorithm will be terminated regardless of convergence. e) The step size α of the fixed-point algorithm; see Eq. (4.2).
Note that α needs to be adjusted for different networks and demand levels, as it has a major impact on convergence with respect to the criterion (5.1). Indeed, a very small step size α tends to terminate the fixed-point iterations prematurely while compromising the quality of the approximate DUE solution.
On the other hand, a very large α causes undesirable oscillations among the iterates; and it normally takes more iterations to reach convergence. Another consideration, arising from Eq. (4.2), is that α p (t, h k ) should be numerically comparable to h k p (t). f) Set the file path for the pre-processed network data (e.g. '.../Braess8 pp.mat').
3.2. Run the script 'MAIN.m'. 3.3. The output file, named 'DUE out.mat', will be saved to the working directory. The output includes the following variables:  Figure 14 illustrates the DUE solution on the Chicago Sketch network obtained from the MATLAB solver. This is done by loading the DUE path flow vector into the DNL solver (script file) before running 'OPEN NETWORK GUI.m'.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.