Dynamic origin–destination demand flow estimation under
congested traffic conditions
Chung-Cheng Lu
a, Xuesong Zhou
b,⇑, Kuilin Zhang
c aDepartment of Transportation and Logistics Management, National Chiao Tung University, Hsinchu 30050, Taiwan
b
Department of Civil and Environmental Engineering, University of Utah, Salt Lake City, UT 84112, USA
c
Transportation Research and Analysis Computing Center, Energy Systems Division, Argonne National Laboratory, Argonne, IL 60439, USA
a r t i c l e
i n f o
Article history: Received 22 July 2012
Received in revised form 10 May 2013 Accepted 13 May 2013
Keywords:
OD demand estimation Path flow estimator Lagrangian relaxation
Newell’s simplified kinematic wave theory
a b s t r a c t
This paper presents a single-level nonlinear optimization model to estimate dynamic ori-gin–destination (OD) demand. The model is a path flow-based optimization model, which incorporates heterogeneous sources of traffic measurements and does not require explicit dynamic link-path incidences. The objective is to minimize (i) the deviation between observed and estimated traffic states and (ii) the deviation between aggregated path flows and target OD flows, subject to the dynamic user equilibrium (DUE) constraint represented by a gap-function-based reformulation. A Lagrangian relaxation-based algorithm which dualizes the difficult DUE constraint to the objective function is proposed to solve the model. This algorithm integrates a gradient-projection-based path flow adjustment method within a column generation-based framework. Additionally, a dynamic network loading (DNL) model, based on Newell’s simplified kinematic wave theory, is employed in the DUE assignment process to realistically capture congestion phenomena and shock wave propagation. This research also derives analytical gradient formulas for the changes in link flow and density due to the unit change of time-dependent path inflow in a general network under congestion conditions. Numerical experiments conducted on three different networks illustrate the effectiveness and shed some light on the properties of the proposed OD demand estimation method.
2013 Elsevier Ltd. All rights reserved.
1. Introduction
Time-dependent origin–destination (OD) demand matrices are fundamental inputs for dynamic traffic assignment (DTA) models to describe network flow evolution as a result of interactions of individual travelers. Moreover, many emerging intel-ligent traffic management applications call for reliable estimates of dynamic OD demand, in order to generate proactive, coordinated traffic information provision and flow control strategies based on reliable traffic state estimates. Transportation authorities and practitioners have long been concerned about the unavailability of high quality time-dependent OD demand estimates which limits the potential for DTA deployments to analyze and alleviate traffic congestion. In the past decades, a rich body of literature, to be presented as follows, has been devoted to the methods of estimating static or time-dependent OD demand tables. However, the development of theoretically sound and practically deployable approaches for time-depen-dent OD demand estimation, particularly under congested conditions, remains a critical and challenging problem that is attracting significant attention from transportation researchers.
0968-090X/$ - see front matter 2013 Elsevier Ltd. All rights reserved.
http://dx.doi.org/10.1016/j.trc.2013.05.006
⇑ Corresponding author. Tel.: +1 801 585 6590; fax: +1 801 585 5477.
E-mail addresses:[email protected](C.-C. Lu),[email protected](X. Zhou),[email protected](K. Zhang). Contents lists available atSciVerse ScienceDirect
Transportation Research Part C
j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / t r c1.1. Literature review
To capture congestion effects in traffic networks, many researchers have attempted to integrate equilibrium assignment into the static OD demand estimation process.Nguyen (1977) and Leblanc and Farhangian (1982)incorporated link count observations into a variable demand user equilibrium (UE) assignment program as equality side-constraints so that the esti-mated link flows can reproduce observed link counts.Fisk (1989)combined the maximum entropy model with an UE assign-ment program to construct a bi-level mathematical programming problem.Yang et al. (1992) and Florian and Chen (1995)
further presented a more flexible bi-level framework to estimate consistent OD demand, where the upper level is a general-ized least squares (GLS)-based OD estimation model and the lower level is an UE assignment program.
Extending the concepts and solution methodologies of the static OD estimation problem,Cascetta et al. (1993)proposed a GLS estimator for dynamic OD demand in a general network. A simplified assignment model was used in their study; that is, path choice fractions are first calculated using a route choice model and then the resulting path flows are propagated to link flows based on link travel times.Tavana (2001)proposed a bi-level GLS optimization model and an iterative solution frame-work to estimate dynamic OD demand, while seeking to maintain internal consistency between the upper-level demand estimation problem and the lower-level DTA problem. Along this line, Zhou et al. (2003) and Zhou and Mahmassani (2006)extended this bi-level dynamic OD estimation approach to utilize multi-day traffic counts and automated vehicle identification (AVI) data, respectively.Van der Zijpp (1997)also addressed the dynamic OD demand estimation problem using data from induction loops and automated vehicle identification (AVI) equipment in the network. Based on a least-square modeling approach,Bierlaire and Crittin (2004)proposed an algorithm for sparse least squares that is computation-ally efficient for real-time estimation and prediction of dynamic OD demand tables.Zhou and Mahmassani (2007)developed a structural state space model for real-time traffic origin–destination demand estimation and prediction in a day-to-day learning framework.
In light of the above review, a typical bi-level dynamic OD demand estimation model needs to iteratively solve two opti-mization subproblems, namely upper-level and lower-level problems. The upper-level problem is the constrained ordinary/ generalized least squares (OLS/GLS) problem, with time-dependent OD flows as decision variables. This approach aims to minimize the following two deviation functions: (i) the deviation between observed and estimated link flows over all time intervals and (ii) the deviation between the target or historical demand and estimated demand matrices. The lower-level problem is the UE DTA problem, which determines a time-dependent network flow pattern that satisfies dynamic user equi-librium (DUE) conditions. However, it has been widely recognized that, under congested conditions, the mapping between demand inflow from the origin and link measurements is not a linear relationship as observed in the static case.
Yang (1995)provided two heuristic solution approaches for solving the general bi-level OD estimation problem, the iter-ative estimation-assignment (IEA) algorithms and sensitivity-analysis based algorithms (SAB).Tavana (2001)suggested that the IEA algorithm still provides solutions to the Cournot–Nash game, rather than the Stackelberg game in a bi-level program, because the upper-level optimization model in IEA does not consider the dependence of link-flow proportions on the OD flows. Based on the SAB approach, an alternative nonlinear least squares formulation was proposed byTavana (2001)to explicitly consider the changes of link flow proportions due to the adjustments in dynamic OD flows. Additionally, numerical derivatives of link flow proportions with respect to OD flows are obtained from a mesoscopic DTA simulation program ( Jaya-krishnan et al., 1994; Mahmassani, 2001). On the other hand, a standard SAB algorithm needs to approximate the derivatives through simulation for each OD pair and each time interval at every iteration, which is computationally intensive, especially for large-scale networks. Recently,Balakrishna et al. (2008) and Cipriani et al. (2011)introduced gradient approximation methods within a simultaneous perturbation stochastic approximation (SPSA) framework in order to reduce the number of simulation runs when calculating numerical derivatives or gradients. Artificial intelligence techniques have recently been adopted to the context of dynamic OD demand estimation to replace iterative procedures. For instance,Kattan and Abdulhai (2006)proposed a non-iterative approach to dynamic OD demand estimation based on a machine-learning technique using advanced parallel evolutionary algorithms. More recently,Huang et al. (2012)developed an approach to estimate travel de-mand in large-scale microscopic traffic simulation models based on a Guided Genetic Algorithm with a distributed imple-mentation to improve computational efficiency and reduce memory requirements.
Intending to develop an internally consistent approach for the dynamic OD demand estimation problem, single-level path flow estimators (PFEs) have been proposed for the static OD estimation problem; e.g., the linear programming PFE bySherali and Park (2001)on estimating UE path flows, and the nonlinear programming PFE byBell et al. (1997)on estimating stochas-tic UE path flows. Recently,Nie and Zhang (2008)formulated a novel single-level formulation based on variational inequal-ities (VI), which utilizes the dynamic link-path incidence relationships in a generic projection-based VI solution framework. Based on a single level reformulation,Lundgren and Peterson (2008)proposed a descent heuristic, which is an adaptation of the projected gradient method, to the dynamic OD demand estimation problem. The search direction was determined by approximating the Jacobian matrix representing the derivatives of the link flows with respect to a change in the OD-flows. By adapting the analytical approach ofGhali and Smith (1995)for evaluating the local link marginal travel times,Qian and Zhang (2011)further incorporated the travel time gradients into the single-level OD estimation framework proposed byNie and Zhang (2008), in order to utilize travel time measurements, while linear mappings between OD flow and link flows are still assumed in their upper-level OD estimation problem. In addition,Nie et al. (2005), Nie and Zhang (2010), and Shen and Wynter (2011)integrated the integral term of UE formulation (Beckmann et al., 1956) with the measurement deviations (in
terms of the GLS objective function) to develop alternative single-level path flow formulations for static OD estimation, and their models can be viewed as a special case of UE assignment with elastic demand.
To fully capture the nonlinear dependency between dynamic OD demand and heterogeneous traffic measurements, such as link flow, density and travel time, in a general network, it is necessary to incorporate a reliable dynamic network loading (DNL) component into the OD demand flow estimation process. The existing DTA methodologies are classified into two ma-jor groups: analytic approach and simulation-based approach. The former includes three types of formulations that have the potential for deriving theoretical insights: mathematical programming, optimal control and VI. Earlier studies (e.g.,Merchant and Nemhauser, 1978; Carey, 1987; Friesz et al., 1989) used link/node exit flow constraints to propagate traffic flows and link performance functions to determine path travel costs. While using well-defined exit constraints and cost functions makes it theoretically possible to establish the mathematical properties of solutions (e.g., existence and uniqueness), these analytical models suffer from several limitations in terms of representing the dynamics and complexity of real-world traffic flow sys-tems. For example, the models are still subject to difficulty in ensuring the first-in–first-out (FIFO) queueing discipline and capturing spillback queues.
In the pioneering works of capturing shock waves and congestion (i.e., queue formation, spillback and dissipation),
Lighthill and Whitham (1955) and Richards (1956)proposed the kinematic wave (KW) theory, which rigorously describes traffic flow dynamics by integrating flow conservation constraints and traffic flow models to a set of partial differential equa-tions (PDE). Because it is difficult to obtain analytical soluequa-tions for these PDEs, many researchers have presented various finite-difference approximations to solve the equations numerically. Based on a triangular flow-density relation, Newell pro-posed a simplified KW model (Newell, 1993) to better keep track of shock waves and queue propagation using cumulative flow counts on links. By discretizing a link into a set of homogeneous unit cells,Daganzo (1994, 1995)presented a cell trans-mission model that adopts a ‘‘supply–demand’’ or ‘‘sending–receiving’’ framework to model flow dynamics between two successive discretized cells in the link.
1.2. Overview and contribution of proposed method
This paper contributes to the growing body of literature on dynamic OD demand estimation as follows:
Instead of working on the commonly-used OD flow variables, this study presents a new path flow-based optimization model for jointly solving the complex OD demand estimation and UE DTA problems. This model simultaneously mini-mizes (i) the deviation between measured and estimated traffic states and (ii) the deviation between aggregated path flows and target OD flows, subject to a dynamic user equilibrium (DUE) constraint, which is reformulated using an equiv-alent gap function. Working in this path-flow dimension, our formulation can directly aggregate estimated path flows to obtain final OD flow patterns, and obviate explicit dynamic link-path incidences, as opposed to the majority of previous studies.
By dualizing the difficult DUE constraint into the objective function, this research proposes an effective Lagrangian relax-ation-based solution framework. The relaxed problem can be viewed as a simultaneous route and departure time user equilibrium (SRDUE) problem with elastic demand The final solution is a set of path flows satisfying ‘‘tolled user equilib-rium’’ (Lawphongpanich and Hearn, 2004), where the deviation with respect to traffic measurements can be viewed as an additional penalty for over-estimated or under-estimated path flows. By incorporating heterogeneous real-world mea-surements in the objective function, such as link densities from video surveillance and road side detectors, the proposed estimation model fully utilizes available information to reflect route choices in a congestion network.
A DNL model that encapsulates Newell’s simplified KW model in a mesoscopic traffic flow simulation framework is pro-posed to describe congestion phenomena, such as queues formation, spillback, and dissipation. Explicitly using the cumu-lative arrival and departure curves, Newell’s traffic flow model provides a rigorous mathematical formulation to realistically represent traffic dynamics and capture the impact of shock waves on various macroscopic traffic measures. Based on the proposed DNL model, this research derives analytical, local gradients of different measurement types, such as link flow and density, with respect to path inflows. This valuable gradient information not only considers the depen-dences of link flow/density changes on the OD/path flow, but also allows for computing feasible descent directions in an efficient gradient-projection-based method embedded in the Lagrangian relaxation-based solution framework.
While the general relaxation scheme and some common features (e.g., path flow variables) are shared by our work and some previous research (e.g.,Nie and Zhang, 2008), our approach differs significantly from the existing methods in the fol-lowing aspects.
The proposed dynamic OD demand estimation model is a GLS-based bi-level formulation that adopts a gap function to measure the deviation from the DUE conditions, while previous research used complementary equations or VIs to describe optimality or equilibrium conditions of the problem.
The proposed Lagrangian relaxation-based algorithm dualizes the gap function-based DUE constraint into the objective function, and solves the single-level relaxation problem by reducing the difference between the upper and the lower bounds. On the other hand, previous research developed different (e.g., column generation) algorithms to solve the VI-based single-level model.
Based on the proposed DNL model, this research derives analytical, local gradients of different measurement types with respect to path inflows. These gradients are different from the link and path marginal travel time or cost typically used in system optimal DTA problems. The proposed algorithm explicitly utilizes the gradient information to compute feasible descent directions in an efficient gradient-projection-based method embedded in the Lagrangian relaxation-based solu-tion framework. Although previous research had applied basic projecsolu-tion algorithm (BPA) and extra-gradient method (EGM) to the dynamic OD demand estimation problem, the issue of gradient derivations has not been explicitly addressed.
This paper is organized as follows. The single-level time-dependent path flow estimation models and the mesoscopic DNL model based on Newell’s simplified KW model are delineated in Section2. Section3presents the Lagrangian relaxation-based solution algorithm. Section 4describes the derivation of the gradients of typical traffic measurements. Numerical studies on a simple corridor and a real-world network with time-dependent sensor data are presented in Section5. Section6
concludes this study.
2. Single-level time-dependent path flow estimation models 2.1. Notation and problem statement
Notation and variables for the proposed single-level dynamic OD demand estimation framework are listed as follows. Set:
A: set of links W: set of OD pairs P: set of paths
S: set of links with sensors, S # A
Hd: set of discretized departure time intervals Ho: set of discretized observation time intervals
Index:
t: index of simulation time intervals, t = 0, . . . , T. This paper refers to any particular time interval t as the time t.
s
: index of departure time intervals,s
2 Hd w: index of OD pairs, w 2 Wp: index of paths for each OD pair, p 2 P l: index of links, l 2 A
Input parameters of link attributes length(l): length of link l
nlanes(l): number of lanes of link l kjam(l): jam density on link l
nmax(l): number of vehicles that can be stored on link l, nmax(l) = k
jam(l) nlanes(l) length(l)
v
f(l): free-flow speed on link lwb(l): backward wave speed on link l
FFTT(l): free-flow travel time on link l, FFTT(l) = length(l)/
v
f(l)BWTT(l): backward wave travel time on link l, BWTT(l) = length(l)/wb(l) capout(l, t): maximal outflow capacity of link l at time t
capin(l, t): maximal inflow capacity of link l at time t
Traffic measurements inputs
qðl; tÞ: observed number of vehicles passing through an upstream detector on link l during observation interval t
kðl; tÞ: observed density on link l during observation interval t
dðwÞ: target demand, which is the total traffic demand for OD pair w over a planning horizon Estimation variables
r(w,
s
, p): estimated path flow on path p of OD pair w and departure time intervals
c(w,
s
, p): estimated path travel time on path p of OD pair w and departure time intervals
p
(w,s
): estimated least path travel time of OD pair w and departure time intervals
q(l, t): estimated number of vehicles passing through an upstream detector on link l during observation interval t k(l, t): estimated density on link l during observation interval t
d(w,
s
): estimated demand of OD pair w and departure time intervals
Given a road network G with a set of nodes and a set of links, the input data to the dynamic OD demand estimation prob-lem include a set of time-dependent observation data (e.g. flows and densities) on a subset of links with sensors for a set of discretized observation time intervals, Ho, and target OD demand matrices for a set of departure time intervals, Hd. The prob-lem aims to find a vector of time-dependent path flows r = {r(w,
s
, p), "w,s
, p} such that the estimated flow pattern can both match the time-dependent observation data and satisfy the DUE conditions (i.e., for each OD pair w and each departure time intervals
, the travel time of a path with a positive path flow is the same and equal top
(w,s
), while the travel time of a path with zero flows is larger thanp
(w,s
)). Then, the time-dependent OD demand estimates can be obtained by aggregating path flows for each OD pair and each departure time interval.2.2. A path flow-based nonlinear program
Given sensor data (i.e. observed link flows and densities) and target (aggregated historical) OD demands, the proposed single-level time-dependent path flow estimation model is a nonlinear program with the path flows r(w,
s
, p), "w,s
, p and least path travel timesp
= {p
(w,s
), "w,s
} as the decision variables. Denote c = {c(w,s
, p), "w,s
, p}, q = {q(l, t), "l, t} and k = {k(l, t), "l, t}. The nonlinear program is presented as follows.P1: Nonlinear program Min Z ¼ bd X w X s2Hd X p rðw;
s
;pÞ dðwÞ " #2 þX l2S X t2Ho bq½qðl; tÞ qðl; tÞ 2 þ bk kðl; tÞ kðl; tÞ 2 n o ; ð1Þ Subject to ðc; q; kÞ ¼ DNLFðrÞ; ð2Þ gðr;p
Þ ¼R
wR
sR
pfrðw;s
;pÞ½cðw;s
;pÞp
ðw;s
Þg ¼ 0; ð3Þ cðw;s
; pÞp
ðw;s
Þ P 0;8
w;s
; p; ð4Þp
ðw;s
Þ P 0;8
p 2 Pðw;s
Þg8
wg;8
w;s
;s
ð5Þ rðw;s
;pÞ P 0;8
w;s
; p: ð6ÞThe objective function, Eq.(1), minimizes the weighted sum of the deviation between estimated time-dependent OD de-mands (or aggregated path flows) and target dede-mands and the deviation between estimated and observed link flows and densities, where bd, bqand bkare the weights reflecting different degrees of confidence on target OD demands and observed link flows and densities, respectively. In a GLS framework, those weights can be viewed as the inverses of the variances of the distinct sources of measurements. In the objective function, OD demand d(w,
s
) is substituted by its corresponding aggre-gated time-dependent path flows, so the OD demand and path flow balance constraints are not included in the model, P1. Eq.(2)states that estimated path travel times (c), and link flows (q) and densities (k) are obtained by a given DNL function of path flows, DNLF(r). This research proposes the DNL model based on Newell’s simplified KW model (Newell, 1993). The details of the DNL model are presented in Sections2.3 and 2.4of the current paper.The gap function-based constraint in Eq.(3)is adopted to satisfy the DUE conditions in the proposed model, whereas other existing single-level OD demand estimation models incorporated VI (e.g.,Nie and Zhang, 2008), or Beckmann’s integral term (e.g.,Shen and Wynter, 2011) to represent the DUE conditions. The non-negative gap function g(r,
p
) is used to measure the deviation from the DUE conditions, and it vanishes when the DUE conditions are satisfied (e.g.,Lu et al., 2009). Eqs.(4) and (5)are definitional constraints of the least travel timep
(w,s
) for each (w,s
). Eq.(6)are non-negative constraints for path flow variables.2.3. Modeling queue spillback through Newell’s simplified kinematic wave model
The proposed time-dependent path flow estimation models, P1, encapsulats the DNL model based on Newell’s simplified KW model (Newell, 1993) to describe traffic congestion propagation (e.g., queue build-up, spillback, and dissipation) in a road traffic network. As shown inFig. 1, Newell’s model is concerned about three state variables on each link l: (i) cumulative flow count A(l, t) for vehicles moving into link l through the upstream node, (ii) cumulative flow count V(l, t) for vehicles wait-ing at the vertical queue of the downstream node of link l at time t, and (iii) cumulative flow count D(l, t) for vehicles movwait-ing out of link l through the downstream node.
Let x be the location along the corridor, and N(x, t) the cumulative flow count at location x and time t of a link. The change of N(x, t) along a characteristic line (wave) is represented as follows.
dNðx; tÞ ¼@N @xdx þ
@N
A wave represents the propagation of a change in flow and density along the roadway, and the wave speed is the slope of the characteristics line w ¼@q
@k¼ dx
dt. Along the movement of a wave, we substitute dt ¼ dx
winto Eq.(7), so that we can link the difference of cumulative flow counts together through
dNðx; tÞ ¼ qdt kdx ¼ k þq w
dx: ð8Þ
For the triangular shaped flow-density relation with constant forward and backward wave speeds, it is easy to verify that, when the speed of the forward wave is
v
f, the general cumulative flow count updating formula, Eq. (7), reduces to k þvqf¼ k þ k ¼ 0. Under congested traffic conditions with a constant backward wave speed wb, we have
k þ q
wb¼ kjam, and this equation can be rewritten as
dN ¼ k þ q wb
dx ¼ kjamðlÞ lengthðlÞ nlanesðlÞ: ð9Þ
Eq.(9)is used to describe how a backward wave travels through the link. As shown inFig. 1, when a queue spills back from the downstream to the upstream, the arrival and departure cumulative flow counts at two ends of a link (at timestamps t and time t-BWTT(l)) need to ensure a constant difference of dN = kjam(l) length(l) nlanes(l), and the capacity restriction is propagated throughout the link using a time duration of BWTT(l) = length(l)/wb(l).
The space–time plot inFig. 2further demonstrates how Newell’s model uses cumulative counts to model the forward and backward waves in describing queue phenomenon on two links a and b. Exactly at time t 1, the tail of the queue at the downstream link b propagates to its upstream node. That is, if the cumulative outflow count at a lagged time stamp t-BWTT(b) 1 on link b equals to the cumulative inflow count at time t 1, then a queue spillback occurs as a backward wave is able to propagate through the congested time–space ‘‘mass’’ of the link. Compared to the conventional vertical queue mod-el, where the inflow rate of a link is not constrained by the available physical length, the inflow rate in Newell’s model is governed by the discharge flow of a link through the backward wave propagation process. Furthermore the number of vehi-cles on a lane should not be greater than kjam length(b). As a result, the proposed model is able to detect and capture queue spillbacks to upstream link(s).
2.4. Traffic flow dynamics constraints and states updating
The main building block of the DNL model is a set of demand, path and link flow balance constraints. Firstly, the demand flow balance constraint for each (w,
s
) pair over its path set is as follows.X
p2Pðw;sÞ and Uðt;sÞ¼1
rðw; t; pÞ ¼ qðw;
s
Þ;8
w;s
; ð10ÞwhereU(t,
s
) = 1 indicates that time interval t is inside departure time intervals
. Then, the DNL needs to load path flow r(w, t, p) to the first link of path p 2 P(w,s
).Aðp; lðp; 1Þ; tÞ ¼ rðw; t; pÞ;
8
p 2 Pðw;s
Þ; w; t: ð11Þwhere l(p, 1) represents the first link of path p, and A(p, l, t) is the cumulative number of vehicles that are assigned to path p and have arrived at the upstream end of link l in time t. Cumulative arrival counts are moved from the beginning of link l to its vertical queue for each p, l, and t. That is,
Vðp; l; tÞ ¼ Aðp; l; t FFTTðlÞÞ;
8
p; l; t; ð12Þwhere V(p, l, t) is the cumulative number of vehicles that are assigned to path p and have waited at the vertical queue of link l at time t.
The cumulative departure flow is less then the cumulative number of vehicles in the vertical queue on a link.
Dðp; l; tÞ 6 Vðp; l; tÞ;
8
p; l; t; ð13Þwhere D(p, l, t) is cumulative number of vehicles that are assigned to path p and have departed from link l at time t. The cumulative departure flow count of link l is transferred to the cumulative arrival count on the next link l + 1 along the path p,
Aðp; l þ 1; tÞ ¼ Dðp; l; tÞ;
8
p; l; t: ð14ÞThe balance constraint at the super sink for each OD pair is
Dðp; L; TÞ ¼ rðw; t; pÞ;
8
w; p; t: ð15Þwhere L is the last link of path p and T is the end of the planning horizon. By summing up all cumulative arrival, vertical queue and departure counts across different paths, we can obtain cumulative link arrival, vertical queue and departure counts as follows. Aðl; tÞ ¼X p Aðp; l; tÞ;
8
l; t: ð16Þ Vðl; tÞ ¼X p Vðp; l; tÞ;8
l; t: ð17Þ Dðl; tÞ ¼X p Dðp; l; tÞ;8
l; t: ð18ÞIn addition to the above set of flow balance constraints, when vehicles are discharged from the queue, Eq.(19)is used to satisfy the outflow capacity constraints and capture the effect of queue spillback.
Dðl; tÞ Dðl; t 1Þ ¼ min Vðl; t 1Þ Dðl; t 1Þ; q maxðl; tÞ; capoutðl; tÞ ;
8
l; t: ð19Þwhere V(l, t 1) D(l, t 1) is the number of vehicles waiting at the vertical queue of link l at time t. To determine the max-imum outflow capacity capout(l, t) in Eq.(19), we consider a corridor with two consecutive links l and l + 1, where a lane drop bottleneck at the downstream end of link l + 1 causes congestion and the queue spills back to link l (i.e., link l + 1 is fully
congested while link l is partially congested). Under congested conditions, capout(l, t) is determined by the bottleneck inflow rate of link l + 1.
capoutðl; tÞ ¼ capinðl þ 1; tÞ: ð20Þ
As shown inFig. 3, the inflow capacity capin(l + 1, t) is defined in terms of the difference of the cumulative arrival flow counts between two consecutive time stamps t 1 and t as follows:
capin
ðl þ 1; tÞ ¼ Amaxðl þ 1; tÞ Aðl þ 1; t 1Þ; ð21Þ
where (according to the backward wave propagation constraint in Eq.(7))
Amaxðl þ 1; tÞ ¼ Dðl þ 1; t BWTTðl þ 1ÞÞ þ kjamðl þ 1Þ lengthðl þ 1Þ nlanesðl þ 1Þ: ð22Þ
Finally, the queue spillback constraint based on Newell’s simplified kinematic wave mode is
Aðl; tÞ 6 Dðl; t BWTTðlÞÞ þ kjamðlÞ lengthðlÞ nlanesðlÞ;
8
l; t: ð23ÞFig. 3shows that if the cumulative outflow count at a lagged time stamp t BWTT(l + 1) 1 on link l + 1 equals to the cumulative inflow count at time t 1, then a queue spillback occurs as a backward wave propagates through the congested link.
The numerical computation scheme of the above DNL model can be implemented as an event-based simulation process which does not keep track of vehicles’ positions on links at each time stamp. Specifically, if a vehicle j is moved into link l at its arrival time Arr(j, l), the time entering the vertical queue at the stop bar is tvq= Arr(j, l) + FFTT(l). When the simulation time clock advances to tvq, if the out-flow link capacity at time tvqis still available for vehicle j to exit from the vertical queue, then this vehicle moves to the next link of its path and its departure time stamp Dep(j, l) = tvq; otherwise vehicle j has to wait in the vertical queue for the available outflow capacity. As this proposed model strictly satisfies the first-in-first-out (FIFO) constraint, if a vehicle j in the beginning of the vertical queue is blocked due to unavailable outflow capacity, then vehicles behind j in the queue will be also blocked as well. To model complex geometric features, such as short left-turn bays on a multi-lane facility, one needs to decompose a link into multiple connected cells with each cell satisfying the FIFO constraint.
To describe traffic flow dynamics with multiple OD pairs and paths at merge and diverge junctions, many existing DNL models (e.g.?Jayakrishnan et al., 1994; Daganzo, 1995; Tampère et al., 2011) proposed various modeling approaches to determine outflow capacities for merging vehicular flows moving out from incoming links to a link, or outflow proportions for diverging vehicular flow moving out from a link to its outgoing links. In this paper, the available outflow capacity of an incoming link at a merge is assumed to be proportional to the number of lanes on that link. Similar to a mesoscopic DNL model(Jayakrishnan et al., 1994), each vehicle carriers its own OD and path information in our proposed model, so the out-flow proportions from a link are indirectly determined by the path attributes associated with vehicles waiting at that link’s vertical queue, which is implemented as a FIFO queue that explicitly enables the FIFO constraint to be satisfied.
3. Solution algorithm
This section describes the Lagrangian relaxation-based heuristic for solving the single-level time-dependent path flow estimation model, presented in Section2. We propose the following heuristic solution method to efficiently obtain good solutions for problem instances on road networks of practical sizes. The heuristic integrates Lagrangian relaxation and col-umn generation methods to solve the time-dependent path flow estimation model, P1. The gap function constraint Eq.(3)is relaxed to the objective function Eq.(1)with a non-negative Lagrange multiplier k. The resulting Lagrangian subproblem is given as follows:
P2 : Min
r;p Lðr;
p
;kÞ ¼ Z þ kfgðr;p
Þg ð24ÞSubject to constraints(2), (4), (5) and (6), where r and
p
are the vectors of path flows and least path times respectively. For a given k, the solution to P2 provides a lower bound to P1. The Lagrangian dual problem is given as follows.P3 : Max
k Minr;p Lðr;
p
;kÞ ð25ÞSubject to k P 0.
The heuristic consists of two major algorithmic steps: at each iteration n, (i) given a Lagrange multiplier k(n), find an opti-mal path assignment r(n)and least path travel times
p
(n)by solving the Lagrangian subproblem, P2 and (ii) given a vehicle path assignment r(n)and least path travel timesp
(n), update the Lagrange multiplier k(n+1)by using the following rule.kðnþ1Þ¼ max 0; kðnÞþ
a
ðnÞR
w
R
sR
prðw;s
;pÞ½cðw;s
;pÞp
ðw;s
Þ
n o
; ð26Þ
where
a
(n)is the step size for updating the Lagrange multiplier.Accordingly, this heuristic has two loops (seeFig. 4). The outer loop is for updating the Lagrange multiplier using the rule described in Eq.(26). For each outer loop iteration n (i.e., corresponding to a given Lagrange multiplier k(n)), a column gen-eration-based approach is used to solve the Lagrangian subproblem P2. This approach forms an inner loop for solving a DUE assignment problem under a restricted feasible solution space. In each inner loop iteration m, a time-dependent shortest path algorithm (Ziliaskopoulos and Mahmassani, 1993) is adopted to generate time-dependent least time paths and to aug-ment the restricted path set. In light of the time-dependent shortest path algorithm, the least path travel times
p
(m)are ob-tained to satisfy the constraints Eqs.(4) and (5), thus, these definitional constraints of the least travel times can be droppedin solving the restricted Lagrangian subproblem. To solve the restricted subproblem, a gradient-projection-based descent direction method (Lu et al., 2009) is used to update path flows r(m+1), while maintaining the feasibility of non-negativity con-straints Eq.(6). Specifically,
rðw;
s
;pÞmþ1¼ Max 0; rðw;s
;pÞmc
ðmÞ b dr
h d ðrÞjr¼rðmÞþ bqr
h q ðrÞjr¼rðmÞþ bkr
h k ðrÞjr¼rðmÞþ kðnÞr
gðr;p
Þjr¼rðmÞ h i n o ð27Þwhere
c
(m)is the step size, and the gradients, which consist of the first-order partial derivatives with respect to a path flow variable r(w,s
, p), can be derived as follows.r
hdðrÞ ¼ @ Ps2H d;p2Pðrðw;s;pÞÞ dðwÞ n o @rðw;s
;pÞ ¼ 2 X s2Hd X p2P rðw;s
;pÞ dðwÞ ! ð28Þr
hqðrÞ ¼ @Pl2S;t2H ot qðl;tÞðrÞ qðl;tÞ h i @rðw;s
;pÞ ¼ 2 X t2Ho X l2S ½qðl; tÞ qðl; tÞ @qðl; tÞðrÞ @rðw;s
;pÞ ð29Þr
hkðrÞ ¼@ P l2S;t2Hot kðl;tÞðrÞ kðl;tÞ @rðw;s
;pÞ ¼ 2 X t2Ho X l2S kðl; tÞ kðl; tÞ @kðl; tÞðrÞ @rðw;s
;pÞ ð30Þr
gðr;p
Þ ¼ @gðr;p
Þ @rðw;s
;pÞ¼ cðw;s
;pÞp
ðw;sÞþ rðw;s
;pÞ @cðw;s
;pÞ @rðw;s
;pÞ ð31ÞEstimated link flows, densities, and link/path travel times and the corresponding partial derivatives, namelyrhq(r),r hk(r) and @c(w,
s
, p)/@r(w,s
, p) are obtained from the DNL model presented in Section2.The steps of this algorithm are presented as follows. Algorithm 1. Lagrangian relaxation-based heuristic
Step 1: Initialization
Step 1.1: Set outer loop iteration counter n = 0. Input a historical demand table, dðwÞ; 8w, its corresponding tem-poral distribution profile, and time-dependent link observations (densities, and flows). Initialize k(n)to a positive value, such as 1.0.
Step 1.2: Perform DUE traffic assignment with time-dependent OD demand matrix, d(n), based on dðwÞ; 8w and a temporal profile, and build a feasible path set bP.
Step 1.3: According to the assignment results (i.e., estimated link densities, flows, and travel times), compute the upper bound of the optimal solution to the primal problem P1 as follows:
ZUB ¼X l2S X t2Ho fbq½qðl; tÞ qðl; tÞ 2 þ bk½kðl; tÞ kðl; tÞ 2 g þ knfgðr;
p
Þg: ð32Þ Note thatPw P s2Hd P prðw;s
;pÞ dðwÞ h i2¼ 0, as the DUE assignment loads the target demands. Step 2: [Inner loop] Solve the Lagrangian subproblem, P2, to find the optimal path flows corresponding to k(n).
Step 2.1: Set inner loop iteration counter m = 0. Read estimated link flows, densities, path flows, path travel times and corresponding gradientsrhd(r),rhq(r) andrhk(r) from the last outer iteration n.
Step 2.2: Identify the least time paths and determine
p
(m) . Step 2.3: Calculate the gradients according to Eqs.(28)–(31).Step 2.4: Determine the step size
c
(m)according to the method of successive averages (MSA),c
(m)= 1/(m + 1). Step 2.5: Update path flows r(m+1), according to Eq.(27).Step 2.6: Load path flow assignment r(m+1)to the DNL model based on Newell’s simplified KW model, to obtain estimated link densities and flows, and path travel times.
Step 2.7: Convergence checking for solving the Lagrangian subproblem, P2 (i.e., the inner loop). Update the objec-tive function, L(r(m+1),
p
(m+1), k(n)). If m < Mmax(maximum number of inner loop iterations) or L(r(m+1),
p
(m+1), k(n)) is improved, then m = m + 1, and go to Step 2.2; otherwise, go to Step 2.8.Step 2.8: Update the tightest lower bound, ZLB, as the following:
ZLB¼ Max ZLB; bd h d rðmÞ þ bqh q rðmÞ þ bkh k rðmÞ þ kðnÞg r ðmÞ;
p
ðmÞ n o : ð33ÞStep 3: Update the upper bound.
Step 3.1: Construct a time-dependent OD demand matrix by aggregating the time-dependent path flows obtained from the solution to the Lagrangian subproblem P2, for each OD pair and each departure time interval.
dðw;
s
Þ ¼R
p2Prðw;s
; pÞ;8
w;s
: ð34ÞStep 3.2: Perform DUE traffic assignment for the constructed OD demand matrix d = {d(w,
s
), "w,s
}. According to the assignment results (i.e., estimated link densities, flows, and travel times), obtain the equilibrium path flows and update the tightest upper bound as:ZUB¼ Min ZUB;bdh d ðrðmÞÞ þ b qh q ðrðnÞÞ þ b kh k ðrðnÞÞ þ kðnÞ½gðrðnÞ;
p
ðnÞÞ n o : ð35ÞStep 4: Convergence checking for solving the Lagrangian dual problem (i.e., the outer loop). If n > Nmax(maximum num-ber of outer loop iterations) or ZUB ZLB< / (a preset threshold), then stop; otherwise, set n = n + 1.
Step 5: Path generation. Compute least time path for each departure time and OD pair using the time-dependent short-est path algorithm, developed byZiliaskopoulos and Mahmassani (1993), and add newly generated paths to the feasible path set bP.
Step 6: Update Lagrange multiplier k.
Step 6.1: Obtain estimated link flows, path flows, and path travel times from the final iteration of the last inner loop (i.e., solving P2 in Step 2). Update the gap value, g(r(n),
p
(n)) defined in Eq.(3), using these estimates. Step 6.2: Determine the step size,a
(n), by MSA, i.e.a
(n)= 1/(n + 1).Step 6.3: Update Lagrange multiplier, k(n+1), according to Eq.(26). Return to Step 2 (the inner loop).
Typically, a Lagrangian solution framework requires obtaining exact solutions to relaxed subproblems. It should be re-marked that, analyzing existence and uniqueness of solutions to the DUE problem for multiple OD pairs are very challenging, and the gradient-based algorithm through Eqs.(27)–(31)cannot guarantee that the relaxed (nonlinear) problem P2 is solved to its optimality. Thus, when no global optimum solution is available for P2, the proposed overall Lagrangian solution algo-rithm is still a heuristic method in nature.
4. Evaluation of partial derivatives with respect to path flow perturbation
Solving the proposed single-level dynamic OD estimation model requires the evaluation of the partial derivatives with respect to time-varying path flows, i.e.,@qðl;tÞðrÞ
@rðw;s;pÞ;
@kðl;tÞðrÞ
@rðw;s;pÞ, and
@cðw;s;pÞðrÞ
@rðw;s;pÞ in Eqs.(29)–(31). These partial derivatives represent the
marginal effects of an additional unit of path inflow on link flow and density and path travel time. This section delineates the evaluation of these partial derivatives due to path flow perturbation in a congested network, based on cumulative link inflow and outflow curves. The following notation is used throughout this section.
L: the number of links on the path, l: link index l = 1, 2, . . . , L,
t0
l: the time when an additional unit of perturbation flow arrives at link l, t00
l: the time when an additional unit of perturbation flow departs from link l, tqsl : the time when the queue starts to form on link l,
tB
l: the time when the queue vanishes on link l, tA
l :¼ tBl FFTTðlÞ,
tql : the time when the queue on link l starts to spillback to its upstream link l 1, nA: cumulative arrivals at time tA
l, n0: cumulative arrivals at time t0
l.
4.1. Evaluation of link partial derivatives on a congested link
In this study, the link partial derivatives are referred to as the changes in link flow and density due to an additional unit of link/path inflow. In their pioneering work,Ghali and Smith (1995)presented an analytical approach to evaluate (local) link marginal travel time (or delay) on a congested link, based on link cumulative flow curves. An illustration of the approach is depicted inFig. 5. In the figure, the two solid lines represent the cumulative arrival and departure curves, while the dashed line represents the cumulative vertical queue. The (outflow) capacity of the link is c. The key result of their approach is that the link marginal delay equals the grey area. While link and path marginal delays were investigated in the context of system optimal DTA (e.g.,Ghali and Smith, 1995; Peeta and Mahmassani, 1995; Shen et al., 2007; Qian and Zhang, 2011), to the authors’ current knowledge, the study on the partial derivatives of link flow and density and path travel time with respect to path flow was nonexistent.
(For example, the queue starts at tqs
l ¼ 7 : 00 AM, an additional unit of vehicle, n’=1000, enters link l at time t0l¼ 7 : 10 AM, leaves the link at time t00
l ¼ 7 : 50 AM, the queue dissipates on link l at tBl ¼ 8 : 30 AM, and tAl ¼ 8 : 15 AM, assuming FFTT (l) = 15 min.)
The following propositions can be directly induced fromFig. 5for deriving the marginal effects on link flow (inflow and outflow) and density.
Proposition 1. Under free-flow conditions, an extra unit of flow arriving at the upstream end of link l at time t0
l results in the
following: (i) the link inflow and outflow increase by 1 at times t0
land t00l, respectively, and the flow rates at other time intervals do
not change; (ii) the link density increases by 1 from t0
l to t00l; and (iii) the individual travel times are not changed, and
t00
l ¼ t0lþ FFTTðlÞ.
Proposition 2. Under partially congested conditions and constant link (outflow) capacity c, an extra unit of flow arriving at the upstream end of link l at time t0
lresults in the following: (i) the link inflow and outflow increase by 1 at times t0land tBl, respectively, and the flow rates at other time intervals do not change; (ii) the link density increases by 1 from t0
lto tBl; and (iii) the flows arriving between t0
land tAl experience the additional delay 1/c, because it takes 1/c to discharge this perturbation flow.
Proposition 3. Under fully congested condition, an extra unit of flow arriving at the upstream end of link l at time t0
lwould not change the link inflow, outflow, and capacity. Both link inflow and outflow remain the same at the maximum link flow rate, while the link density is equal to the jam density. The flows arriving between t0
land tAl experience the additional delay 1/c, because it takes 1/c to discharge this perturbation flow.
With the proposed DNL model based on Newell’s simplified KW model, we can adapt Eq.(9)to detect if the queue spills back to the upstream end of link l using the following equation.
Aðl; tÞ < Dðl; t BWTTðlÞÞ þ kjamðlÞ lengthðlÞ nlanesðlÞ: ð36Þ
Specifically, if the above strict inequality holds, then the queue has not propagated back to the upstream end of the link (or the link is partially congested). Otherwise, the link is fully congested and the queue propagates throughout the link along the backward wave line, as illustrated inFig. 5. Note that, under fully congested conditions, the prevailing link density k(l) at time t might be smaller than kjam(l).
A common pitfall for deriving the partial derivative of link density under congested conditions, is to record the increase in density by 1 from t0
lto t00l (e.g., from 7:10 AM to 7:50 AM with the duration of 40-min inFig. 5).Proposition 2, induced from
Fig. 5, clarifies that the actual change in link density would last until the queue vanishes at time tB
l, so the change in link density should cover a duration of 80-min from 7:10 AM to 8:30 AM in our example.Proposition 2also indicates that the change in link outflow, due to an extra unit of flow arriving under congested conditions, occurs at the time tB
l (e.g. 8:30 AM), rather than t00
l (e.g., 7:50 AM). Note that, by multiplying the (individual) extra delay, 1/c, by the number of vehicles arriving between t0
land tAl (i.e., n A
n0), we can obtain that, if t0
lis between t qs
l and tAl, the local link marginal delay is equal to tB
l t0l, which is the sum of the traversal time of the perturbation flow, t00l t0l, and the additional delay imposed by the per-turbation on others is equal to tB
l t00l. This evaluation approach for link marginal travel times was also adopted by,Shen et al.
(2007) and Qian and Zhang (2011)in the context of system optimal dynamic traffic assignment. 4.2. Evaluation of the impact of path flow perturbation on two sequential links
In static transportation networks, the path travel time marginal, which is the partial derivative of the total system-wide travel time, can be obtained as the sum of the link marginals of a path’s constituent links. However, in dynamic and con-gested transportation networks, this additivity assumption may lead to a severe deficiency in the evaluation of path margin-als, as indicated byShen et al. (2007). They showed that it is necessary to explicitly trace the propagation of path flow perturbation in evaluating path marginal travel times and proposed an evaluation method of path marginals. Thus, this study evaluates the impact of path flow perturbation (i.e., the partial derivatives) in the individual link-time level by tracing the changes in link flow and density on a sequence of links (i.e., a path) and over different time intervals, due to the addition of unit flow to a path.Qian and Zhang (2011)conducted a similar analysis for individual path marginal travel times.
Let’s first consider a freeway or an arterial segment with two sequential links, without merges and diverges, say link l 1 and link l. Under congested conditions, there are three basic cases of interest, when the additional unit of flow arrives at this segment.
(i) There is a bottleneck on the downstream link l and the queue on link l does not spill back to link l 1; that is, link l 1 is in free-flow condition while link l is partially congested.
(ii) There is a bottleneck on the downstream link l and the queue on link l spills back to link l 1; that is, link l 1 is par-tially congested while link l is fully congested.
(iii) There is a bottleneck on each of the two links, and the two bottlenecks are independent, assuming that both links are sufficiently long so that the queue in the downstream does not spill back to the upstream. This is in fact the case in which both links l 1 and l are partially congested.
For cases (i) and (iii),Proposition 2can be applied to determine the marginal effects of the additional unit of flow on link flow and density. For case (ii), there are two possible scenarios. As depicted inFig. 6a, one scenario is that the additional unit of flow does not encounter the queue on link l 1, so it can enter link l at time t0
l¼ t00l1¼ t0l1þ FFTTðl 1Þ, when link l is partially congested, i.e.,t0
l<t q l or t0l>t
B
l1. The other scenario, as depicted inFig. 6b, is that the perturbation flow encounters the queue on link l 1, i.e. tqsl1<ðt0
l1þ FFTTðl 1ÞÞ < tBl1, so it cannot enter link l until time t0l¼ t00l1¼ tBl1. Note that, tqs
l1¼ t q
l , the time at which the queue on link l start to spill back to link l 1. While it is possible to detect that the queue spills back based on Eq.(36), tracing the propagation of the perturbation flow in the case of queue spillback is still very dif-ficult. In order to strike a balance between computational efficiency and numerical accuracy, this study appliesProposition 3
to the analysis of link marginal effects on the current link l under fully congested states, and does not trace back the flow perturbation from link l to link l 1.
Next, we discuss the evaluation of partial derivatives of two sequential links in merge or diverge junctions using the illus-trative network depicted inFig. 7. Assume that link b is sufficiently long or consists of several links, so the two junctions are distantly separated. The additional perturbation flow goes through links a1, b, and c1 in order. In the merge junction, the inflow capacity of link b is assumed to be distributed according to the number of lanes of the incoming links (i.e., links a1 and a2), so the perturbation flow from link a1 to link b does not affect the density and flow on link a2. Therefore, the anal-ysis results presented above for two sequential links (i.e., Cases (i)–(iii)) can be directly employed to determine the path mar-ginal effects of link flow and density.
When the perturbation flow arrives at the diverge junction, if there is a bottleneck active on link c1 or a bottleneck active on each of the two links, b and c1, the results discussed above for two sequential links can be applied to evaluate the path marginal effects of interest. Besides, a bottleneck active on link c2 will not affect the evaluation of the path partial derivatives on links b and c1, unless the queue on link c2 spills back to link b. In this case (i.e., there is no inflow capacity on link c2 to accommodate the incoming flow), due to the FIFO constraint, the perturbation flow cannot enter link c1 until the vehicles in front of it in the vertical queue of link b get discharged. We may consider this case as there is a bottleneck at the downstream end of link b, and adopt the results inProposition 2to determine the link/path partial derivatives.
Lastly, we illustrate the evaluation of the partial derivative of path travel time with respect to an unit of additional flow in the same time interval, i.e.,@cðw;@rðw;ss;pÞ;pÞin Eq.(31), based on the case of two partially congested links (i.e., case (iii)), depicted in
Fig. 8. Note that the analysis result can be applied to the other two cases (i and ii). Let link l 1 be the first (congested) link of a path under evaluation. Consider a vehicle veh that enters link l 1 at time t0
l1and leaves at time t00l1. According to
Prop-osition 2, the extra delay experienced by this vehicle on link l 1 is 1/c, due to the perturbation flow entering the link in the same time interval, where c is the capacity of link l 1. Then, veh enters link l at time tin
l and leaves at time toutl , but the im-pact of path flow perturbation on link l occurs at a later time stamp t0
l¼ tBl1>tinl ¼ t00l1. Because @cðw;s;pÞ
@rðw;s;pÞrefers to the marginal effect of the perturbation flow on the path travel time, c(w,s,p), of vehicle veh, c(w,s,p)is affected by the perturbation flow only on link l 1. The perturbation flow, which departs in the same time interval
s
as vehicle veh, will not affect the travel time of veh after link l 1. Therefore, the partial derivative can be approximated as@cðw;s;pÞ@rðw;s;pÞ¼ 1=c.
(a)
(b)
4.3. Computational method for evaluating partial derivatives due to path flow propagation
Based on the above analyses, we now presentAlgorithm 2, which traces the propagation of an additional unit of flow and calculates the gradients in Eqs.(29) and (30)of a given path at departure time
s
. Let tsland telbe the starting and end times of (event) impact period on link l, respectively. Under the event-based mesoscopic traffic simulation framework, te
l can take on one of the three values: ts
lþ FFTTðlÞ; tBl, or t q
l , corresponding to the free-flow condition, congested without queue spillback (or partially congested), and congested with queue spillback (or fully congested), respectively.
Algorithm 2. Computational method for path partial derivatives Initialize link entry time ts
l ¼ departure time
s
, and link index l = 1.Do while link index l 6 L Given ts
l, determine the end time of impact period t e
l as follows:
Step 1: Set the tentative entrance time to the vertical queue as ttemp
l ¼ t
s
lþ FFTTðlÞ. Step 2: (Determining the current traffic condition)
Step 2.1: If there is no vertical queue at time ttempl and ttempl is not equal to the beginning of the queue (tqsl Þ, then the perturbation flow enters link l under uncongested conditions: perform Step 3; otherwise, time ttempl is under con-gested conditions.
Step 2.2: Next, if there is no queue spillback between time ts
lþ FFTTðlÞ and t B
l, then perform Step 4; otherwise, perform Step 5.
Step 3: (Free-flow condition) Step 3.1: te
l ¼ tslþ FFTTðlÞ.
Step 3.2: Update link inflow, outflow and travel time partial derivatives according toProposition 1.
Step 3.3: If link l has the measurements (e.g., q and k), then cumulaterhq(r) andrhk(r), according to Eqs.(29)
and (30).
Step 3.4: Set l = l + 1 and ts l¼ tel1.
Step 4: (Congested condition without queue spillback) Step 4.1: Set te
l ¼ tBl.
Step 4.2: Update link partial derivatives according toProposition 2for impact period ts l; tel
.
Step 4.3: If link l has the measurements (e.g., q and k), then cumulaterhq(r) andrhk(r), according to Eqs.(29) and
(30).
Step 4.4: Set l = l + 1 and ts l ¼ tel1.
Step 5: (Congested condition with queue spillback):
Step 5.1: Find the queue spillback time (in terms of link entrance time) tq
l , and set tel ¼ t q l . Step 5.2: Update link partial derivatives according toProposition 3for impact period ts
l; tel
.
Step 5.3: If link l has the measurements (e.g., q and k,), then cumulaterhq(r) andrhk(r), according to Eqs.(29) and (30).
Step 5.4: Keep l unchanged, set ts l¼ t e l, go back to Step 4. End Do a1 a2 b c1 c2
Fig. 7. An illustrative network with merge and diverge junctions.
Note that, because of the limitation of using local (instead of global) link partial derivatives and the lack of an exact ap-proach to keep track of the propagation of perturbation flows, the gradients obtained usingAlgorithm 2are still numerical approximates. Thus, the proposed solution algorithm, based on the approximate gradients, is a heuristic that cannot be guar-anteed to find (globally) optimal solutions.
5. Numerical experiments
Three sets of numerical experiments were conducted on three different networks to systematically evaluate the perfor-mance of the proposed dynamic OD demand estimation algorithm under different scenarios. The first network is a simple two-link corridor with steady state travel time functions. The second network is a simple freeway corridor with time-depen-dent sensor data. The last network is a real-world transportation network with time-depentime-depen-dent sensor data.
5.1. Experiments on a simple two-link corridor with steady state travel time function
In the first set of experiments, we aim to examine the convergence pattern of the proposed algorithm on a simple corridor with a single OD pair connected by two parallel links (or paths). A simple linear travel time function, Eq.(37), is used in per-forming the traffic assignment which loads a total peak-hour demand, 8000 vehicles/hour (or veh/h) to those two paths. Then, the resulting UE assignment results, shown inTable 1, are used as the ground-truth condition to evaluate the path flow estimation performance under various testing conditions.
TðlÞ ¼ FFTTðlÞ þ rl=capðlÞ; ð37Þ
where T(l) and FFTT(l) are the travel time and free-flow travel time on link/path l, respectively. rland cap(l) are the flow vol-ume and capacity of link/path l, respectively.
First, we start with an initial path flow distribution that loads 3000 veh/h to each link. The ground-truth demand of 8000 veh/h is set as the target demand, and the error-free flow counts (r1= 5400, r2= 2600) are used as the observations. In this case, the objective function of problem P1 can be simplified as follows:
Min LðrÞ ¼ bdf ðrÞ þ bqh q
ðrÞ þ k½gðr;
p
Þ ð38ÞFor simplicity, the weight parameters bd= 1, and bq= 1, which means that the decision maker has the same level of con-fidence on target demand and flow observations. By further considering the Lagrange multiplier, k = 1, the first order partial derivative of the objective function with respect to path flow on path 1 is
@LðrÞ=@r1¼ 2ðr1þ r2 8000Þ þ 2ðr1 5400Þ þ Tð1Þ
p
þ r1=capð1Þ; ð39Þwhere the minimum path travel time
p
is considered as an exogenous variable for this restricted optimization problem, and it can be obtained asp
= Min{T(1), T(2)} in each iteration. The second order partial derivative for path flows r1and r2are 4 + 2/cap(1) and 4 + 2/cap(2), respectively. Similar to the Newton-type algorithm, if the inverse of the second order gradient is used as the step size in Eq.(26), then the step sizec
min inner iteration m for updating the flow on path l can be approx-imated as 1/4, as 2/cap(l) reduces to a very small value. This leads to the following approximate gradient-projection-based flow updating formularmþ1 l ¼ max 0; r m l 1=4 @LðrÞ @rl ð40Þ Figs. 9 and 10demonstrate the convergence patterns of the proposed path flow estimation algorithm in the first 20 inner iterations. We can observe that, after 3 or 4 inner iterations, the total estimated demand is quickly adjusted to a level very close to the ground-truth demand, while the equilibrium processes of path flow distribution and path travel times are relatively slow.
The following experiment was conducted to examine the impact of different values of the Lagrange multiplier k on the solution quality. We consider two different criteria for evaluating the solutions, (i) the total gap g(r,
p
) that measures the distance from the UE conditions and (ii) the value of bdf(r) + bqhq(r) that measures the total deviation from the target demand and observed path flows. In this experiment, a slightly biased target demand with 7000 veh/h and flow observations (r1= 5500, r2= 2500) are adopted. The weight parameters remain as bd= 1, and bq= 1, but the Lagrange multiplier k was var-ied from 0.1 to 10 to obtain different solutions. By solving the optimization problem in Eq.(38)using Microsoft’s Excel Solver,Table 1
User equilibrium traffic assignment results on the two-link corridor.
Path FFTT (min) Capacity (veh/h) Assigned flow (veh/h) Travel time (min)
Path 1 20 3000 5400 56
we can construct the resulting trade-off plot between user equilibrium gap and total deviation, as shown inFig. 11. Essen-tially, a larger Lagrange multiplier/penalty associated with the gap function leads to a solution closer to the UE conditions at the expense of an increased total deviation from the target demand and the observations.
It should be noted that, in order to iteratively search for the maximum value of the Lagrangian dual problem, one can use the rule, described in Eq.(26), to determine the step size
a
(n)in an outer iteration n for updating the Lagrange multiplier k(n+1). However, because the gap function value could vary significantly, e.g., from 0 to 14,000 in the above example, we can use a simple line search method to update the step k(n+1)in this experiment to obtain stable improvement in the search process.As described in Step 2 ofAlgorithm 1, in each outer iteration n, the estimated OD demand and path and link flows (cor-responding to a particular value of the Lagrange multiplier) are used to compute the lower bound (LB) of the primal problem, P1, based on Eq.(33). On the other hand, the upper bound (UB) is obtained by solving the corresponding UE problem. Then, in Step 3, the solution with the smallest gap between the UB and LB values is selected as the final solution. As shown inFig. 12,
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 0 5 10 15 20 Flow Volume Iteration Ground-truth total demand Estimated total demand Ground-truth flow on route 1 Estimated Flow on route 1 Estimated Flow on route 2 Ground-truth flow on route 2
Fig. 9. Path flow volume convergence pattern as a function of inner iteration number.
40 45 50 55 60 65 70 0 5 10 15 20
Travel Time (min)
Iteration
Estimated travel time on route 1 Estimated travel time on route 2 User equilibrium travel time
Fig. 10. Path travel time convergence pattern as a function of inner iteration number.
in this experiment, k = 7 gives the solution with the smallest gap, which corresponds to a total demand of 7353 veh/h, a path flow distribution of (r1= 5012, r2= 2341), and equilibrium travel time of 53.4 min. In this final solution, the resulting relative Lagrangian gap value=(UB-LB)/LB = 0.17%, indicating that a close-to-optimal solution is found to the original primal problem, although the total demand, 7353 veh/h, is still slightly different from the ground-truth OD demand, 8000 veh/h, due to the inherent error in the target demand, 7000 veh/h.
Note that, when varying the Lagrange multiplier from 0.1 to 7, we obtain several similar solutions, with the total demand ranging from 7334 to 7353 veh/h. This demonstrates that the proposed solution framework is able to provide a theoretically rigorous method to quantify the solution quality, generate alternative solutions, and finally identify the optimal solution that can minimize the total measurement deviation and satisfy the UE conditions.
Table 2shows the estimation results under different degrees of information availability, for example, partial vs. complete sensor coverage, slightly biased vs. error target demand, as well with or without travel time observations. In general, more information from either target demand or measurements leads to a demand/path flow estimate closer to the ground-truth value.
5.2. Experiments on a freeway corridor with time-dependent sensor data
In the second set of experiments, we test the performance of the proposed algorithm on a freeway corridor with time-dependent real-world sensor data. As shown inFig. 13, the freeway corridor of interest is a 2-mile section of I-210 west-bound, located in Los Angeles, CA. This corridor includes three on-ramps and one off-ramp. In the network representation constructed for the proposed DNL model, we ignore the HOV lane and only consider four general purpose lanes on the free-way. Traffic speed, flow count and occupancy are measured at 5-min intervals on freeway and ramp links.
In this simple corridor, each OD pair only has a single path, so it does not involve a complex flow equilibration process required for multiple alternative paths. As a result, our focus is on demonstrating how the proposed gradient-based
Fig. 12. Upper bound and lower bound of objective function as a function of Lagrangian multiplier.
Table 2
Estimation results under different degrees of information availability.
Information availability Estimation result
Volume observations on path 1 Volume observations on path 2 Error-free target demand, 8000 veh/h Error-free travel time on path 1 Flow on path 1 Flow on path 2 Total estimated demand Equilibrium travel time (min)
X 5051.7 2367.8 7419.5 53.7
X 4967.7 2311.8 7279.4 53.1
X X 5011.8 2341.2 7353.0 53.4
X X X 5387.9 2592.0 7979.9 55.9
X X X X 5401.1 2600.7 8001.8 56.0
adjustment algorithm adjusts the incoming demand pattern to capture the observed queue formation, propagation and dis-sipation. We first describe the details for preparing the input data for the dynamic OD demand estimation problem under consideration.
5.2.1. Traffic measurements
The traffic measurements were obtained from the PErformance Measurement System (PeMS) database (Varaiya, 2002). This study considers a planning horizon between 6AM and 10AM, and uses the flow count and density data on April 28, 2008 as the measurements, while the density data are constructed based on the observed occupancy data on the same day and estimated average vehicle length. Due to possible sensor malfunction, clean ramp sensor data are not always avail-able for the entire planning horizon, so we use the historical average time series to fill out the data gaps on-ramps. 5.2.2. Traffic flow model
Using the historical data from April 22 to 25, 2008 (i.e., 4 weekdays) during the early peak hours between 6 AM and 10 AM, we calibrate the triangular traffic flow-density model for each link. For instance,Fig. 14shows the flow-density rela-tionship for sensor b with jam density Kjam= 120 veh/ml/lane, free-flow speed Vf= 76 mph, backward wave speed w = 20 mph, and a maximum capacity of 1900 veh/h/lane.
5.2.3. Boundary outflow capacity
One of the critical boundary conditions in this OD estimation problem is the time-dependent maximum outflow dis-charge rate (i.e. capacity) on the downstream sensor d, which is a result of the queue spillback from links further down-stream and cannot be captured internally in the study network. Thus, we use the historical average flow data to construct the link outflow apacity on station d and set it as the fixed input for the DNL model.
5.2.4. Historical demand
Based on the historical flow counts on the boundary sensors, namely, the stations on freeway location a and all on-ramps, we estimate the total origin volume and then apply an estimated destination split to setup a historical OD demand table, where the destination split assigns the majority of flow toward the freeway destination d. This historical demand table serves as the target demand table in the experiment.
5.2.5. Initial demand table
To evaluate if and how the gradient-based algorithm can adjust a biased, initial OD demand table toward the target de-mand table, we start with an OD dede-mand table with only 70% of historical dede-mand volume.
Fig. 15shows the observed time series of the traffic speeds at stations b, c and d. After 6:30 AM, the traffic congestion from downstream location d is propagated to upstream sensors b and c, leading to a slow-moving queue on freeway mainline seg-ments between 7:00 AM and 9:00 AM. The traffic congestion starts to dissipate afterwards, but the overall speed is still sig-nificantly below the free-flow level. The observed flow pattern at entrance station a can be viewed as the total origin demand flow from station a to different destinations. The root mean squared absolute errors (RMSEs) of the estimation results for speed measurements and for link count measurements are 13.2 mile/h and 159.1 veh/h/lane, respectively.Fig. 16shows the estimated and observed flow patterns on entry link a, and the corresponding average relative estimation errors are less than 10%. This indicates that the proposed flow estimation algorithm can adjust a biased, initial demand pattern to match the target demand volume at the entry point. The estimated space speeds and the observed point speeds at station c are plotted inFig. 17, which demonstrates that the DNL model is able to accurately reproduce the queue spillback phenomenon along the corridor.
5.3. Experiments on large real-world traffic networks
The third set of numerical experiments was performed on two real-world networks. The first is a subarea network within the Portland, Oregon metropolitan area, which includes 858 nodes, 2000 links, and 208 OD zones shown inFig. 18. We first loaded a subarea OD demand matrix generated from a recent study (Kittelson et al., 2011) as the ‘‘historical’’ OD demand, and obtain 5-min flow counts from 14 loop detectors and peak-hour flow counts converted from 392 annual average daily traffic (AADT) counts. To assess the estimation performance of the proposed approach, we used the root mean squared abso-lute error (RMSE) of peak-hour link volume as the estimation criterion.Fig. 19shows a decreasing pattern of the estimation error as the iterative algorithm proceeds. The final RMSE value of peak-hour link volume is about 55.6 veh/link, correspond-ing to 15% error compared to the total link volume and R2= 0.91. The average travel time gap reduces to 0.33 min/veh on a network with an average travel time of 7 min, which demonstrates that the proposed algorithm is able to produce accurate estimation results for this medium-scale network.
The second real-world network is the Triangle Regional Model network, North Carolina, USA. This large-scale regional network has 2389 zones, 20,259 links and about 2000 signalized intersections. Provided by the local metropolitan planning agency, the morning peak-hour demand matrix has about 1.06 million vehicles, covering a time period from 6 AM to 10 AM. There are about 16 and 14 sensors, respectively, on freeway and arterial links, producing a total of 120 hourly link count observations. The experiments were performed on a PC with 16 GB memory and 8-core processors running at 2.70 GHz. We implemented the proposed path flow adjustment algorithm in the open-source dynamic traffic assignment package DTA-Lite (Zhou et al., 2012), which has a mesoscopic traffic network loading model based on Newell’s simplified kinematic wave theory. The starting and ending timestamps of vehicles on each link are used to construct link-based cumulative arrival and departure count curves, which are further used to track the queue spillback phenomena and shock wave propagation.
Fig. 15. Observed speed time series on freeway stations (demonstrating a congested period from 7:00 AM to 9:00 AM).
Fig. 16. Observed lane volume on station a vs. estimated lane volume on entrance link.
