Statistical inference of probabilistic origin-destination demand using day-to-day traffic data

TL;DR

This paper introduces a novel statistical framework for estimating probabilistic origin-destination demand and network flow distributions using extensive day-to-day traffic data, accounting for various sources of variability.

Contribution

It develops an iterative generalized least squares method with Lasso regularization to estimate O-D demand mean and covariance, improving interpretability and computational efficiency.

Findings

Framework accurately estimates O-D demand mean and variance.

Method converges quickly on multiple network scales.

Lasso regularization balances bias and variance in covariance estimation.

Abstract

Recent transportation network studies on uncertainty and reliability call for modeling the probabilistic O-D demand and probabilistic network flow. Making the best use of day-to-day traffic data collected over many years, this paper develops a novel theoretical framework for estimating the mean and variance/covariance matrix of O-D demand considering the day-to-day variation induced by travelers' independent route choices. It also estimates the probability distributions of link/path flow and their travel cost where the variance stems from three sources, O-D demand, route choice and unknown errors. The framework estimates O-D demand mean and variance/covariance matrix iteratively, also known as iterative generalized least squares (IGLS) in statistics. Lasso regularization is employed to obtain sparse covariance matrix for better interpretation and computational efficiency. Though the probabilistic O-D estimation (ODE) works with a much larger solution space than the deterministic ODE, we show that its estimator for O-D demand mean is no worse than the best possible estimator by an error that reduces with the increase in sample size. The probabilistic ODE is examined on two small networks and two real-world large-scale networks. The solution converges quickly under the IGLS framework. In all those experiments, the results of the probabilistic ODE are compelling, satisfactory and computationally plausible. Lasso regularization on the covariance matrix estimation leans to underestimate most of variance/covariance entries. A proper Lasso penalty ensures a good trade-off between bias and variance of the estimation.

Figures

Tables

Table 1. Table 1 : List of notations

Network Variables
$A$	The set of all links
$A^o$	The set of links with flow observations
$N$	The number of links
$K_q$	The set of all O-D pairs
$K_{rs}$	The set of all paths between O-D pair $r-s$
$\Delta$	Path/link incidence matrix
${\Delta }^o$	Path/observed link incidence matrix
$M$	Path/O-D demand incidence matrix
Random Variables
$Q_{rs}$	The demand of O-D pair $r-s$
$Q$	The vector of O-D demands
$X_a$	Link flow on link $a$
$X$	The vector of link flow
$X^o$	The vector of link flow on observed links
$X_m$	The vector of measured link flow
$F_{rs}^k$	Path flow on path $k$ between O-D pair $r-s$
$F_{rs}$	The vector of path flow between O-D pair $r-s$
$F$	The vector of path flow between all O-D pairs
$C$	The vector of path costs between all O-D pairs
$E$	The vector of unknown error for links
Parameters for Probability Distributions
$q=𝔼\left(Q\right)$	The vector of means of O-D demands
$q_{rs}=𝔼\left(Q_{rs}\right)$	The mean of O-D demand $Q_{rs}$
${\Sigma }_q$	Covariance matrix of O-D demands
$p_{rs}^k$	The probability of choosing path $k$ in all paths between O-D pair $r-s$
$p_{rs}$	The vector of route choice probabilities of all paths between O-D pair $r-s$
$p$	Route choice probability matrix, consisting of all $p_{rs}^k$
$c=𝔼\left(C\right)$	The vector of means of path costs for all O-D pairs
${\Sigma }_c$	Covariance matrix of path costs vector
$x=𝔼\left(X\right)$	The vector of means of link flow
$x_a=𝔼\left(X_a\right)$	The mean of link flow $X_a$
${\Sigma }_x$	Covariance matrix of link flow
$f=𝔼\left(F\right)$	The vector of means of path flow
$f_{rs}^k=𝔼\left(F_{rs}^k\right)$	The mean of path flow $F_{rs}^k$
${\Sigma }_f$	Covariance matrix of path flow
${\Sigma }_e$	Covariance matrix of unknown error
Parameters for Conditional Probability Distributions
${\Sigma }_{f|Q}$	Covariance matrix of path flow conditional on O-D demands $Q$
${\Sigma }_{f|q}$	Covariance matrix of path flow conditional on O-D demands $Q=q$, ${\Sigma }_{f|q}={\Sigma }_{f|Q=q}$
Observed Link Flow Data
${𝗑}^o$	Observed link flow matrix
$n$	Number of observed link flow vectors, namely number of observed days
${𝗑}_i^o$	The $i$-th observed link flow vector (on the $i$-th day), $1\le i\le n$
$S_x^o$	Empirical covariance matrix of the observed link flow

Table 2. Table 2 : Results of probabilistic ODE on the four-link toy network (no historic O-D demand information is used)

True $\rho$	Settings	${\hat{q}}_{1\to 3}$	${\hat{q}}_{2\to 3}$	${\hat{\sigma }}_{1\to 3}^2$	${\hat{\sigma }}_{2\to 3}^2$	$\hat{\rho }$	RMPSE	KL-distance
	True value	$700$	$500$	$175$	$125$	NA	NA	NA
$0.5$	w/o EC - w/o Lasso	$722.17$	$500.41$	$186.69$	$134.21$	$0.56$	$3.62\%$	$3.64$
	Logit - w/o Lasso	$682.36$	$499.63$	$207.94$	$134.21$	$0.50$	$2.08\%$	$1.17$
	Probit - w/o Lasso	$699.50$	$499.63$	$200.94$	$134.21$	$0.52$	$0.07\%$	$0.01$
$0$	w/o EC - w/o Lasso	$715.91$	$500.46$	$143.05$	$138.74$	$0.03$	$1.87\%$	$0.74$
	Logit - w/o Lasso	$681.28$	$500.46$	$162.49$	$138.75$	$0.02$	$2.21\%$	$1.01$
	Probit - w/o Lasso	$700.30$	$500.46$	$152.15$	$138.75$	$0.03$	$0.06\%$	$0.01$
	Logit - w/ Lasso	$681.28$	$500.46$	$144.52$	$128.75$	$0.00$	$2.21\%$	$1.01$
	Probit - w/ Lasso	$700.02$	$500.46$	$132.27$	$128.75$	$0.00$	$0.05\%$	$0.004$
$-0.5$	w/o EC - w/o Lasso	$703.41$	$499.06$	$173.34$	$132.60$	$-0.41$	$0.43\%$	$0.04$
	Logit - w/o Lasso	$681.05$	$499.06$	$184.13$	$132.60$	$-0.39$	$2.23\%$	$1.47$
	Probit - w/o Lasso	$701.71$	$499.06$	$174.19$	$132.60$	$-0.41$	$0.23\%$	$0.02$

Equations

$$(X_1,X_2{)}^T$$

$$Q$$

$$\text{mbox}\mathrm{Var}(X_1)$$

$$100$$

$$Q=\left[Q^{\prime }\right],Q^{\prime }\sim \mathcal{N}(q,{\Sigma }_q)$$

$$E\sim \mathcal{N}(\text{mbox}0,{\Sigma }_e)$$

$$p_{rs}^k=\mathrm{Pr}\left[C_{rs}^k\le \underset{i\ne k}{\min }{C}_{rs}^i\right]$$

$$Level1:$$

$${ϵ}_e\sim \mathcal{N}(\text{mbox}0,{\Sigma }_e)$$

$$Level2:$$

$$F\sim \mathcal{M}\mathcal{N}(\tilde{p}Q,{\Sigma }_f)$$

$$Level3:$$

$$F\sim \mathcal{N}(f,{\Sigma }_f)$$

$$X\sim \mathcal{N}(x,{\Sigma }_x)$$

$$X_m\sim \mathcal{N}(x,{\Sigma }_x+{\Sigma }_e)$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\min_{q,\Sigma_{q}}&\lx@intercol\displaystyle\alpha_{1}f_{1}(x(q),x^{o})+\alpha_{2}f_{2}(\Sigma_{x}(\Sigma_{q}),\Sigma_{x}^{o})\hfil\lx@intercol\\ \end{array}$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\min_{q,\Sigma_{q}}&\lx@intercol\displaystyle\alpha_{1}f_{1}(x(q),x^{o})+\alpha_{2}f_{2}(\Sigma_{x}(q,\Sigma_{q}),\Sigma_{x}^{o})\hfil\lx@intercol\\ \end{array}$$

$$\tau$$

$$D_H(({\mu }_1,{\Sigma }_1{)}^T,({\mu }_2,{\Sigma }_2{)}^T)$$

$$D_{KL}(({\mu }_1,{\Sigma }_1{)}^T,({\mu }_2,{\Sigma }_2{)}^T)$$

$$l({\mathsf{x}}^o)=\prod _{i=1}^n\frac{1}{\sqrt{(2\pi {)}^{|A|}|{\Sigma }_x^o|}}\exp \left(-\frac{1}{2}(x^o-{\mathsf{x}}_i^o{)}^T{\left({\hat{\Sigma }}_x^o\right)}^{-1}(x^o-{\mathsf{x}}_i^o)\right)$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\min_{x^{o}}&\lx@intercol\displaystyle\frac{1}{2}\sum_{i=1}^{n}(x^{o}-{\mathsf{x}}_{i}^{o})^{T}\left(\hat{\Sigma}_{x}^{o}\right)^{-1}(x^{o}-{\mathsf{x}}_{i}^{o})\hfil\lx@intercol\\ \textrm{s.t.}&x^{o}&\geq&0\\ \end{array}$$

$${\hat{x}}^o\sim \mathcal{N}(x^o,\frac{1}{n}{\Sigma }_x^o)$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\min_{q}&\lx@intercol\displaystyle n\left(\Delta^{o}\tilde{p}q-\hat{x}^{o}\right)^{T}\left(\hat{\Sigma}_{x}^{o}\right)^{-1}\left(\Delta^{o}\tilde{p}q-\hat{x}^{o}\right)\hfil\lx@intercol\\ \textrm{s.t.}&q&\geq&0\\ \end{array}$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\min_{q}&\lx@intercol\displaystyle\left\lVert\Theta_{x}(\Delta^{o}\tilde{p}q-\hat{x}^{o})\right\rVert_{2}^{2}\hfil\lx@intercol\\ \textrm{s.t.}&q&\geq&0\\ \end{array}$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\min_{q}&\lx@intercol\displaystyle n\left(\Delta^{o}\tilde{p}q-\hat{x}^{o}\right)^{T}\left(\hat{\Sigma}_{x}^{o}\right)^{-1}\left(\Delta^{o}\tilde{p}q-\hat{x}^{o}\right)+(q^{H}-q)^{T}\left(\Sigma_{q}^{H}\right)^{-1}(q^{H}-q)\hfil\lx@intercol\\ \textrm{s.t.}&q&\geq&0\\ \end{array}$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\min_{f}&\lx@intercol\displaystyle n\left(\Delta^{o}f-\hat{x}^{o}\right)^{T}\left(\hat{\Sigma}_{x}^{o}\right)^{-1}\left(\Delta^{o}f-\hat{x}^{o}\right)+(q^{H}-Mf)^{T}\left(\Sigma_{q}^{H}\right)^{-1}(q^{H}-Mf)\hfil\lx@intercol\\ \textrm{s.t.}&f&\geq&0\\ \end{array}$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\min_{f}&\lx@intercol\displaystyle n\left(\Delta^{o}f-\hat{x}^{o}\right)^{T}\left(\hat{\Sigma}_{x}^{o}\right)^{-1}\left(\Delta^{o}f-\hat{x}^{o}\right)+(q^{H}-Mf)^{T}\left(\Sigma_{q}^{H}\right)^{-1}(q^{H}-Mf)\hfil\lx@intercol\\ \textrm{s.t.}&f&\in&\Phi^{+}\\ \end{array}$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle Z_{1}(f)=&\lx@intercol\displaystyle\frac{1}{\Theta}\sum_{a}\int_{0}^{x_{a}}t_{a}(w)dw\hfil\lx@intercol\\ \textrm{\text{where}}&x&=&\Delta f\\ &f&\geq&0\\ \end{array}$$

$${\Phi }_{\text{UE}}^{+}=\{f_1|Z_1(f_1)\le Z_1(f_2),\forall f_2\ge 0\text{ such that }M{f}_1=M{f}_2\}$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle Z_{2}(f)=&\lx@intercol\displaystyle\frac{1}{\Theta}\sum_{ts}\sum_{k}f_{rs}^{k}\log(f_{rs}^{k})\hfil\lx@intercol\\ \textrm{\text{where}}&f&\geq&0\\ \end{array}$$

$${\Phi }_{\text{Logit}}^{+}=\{f_1|Z_2(f_1)\le Z_2(f_2),\forall f_2\ge 0\text{ such that }M{f}_1=M{f}_2\}$$

$$\begin{array}[]{rrclcl}\vspace{5pt}\displaystyle\max_{\Sigma_{q}}&\lx@intercol\displaystyle\log\det(\left(\Sigma_{x}^{o}\right)^{-1})-\text{trace}(S_{x}^{o}\left(\Sigma_{x}^{o}\right)^{-1})\hfil\lx@intercol\\ \textrm{s.t.}&\Sigma_{x}^{o}&=&\Delta^{o}\Sigma_{f|q}\left(\Delta^{o}\right)^{T}+\Delta^{o}\tilde{p}\Sigma_{q}\tilde{p}^{T}\left(\Delta^{o}\right)^{T}\\ &\Sigma_{q}&\succeq&0\end{array}$$

Full text

1 Introduction

Origin-destination (O-D) demand is a critical input to system modeling in transportation planning, operation and management. For decades, O-D demand is deterministically modeled, along with deterministic models of link/path flow and travel cost/time in classical traffic assignment problems. Transportation network uncertainty and reliability call for modeling the stochasticity of O-D demand, namely its spatio-temporal correlation and variation. With the increasing quantity and quality of traffic data collected years along, it is possible to learn the stochasticity of O-D demand for better understanding stochastic travel behavior and stochastic system performance metrics. Some studies considered the stochastic features of O-D demand, but few estimated the mean and variance of O-D demand from day-to-day traffic data. What is missing in the literature is the capacity to estimate spatially correlated multivariate O-D demand, in conjunction with a sound network flow theory on probabilistic route choices that can be learned from day-to-day traffic data. In view of this, this paper develops a novel data-friendly framework for estimating the mean and variance/covariance of O-D demand based on a generalized statistical network equilibrium. The statistical properties towards the estimated probabilistic O-D demand are also analyzed and provided. The process of O-D demand estimation (ODE) is further examined in real-world networks for insights.

O-D estimation (ODE) requires an underlying behavioral model, based on which O-D demand is estimated such that it best fits observations. Behavioral models in the static network context, namely route choice models, are also known as static traffic assignment models. The classical traffic assignment models [e.g., 21, 64] deterministically map the deterministic O-D demand $q\in\mathbb{R}^{R\times S}_{+}$ to link flow $x\in\mathbb{R}^{N}_{+}$ or path flow $f\in\mathbb{R}^{K}_{+}$ (where $R$, $S$, $N$ and $K$ are the cardinality of sets of origins, destinations, links, and paths, respectively). In fact, the deterministic O-D demand $q$ is assumed to represent the mean number of vehicles in the same peak hour from day to day. Likewise, link (path) flow is also deterministic, representing the mean number of vehicles on a link (path) in the same hour from day to day. The classical traffic assignment models (such as User Equilibrium and Stochastic User Equilibrium) lay out the foundation of deterministic O-D estimation methods, namely to estimate $q$ in a way to best fit observed data related to a subset of link/path flow $x,f$. Deterministic O-D estimation (ODE) include the entropy maximizing models [81], maximum likelihood models [76, 88], generalized least squares (GLS) models [10, 6, 94, 85], Bayesian inference models [50, 77] and some recent emerging combined models [3, 15, 44, 16, 46, 91]. For more details, readers are referred to the comprehensive reviews by Bera and Rao [7], Castillo et al. [12].

Classical traffic assignment models and ODE methods overlook the variance/covariance of demand and link/path flow, an essential feature for network traffic flow. Recent studies on network reliability and uncertainty model stochasticity of the network flow. Since ODE requires a traffic assignment model, we first review statistical traffic assignment models, followed by ODE models that take into account stochasticity.

One aspect of the stochastic network flow is using probability distributions to represent O-D demand. For example, Waller et al. [84], Duthie et al. [25] sampled O-D demand from given multivariate normal distributions (MVN) and evaluated the network performance under classical User Equilibrium (UE) condition. Chen et al. [18] used a similar simulation-based method to evaluate travelers’ risk-taking behavior due to probabilistic O-D demand. All these studies indicate the O-D variation is of great importance to network modeling and behavioral analysis. Statistical traffic assignment models consider various O-D probability distributions such as Poisson distributions [20], MVN [14], multinomial distributions [52]. Nakayama and Watling [55] summarized different formulations and proposed a unified framework for stochastic modeling of traffic flows. Advantages and disadvantages of modeling traffic with those probabilistic distributions are also discussed by Castillo et al. [17]. Shao et al. [70, 71] proposed a reliability-based traffic assignment model (RUE) and extended it to consider different travelers’ risk taking behavior. Lam et al. [41] further extended the model to consider the traffic uncertainty and proposed reliability-based statistical traffic equilibrium. Zhou and Chen [100], Chen and Zhou [19] proposed a $\alpha$-reliable mean-excess traffic assignment model which explicitly models the travel time distribution and consider the reliability and uncertainty of the travel time on travelers’ route choice behavior. Other studies [31, 83, 25] show that the variance/covariance matrix of the O-D demand have a significant influence on network traffic conditions.

Though adopting stochastic O-D demand, those traffic assignment models (except Clark and Watling [20]) assumed non-atomic (infinitesimal) players and therefore are unable to capture the stochasticity of route choices that vary from day to day (the proof is shown in Ma and Qian [49]). Classical UE, SUE and RUE are all deterministic route choice models where the number of (infinitesimal) players assigned to each route is fixed, rather than being stochastic. Thus, those models are unable to explain the day-to-day variation of observed traffic counts at the same location and the same time of day. To further see how classical equilibrium models overlook the route choice stochasticity, suppose there is $Q$ travelers where $Q$ is a random variable to be realized on each day. Given the probability of choosing a route $p$, the route flow is deterministically identified by the number of infinitesimal users who take this route, $F=pQ$. Even if the route choice probability $p$ is determined by stochastic choice models (such as Probit, Logit, etc.), these models still assume that a fixed number $Qp$ of travelers take this route on each day, as a result of non-atomic equilibrium. This does not, theoretically, allow to model the day-to-day variations of travelers’ choices. Recent studies on statistical traffic assignment models indicate that the route flow is the aggregation of random choices of O-D demand, and thus also random [86, 87, 53, 54, 55]. Travelers’ route choice follows a multinomial distribution with the probability obtained from route choice models, $f\sim\text{Multinomial}(Q,p)$. To distinguish to what extent stochasticity is modeled for route choices, we refer to the former (classical) models as “*fixed portions with stochastic route choice models *” and the latter models as “probabilistic distributions with stochastic route choice models”. Though route choices are stochastic, these studies did not work directly with the covariance of demand among all O-D pairs. A detailed comparison of those assignment models is further illustrated in Ma and Qian [49].

Given a route choice probability, we can derive probability distributions of path/link flow. However, it is non-trivial to establish a statistical network equilibrium where the route choice probability is determined endogenously as a result of stochastic O-D demand, path/link flow and network conditions. In the deterministic settings, UE, SUE or RUE simultaneously determines the mean path/link flow, and the route choice probability [e.g., 94]. Very few studies examined the statistical network equilibrium. Davis and Nihan [23] proposed a Markov process to model the day-to-day variation of traffic flow given a fixed regional population. The stochastic route choice on a particular day is assumed to be related to the stochastic network condition of the previous days [11]. When the network evolves from day to day, there exists a network equilibrium where stabilized probability distributions of path/link flow are reached. Different from this approach, Ma and Qian [49] proposed a generalized statistical equilibrium where each traveler makes stochastic choices based on his/her entire past experience, namely the probability distributions of the equilibrated network conditions. Ma and Qian [49] integrates multivariate probability distributions of O-D demands and link/path flow into the stochastic route choice models, and ultimately solved for the probabilistic network flow. It also analytically decomposes the variance of link/path flow into three sources, O-D demand variation, route choice variation and measurement error.

With little work on statistical traffic assignment models, probabilistic ODE can be challenging. An ideal probabilistic ODE should posses three critical features: estimating variance/covariance of O-D demand in addition to its mean, a statistical network equilibrium that can fit massive data collected years along, and consideration of day-to-day route choice variation. To our best knowledge, few literature works with ODE models that take into account any of the three features. We summarize existing representative ODE models in Figure 1.

To further distinguish our work from existing literature, Vardi [82], Hazelton [32, 33, 34], Li [43], Parry and Hazelton [62] assumed O-D demand follows Poisson distributions that are independent among O-D pairs, and Hazelton [32, 33, 34], Li [43], Parry and Hazelton [62] further considered day-to-day route choice variation. They formulated maximum likelihood estimator and Bayesian inference method for O-D demand, whereas the demand covariance among O-D pairs is not considered. In addition, using uncongested networks simplifies the route choice model, and thus no network equilibrium under congestion is proposed. Shao et al. [68], Yang et al. [96] proposed a generalized model to estimate the mean and variance of O-D demand with MVN. Shao et al. [69] further extended the model to estimate multi-class O-D demand and used $L^{1}$ regularizer to shrink the model dimensions. An equilibrium is used to model deterministic route choices (similar to classical UE and RUE), while neither of them considered day-to-day route choice variation.

To our best knowledge, there is a lack of study that estimates correlated multivariate probabilistic O-D demand under a mathematically sound statistical network equilibrium (namely a truly stochastic route choice model) while considering day-to-day random route choices simultaneously. To fill up this gap, this paper builds an ODE model based on the generalized statistical traffic assignment model (GESTA) proposed by Ma and Qian [49]. Any observation (such as link/path flow, or travel time/cost) has a variance that stems from three sources, O-D demand variation, route choice variation and unknown errors, all from day to day. Using GESTA as the underlying behavioral model, we estimate probability distributions of O-D demand using data from various data sources. Furthermore, conventional goodness-of-fit indicators [e.g., 2] is not suitable for probabilistic ODE. This paper also proposes new goodness-of-fit indicators based on probability distributions to evaluate the performance of probabilistic ODE.

Another important issue for ODE is the observability problem [13, 95]. It is well known that the ODE is underdetermined using observed link-based traffic counts [94]. Hazelton [35] suggested to use second-order statistical information of the observed link counts to estimate Poisson distributed O-D demand without covariance among O-D pairs. By utilizing the second-order information, the Poisson distributed O-D demand can be estimated uniquely. Yang and Fan [95] proposed to uniquely determine the O-D demand by properly selecting data among observed link flow data and historical O-D information. Yang et al. [96] also showed that possibly dissimilar multi-day observation improves the observability of OD demand [96]. Hazelton et al. [36] proposed new method for sampling latent route flows conditional on the observed link counts when the network’s link-path incidence matrix is unimodular. In this paper, the observability of the proposed probabilistic ODE is examined. The new ODE guarantees that its estimator for the O-D demand mean is no worse than the conventional deterministic ODE by the order of ${\cal O}(\frac{1}{n})$, provided with sufficient data.

Under the proposed probabilistic ODE framework, we propose to estimate the O-D demand mean and variance-covariance matrix iteratively, which decomposes the complex ODE into two sub-problems. In the sub-problem of estimating O-D demand mean vector $q$, we extend the statistical ODE from Menon et al. [51], and both single level and bi-level ODE formulations are discussed. In the sub-problem of estimating O-D demand variance-covariance matrix $\Sigma_{q}$, we utilize estimated link flow covariance to formulate the estimation problem, and then apply Lasso regularization and convex relaxation on the formulation. How to use traffic speed data in addition to traffic counts data is discussed. Furthermore, the statistical properties of the estimated probabilistic O-D demand are provided for insights. The observability of the ODE problem is also examined.

The main contributions of this paper are summarized as follows:

It proposes a novel theoretical framework for estimating probabilistic O-D demand (namely mean and variance/covariance of O-D demand) considering newly defined generalized statistical network equilibrium. The statistical equilibrium simultaneously integrates probabilistic O-D demand and travelers’ day-to-day route choice variation. 2. 2)

It defines the goodness-of-fit for probabilistic ODE and develops a theory for variance analysis of estimated probabilistic O-D demand and path/link flow. 3. 3)

It discusses the observability issue for probabilistic O-D demand, and shows that the estimated mean of O-D demand using the the probabilistic ODE is no worse than the estimate O-D demand using a deterministic ODE model, provided with sufficient data. 4. 4)

It intensively examines the proposed probabilistic ODE framework on two large-scale real networks using both simulated traffic data and real world traffic counts to gain insights from solutions, as well as to show the computational efficiency of the solution algorithms.

The remainder of this paper is organized as follows. Section 2 presents an illustrative example to further explain how and why it is necessary to consider day-to-day route choice variation. Section 3 discusses the formulation details, followed by section 4 presenting properties of the model and addressing the observability issues. Section 5 proposes the entire probabilistic ODE framework. In section 6, two simple illustrative examples are used to demonstrate the concepts and ODE results. Two large scale networks are used to demonstrate the computational efficiency of the probabilistic ODE method, and its ability to work with real-world data. Finally, conclusions are drawn in section 7.

2 An illustrative example

Observations of link/path flow varies from day to day. In principle, the day-to-day variation is attributed to three sources, O-D demand variation, route choice variation and sensing measurement error. The impact of day-to-day O-D demand variation and measurement error on link/path flow have been thoroughly discussed by Yang et al. [94], Waller et al. [83], Shao et al. [68]. Here we use a toy example to compare probabilistic O-D estimation results with and without the consideration of day-to-day route choice variation. It is intended to illustrate that where the day-to-day route choice variation comes from, and why it is important to not neglect it when estimating O-D demand.

Consider a toy network as shown in Figure 2, on each day $Q$ amount of vehicles depart from node $r$ to node $s$. We assume $Q$ is normally distributed, so is the link flow (i.e. path flow in this example) $X_{1}$ and $X_{2}$, the number of vehicles on link 1 and link 2, respectively. Suppose both links are indifferent. We do not consider measurement error here.

On each day, the demand $Q$ is realized, but cannot be directly measured. Instead, we measure the link flow $X_{1}$ on each day. Suppose after observations of many days, we determine the probability distribution of $X_{1}\sim{\cal N}(50,10^{2})$. If we do not consider the day-to-day route choice variation, the estimated O-D demand follows the rounded normal distribution $Q=\left[Q^{\prime}\right],Q^{\prime}\sim{\cal N}(100,20^{2})$ since $Q\simeq 2X_{1}$ by equilibrium conditions. In fact, regardless of which route choice models (logit, probit, etc.), the distribution of the O-D demand can be directly computed since both links are indifferent to travelers. Both links have the same flow distributions.

Now consider travelers’ route choice that varies from day to day. The probability of any vehicle choosing link $1$ and link $2$ is $p_{1}=p_{2}=0.5$. All $Q$ vehicles make this choice independently on each day. Thus, the link flow follows a multinomial distribution (binomial distribution in this example),

[TABLE]

Where $q$ denotes the mean of the probability distribution of $Q$. Since $q=2{\mathbb{E}}(X_{1})=100$, by law of total variance we have,

[TABLE]

Thus $\mbox{\rm Var}(Q)=300$. Clearly, the probabilistic O-D estimation without considering day-to-day route choice variation overestimates the variance of O-D demand by the term $p_{1}p_{2}\mathbb{E}(Q)$, namely $33.3\%$, a fairly substantial quantity. Intuitively, probabilistic route choices that are made by travelers independently on a daily basis tend to reduce the day-to-day variation of route flow. As shown in Ma and Qian [49], classical Stochastic User Equilibrium theory works with non-atomic (infinitesimal) users, which in theory does not allow day-to-day choice variation. This greatly simplifies the computation. However, it potentially conflicts with real-world observations (such as counts and speeds) that vary from day to day, and thus might be difficult to handle massive data collected over many years. Provided that the route choice variation can be substantial in real-world traffic, our ODE method needs to model atomic users, and take day-to-day route choice variation into consideration while estimating probability distributions of O-D demand.

3 Formulations

In this section, we discuss the probabilistic ODE framework. We first present the notations and assumptions, compare this framework to existing formulations, and finally discuss each component of the framework in details.

3.1 Notations

Please refer to Table LABEL:tab:notation. The superscript $\cdot^{o}$ indicates that the variable is projected onto observed links only. The hat symbol, $\hat{\cdot}$, indicates the variable is an estimator for the true (unknown) variable.

3.2 Assumptions

O-D demands follow a rounded multivariate normal (MVN) distribution with mean $q$ and covariance matrix $\Sigma_{q}$.

[TABLE]

where the random variable $Q$ are generated from a standard continuous MVN $Q^{\prime}$ and then rounded to the nearest integer. When the number of travelers is sufficiently large, the rounding error is negligible [52]. In this paper demand is approximated by a continuous MVN. Because the demand is always non-negative, the MVN is truncated at zero. When the O-D demands are sufficiently large, the effect of the truncation is also negligible. 2. 2.

All travelers between the same O-D pair are homogeneous when making route choices. They perceive the probability distribution of travel time/cost of the network from their entire past experience, similar to Nakayama and Watling [55]. On each day, each atomic traveler independently and identically makes a route choice. 3. 3.

Each traveler makes her route choice decision solely based on the perception of the past traffic conditions, unlike the classical Wardrop UE or atomic Nash Equilibrium where each traveler is fully aware of others’ choices. 4. 4.

The day-to-day variation of observed path/link traffic flow is resulted from O-D demand variation, travelers’ route choice variation and unknown errors. Unknown errors include measurement errors, and other unobserved error (such as traffic incident and non-recurrent traffic behavior). Unknown errors are link-based and have zero mean, which also follow MVN.

[TABLE] 5. 5.

The observed link flow vector ${\mathsf{x}}_{i}^{o}$ on the $i$-th day ($1\leq i\leq n$) is an i.i.d sample from the probability distribution of the observed link flow vector $X_{m}^{o}$.

3.3 Review of the generalized statistical traffic assignment (GESTA) model

As the underlying model, a generalized statistical traffic assignment model (GESTA) is proposed in Ma and Qian [49]. We briefly review the GESTA, and then propose the ODE framework based on GESTA.

Consider a graph with $|N|$ nodes, $|A|$ links, $|K_{q}|$ O-D pairs and $|K_{rs}|$ paths for each O-D pair $r-s$, GESTA maps O-D demands $Q$ to link/path flows $X,F$ and path costs $C$ under a statistical equilibrium. From day to day, the recurrent traffic reaches a statistical equilibrium among $Q,C,F,X$ such that, provided with exogenous demand $Q$ (a random variable vector), $C,F,X$ follow stabilized probability distributions that can be represented by a stabilized route choice probability vector $p$ [49]. The statistical equilibrium is defined as follows,

Definition 1**.**

[49]** A transportation network is under a statistical equilibrium, if all travellers practice the following behavior: on each day, each traveler from origin $r$ to destination $s$ independently chooses a route $k$ with a deterministic probability $p_{rs}^{k}$. For a sufficient number of days, this choice behavior leads to a stabilized probability distribution of travel costs with parameters $\vartheta$. This stabilized probability distribution, in turn, results in the probabilities $p=\psi(\vartheta)$ where $\psi(\cdot)$ is a general route choice function.

Note that the definition of statistical equilibrium is independent of time and describes the equilibrated condition. The definition indicates that the traffic flow/costs on each day are random and by no means identical from day to day. However, the probability distributions of traffic conditions over a course of a number of days (such as months or years) are stabilized. Based on the O-D demand $Q$, we express the probability distributions of $X,F,C$ all by the deterministic probability vector $p$, and construct a fixed point problem regarding $p$ to solve for $X,F,C$.

Note that the route choice model $\psi(\vartheta)$ can be generic. For instance, $\psi(\vartheta)$ can represent the classical UE where $\vartheta\doteq\mathbb{E}(C)=c$. If we use a random utility model as the route choice model, then the route choice probability becomes

[TABLE]

Each traveler perceives his/her stochastic costs on all possible paths between O-D pair $r-s$, and the probability of choosing path $k$ equals the probability of path costs outperforming costs of other paths. For this reason, $\psi(\vartheta)$ can also easily represent the Logit model and Probit model. However, Probit model would be the most convenient and natural representation. This is because $C$ is approximated by MVN under GESTA, and this is consistent with the normally distributed perceived travel cost under Probit. If the Logit model is adopted for $\psi(\vartheta)$, then one needs to assume that travelers’ perception on travel time/cost would add a term that follows Gumbel distribution in addition to the cost mean $\mathbb{E}(C)=c$.

The hierarchical formulation of GESTA is represented in Equation 9.

[TABLE]

Level $3$ represents the O-D demand variation. The Level $2$ formulation indicates that each traveler makes their route choice independently, thus the path flow vector $F$ follows a multinomial distribution. In Level $1$, the unknown errors are applied to the link flow vector. $X_{m}$ is the random variable implying the observations of link flow.

Define $\tilde{p}\doteq diag(p)B$ and $B$ is a transition matrix (with blocks of ones and zeros) defined in Watling [86], Ma and Qian [49]. After some derivations, the path/link flow distributions can be obtained, which are stated in the following two propositions.

Proposition 1**.**

The marginal distribution of $F$ can be approximated by,

[TABLE]

where $f=\tilde{p}q$, $\Sigma_{f}=\Sigma_{f|q}+\tilde{p}\Sigma_{q}\tilde{p}^{T}$, $\Sigma_{f|q}$ is the covariance matrix of path flow conditional on O-D demands $Q=q$ in multinomial distribution, $\Sigma_{f|q}=\Sigma_{f|Q=q}$. The matrix $\Sigma_{f|Q}$ is built upon $p$, which can be found in Ma and Qian [49].

Proposition 2**.**

The marginal distribution of $X$ and $X_{m}$ follows,

[TABLE]

[TABLE]

where $x=\Delta f$, $\Sigma_{x}=\Delta\Sigma_{f}\Delta^{T}$.

3.4 Review of existing probabilistic ODE formulations

In this subsection, we will review existing probabilistic ODE formulations and discuss the research gap. The maximum likelihood formulations proposed by Vardi [82], Hazelton [32] rely on the assumption that the O-D demand follows Poisson distributions that are independent among O-D pairs. The Poisson distributed assumption is actually quite restrictive since it explicitly adds a constraint between the mean and variance, often violated in practice [85]. More importantly, the maximum likelihood formulation cannot be easily extended to the case of multivariate probability distributions for O-D demand, since the likelihood function will become intractable.

Shao et al. [68, 69] proposed a generalized least square formulation to estimate the mean and variance of O-D demand that follows the MVN. The formulation is a bi-level optimization problem. The upper level minimizes the difference between observed and estimated link flow mean/variance. The objective function can be written as,

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ are constant weights and $f_{1}(\cdot,\cdot),f_{2}(\cdot,\cdot)$ measures the difference between the observed and estimated link flow mean and covariance/variance, respectively. $x(q)$ implies that the mean of link flow is a function of O-D demand mean, whereas $\Sigma_{x}(\Sigma_{q})$ implies that the variance/covariance matrix of link flow is a function of O-D demand variance/covariance matrix.

As discussed in the introduction section, Formulation 10 does not consider the route choice stochasticity (namely the day-to-day route choice variation) as a result of non-atomic equilibrium. Thus, $\Sigma_{x}$ is only dependent on the O-D variation and physical properties of the network. The derivative of the objective function can be derived analytically. Formulation 10 can be solved using gradient methods. However, if the route choice stochasticity is considered in link flow for atomic users, then $\Sigma_{x}$ is dependent on both O-D demand mean $q$ and variance $\Sigma_{q}$ by Propositions 1 and 2, the objective function becomes,

[TABLE]

where both $x(q)$ and $\Sigma_{x}(q,\Sigma_{q})$ are determined by GESTA [49]. Since the function $\Sigma_{x}(q,\Sigma_{q})$ is complex, deriving its gradient is challenging. Though gradient-free algorithms can be used to find the optimal solution, those methods can have a hard time being applied to large-scale networks. To effectively solve the new objective function subject to the statistical equilibrium based GESTA, this paper proposes the following iterative method.

3.5 A novel iterative formulation for estimating probabilistic O-D demand

Instead of estimating the mean and variance/covariance matrix of O-D demand simultaneously, we propose an Iterative Generalized Least Square (IGLS) framework to estimate mean and variance/covariance iteratively [29]. A large number of statistical algorithms are developed for IGLS. Del Pino [24] shows that IGLS shares some properties with the Newton-Raphson algorithm, and IGLS can be used for solving Maximum Likelihood Estimator (MLE) and quasi-MLE problems.

IGLS framework mainly contains two sub-problems: estimating the O-D demand mean and estimating the variance/covariance matrix. IGLS runs the two sub-problems iteratively to update estimators $\hat{x},\hat{f},\hat{p},\hat{q},\hat{\Sigma}_{q}$. In the sub-problem of estimating the mean vector $q$, $\hat{\Sigma}_{q}$ and $\hat{\Sigma}_{x}$ are seen as given. In the sub-problem of estimating the variance-covariance matrix $\Sigma_{q}$, $\hat{q}$ is seen as given. As we will show later, the link flow covariance $\Sigma_{x}$ follows the Wishart distribution. We then formulate the ODE as an approximated maximum likelihood estimator of the covariance matrix with Lasso regularization.

Compared to the framework by Shao et al. [68, 69], we argue that the IGLS-based formulation is statistical interpretable and more computationally efficient to solve on large-scale networks. As shown in Del Pino [24], each iteration in the IGLS resembles one gradient descent step in the Newton-Raphson algorithm. Since the convergence rate of the Newton-Raphson algorithm is quadratic, a good convergence rate can be expected for IGLS. Later, we also show that the IGLS framework can theoretically provide a better performance in terms of the solution algorithm.

Before detailed models for each of the two sub-problems are presented, we first discuss the stopping criterion. Stoping criterion is needed for three processes in the IGLS framework: estimating mean, estimating O-D variance/covariance matrix and the overall IGLS iteration. For the two sub-problems, since both optimization problems are convex (as we will show later), the stopping criterion can be easily defined [9]. For the stopping critera for the overall IGLS iteration, note that our ultimate goal is to estimate the probability distribution of the O-D demand $Q$. We keep track of the discrepancy between the mean and variance/covariance matrix along the iterations. If the discrepancy is sufficiently small, then the IGLS iteration can be stopped.

Definition 2** (Stopping Criterion for IGLS).**

If the probability distribution of $Q$ is estimated as ${\cal N}(\hat{q}^{+},\hat{\Sigma_{q}}^{+})$ following the estimator from a previous iteration ${\cal N}(q,\Sigma_{q})$, define

[TABLE]

where $D(\cdot,\cdot)$ is a distance measure between the two estimators. The IGLS iteration terminates when $\tau$ is sufficiently small.

There are a large number of candidates for the choice of $D(\cdot,\cdot)$. In this paper we choose two distance functions: Hellinger distance and Kullback-Leibler distance (KL distance), computed as the following for any two MVNs ${\cal N}(\mu_{1},\Sigma_{1}),{\cal N}(\mu_{2},\Sigma_{2})$ with the same dimension $d$,

[TABLE]

3.6 Estimating the mean of O-D demand

In this sub-problem, we estimate O-D demand mean $q$, assuming the estimated variance/covariance of the link flow and O-D demand $\hat{\Sigma}_{x},\hat{\Sigma}_{q}$ is given. The estimated variance/covariance of the link flow on those observed links $\hat{\Sigma}_{x}^{o}$ is also given.

3.6.1 Estimating probabilistic link flow on observed links

Before estimating the O-D demand mean, we first estimate the mean of link flow on observed links, $x^{o}=\mathbb{E}(X^{o})$, as an intermediate step from the observed data ${\mathsf{x}}_{i}^{o}$ to the unknown demand mean $q$. If the unknown errors $\Sigma_{e}$ can be calibrated exogenously, the probability distribution of $X_{m}$ can be determined by ${\cal N}(x,\Sigma_{x}+\Sigma_{e})$ given the probability distribution of $X$. In most cases, $\Sigma_{e}$ is never known, and it represents the errors that cannot be explained by the model. We then estimate $X=X_{m}$ with the hope that $\Sigma_{e}=\mbox{\bf 0}$. After the estimation process, the mismatch between $\hat{X}$ and the data can be viewed as the unknown errors $E$. This will be discussed later.

The reason why we start from estimating $x^{o}$ is that in the classical deterministic ODE, $\left\lVert\Delta pq-\hat{x}^{o}\right\rVert_{2}^{2}$ is used as the objective function, where $\hat{x}^{o}$ represents the daily average of traffic counts. To interpret and understand the probabilistic ODE, we need to rigorously estimate $\hat{x}^{o}$ before constructing an objective function.

Given $\hat{\Sigma}_{x}^{o}$ to approximate $\Sigma_{x}^{o}$, the likelihood function of observed link flow ${\mathsf{x}}_{i}^{o}$ can be constructed as,

[TABLE]

A Maximum Likelihood Estimator (MLE) is to estimate $x^{o}$ to maximize the likelihood,

[TABLE]

Proposition 3** (MLE).**

Given the data set ${\mathsf{x}}^{o}$ and link flow covariance matrix $\Sigma_{x}^{o}$ on observed links, the estimator is $\hat{x}^{o}=\frac{1}{n}{\mathsf{x}}_{i}^{o}$.

The detailed analysis of the optimization problem can be found in A. The estimator $\hat{x}^{o}$ can be derived in a closed form, $\hat{x}^{o}=\frac{1}{n}{\mathsf{x}}_{i}^{o}$. The proof can also be found in A. It is no surprise that the estimator of $x^{o}$ is the average of observed data ${\mathsf{x}}_{i}^{o}$. This is consistent with the formulation of deterministic ODE that searches the best $\hat{q}$ to minimize the discrepancy between estimated link flow derived from $\hat{q}$ and daily average observed counts (namely $\frac{1}{n}{\mathsf{x}}_{i}^{o}$), both for those observed links.

Having the closed formulation of $\hat{x}^{o}$, we can further derive the probability distribution of $\hat{x}^{o}$.

Proposition 4**.**

If the observed link flow ${\mathsf{x}}_{i}^{o}$ i.i.d follows the probability distribution of ${\cal N}(x^{o},\Sigma_{x}^{o})$, then $\hat{x}^{o}$ follows:

[TABLE]

3.6.2 Estimating demand mean on an uncongested network

For an uncongested network, the travel costs are dependent on free-flow speed and road length, and independent of link/path flow. The route choice probability $\tilde{p}$ can be calculated exogenously. Given the route choice probability is known, we estimate the demand mean $q$ by,

[TABLE]

Formulation 18 can be viewed as a generalized least square (GLS) formulation when regarding $n\left(\Sigma_{x}^{o}\right)^{-1}$ as the weight matrix in [67, 10], which minimizes the weighted discrepancy of the mean link flow on observed links between from data $\hat{x}^{o}$ and from the $q$ estimator, $x=\Delta\tilde{p}q$. The above formulation can also be statistically interpreted as a maximum likelihood estimator (MLE) of the O-D demand $q$, provided with the probability distribution of the mean link flow estimator $\hat{x}^{o}$.

A deterministic ODE often uses Formulation 19 where the diagonal matrix $\Theta_{x}$ denotes the confident level for each of the observed link flow.

[TABLE]

This formulation is similar to the Formulation 18 if the row vectors of $\Delta^{o}$ are fully ranked, except for that the weights on each link flow observation can differ. In practice, when $\Sigma_{x}^{o}$ cannot be derived directly (e.g., due to insufficient data points), we can use the simplified Formulation 19 to estimate O-D demand mean $q$.

Another issue for both Formulations 18 and 19 is that the optimal $\hat{q}$ may not be unique [93]. To address the non-uniqueness issue, history O-D information is usually employed. An extended generalized least square can be built by assuming that history O-D information can be acquired and is independent of observed flow given $q$, then we have following formulation.

If the historical O-D demand mean and covariance matrix is given as $q^{H},\Sigma_{q}^{H}$, respectively, we identify the unique solution $q$ by,

[TABLE]

where the historical O-D covariance matrix $\Sigma_{q}^{H}$ is usually unknown. The identity matrix is used as an alternative choice.

Cascetta [10] proposed a similar formulation to Formulation 20, but the derivation of the inverse of observed link flow covariance matrix ${\Sigma_{x}^{o}}$ is unclear. In our formulation, ${(\Sigma_{x}^{o}})^{-1}$ can be derived analytically given $\Sigma_{q}$ using GESTA.

Based on Formulation 20, the number of observed data $n$ and the quality of $q^{H}$ are two major factors affecting the accuracy of estimated O-D demand $\hat{q}$. Several remarks are made regarding $n$ and $q^{H}$.

Remark 1**.**

Incrementing data quantity does not address the non-uniqueness issue of Formulation 20.

If $\Delta^{o}$ are fully ranked, Formulation 18 has a unique solution for any $n\geq 1$. However, when $\Delta^{o}$ is not fully ranked, Formulation 18 has multiple optimal solutions regardless the value of $n$.

Remark 2**.**

Incrementing data quantity increases the accuracy of estimated O-D demand $\hat{q}$.

To simplify the discussion, we assume $\Delta^{o}$ is fully ranked, namely $\Delta^{o}\tilde{p}$ is invertible. If no observation is obtained, $n=0$, then the best estimation of $q$ would be $q^{H}$. This implies that the accuracy of $\hat{q}$ is solely dependent on the accuracy of $q^{H}$. When observations are available, the estimated O-D demand $\hat{q}$ is the weighted average of $\left(\Delta^{o}\tilde{p}\right)^{-1}\hat{x^{o}}$ and $q^{H}$. As $n$ increases, $\hat{q}$ becomes close to $\left(\Delta^{o}\tilde{p}\right)^{-1}\hat{x^{o}}$, otherwise to $q^{H}$. If we have infinite data, then $\hat{q}=\left(\Delta^{o}\tilde{p}\right)^{-1}\hat{x^{o}}$. Thus, $\mbox{\rm Var}(\hat{q})=\mbox{\bf 0}$ due to the central limit theorem (CLT) when $n=\infty$, implying $\hat{q}$ is perfectly accurate. In summary, incrementing data quantity reduces $\mbox{\rm Var}(\hat{q})$, from relying solely on the accuracy of $q^{H}$ ($n=0$) to being perfectly accurate ($n\to\infty$).

Remark 3**.**

Formulation 20 is not a Maximum a Posterior (MAP) estimator.

Different from Menon et al. [51], Formulation 20 cannot be interpreted as a Maximum a Posterior (MAP) estimator. We first build an MLE according to the probability distribution of $\hat{x^{o}}$, so we need the prior of $\hat{x^{o}}$ to build the MAP estimaor. However, the information used in Formulation 20 is the historical O-D $q^{H}$ rather than historical link flow $x^{H}$. Thus, formulation 20 is not an MAP estimator. Formulation 20 can only be interpreted as a GLS model when the history O-D is independent of observed traffic flow given $q$. Similar arguments can also be found in Yang and Fan [95].

Remark 4**.**

If $q^{H}\neq q$, the estimator of O-D demand mean $\hat{q}$ from Formulation 20 is biased.

To summarize, we conclude that data quantity helps increase the accuracy of the estimated O-D demand mean $\hat{q}$ while historical O-D information $q^{H}$ addresses the non-uniqueness issue. One subtle issue is that even if we observe a large number of data (on a large number of days), uniquely determining $\hat{q}$ may still be impossible.

3.6.3 Estimating demand mean on a congested network

In congested networks, the route choice probability $\tilde{p}$ is endogenously determined by GESTA. The probabilistic ODE needs to estimate the demand mean $q$ and the route choice probability $p$ simultaneously. Instead of estimating both separately, we can estimate the path flow $f=\tilde{p}q$, analogous to the Path Flow Estimator (PFE) in the deterministic ODE settings.

Because the dimension of path flow $f$ is greater than the O-D demand $q$, using historical O-D information does not necessarily address the non-uniqueness issue in this case. In this subsection, we assume the dimension of $f$ is greater than that of link flow on observed locations $x^{o}$. If it is not the case, then we can always enlarge the path set by generating more paths or shrinking the observation size with network consolidation. The basic formulation is proposed in Formulation 21,

[TABLE]

The fundamental problem resulting the non-uniqueness issue is that the number of observed links is far smaller than the total number of paths. Even though all links are covered by surveillance, the path flow estimator $\hat{f}$ can still be non-unique. However, it is possible to restrict the feasible set of path flow $f$ by route choice models (namely equilibrium conditions), such as GESTA. The restricted formulation is presented in 22.

[TABLE]

Where $\Phi^{+}$ is the feasible set of $f$. We adopt GESTA to model the traffic conditions and route choices since it generally works with any specific route choice models. Here we demonstrate the idea using deterministic UE and Logit/Probit-based SUE as the route choice model.

1. i)

UE-based GESTA

UE can be formulated as an optimization program to minimize $Z_{1}(f)$ [72],

[TABLE]

Then $\Phi_{\text{UE}}^{+}$ is defined as:

[TABLE]

$\Phi_{\text{UE}}^{+}$ can also be written as a link-based or path-based variational inequality formulation [75]. Formulation 22 under UE constraints is known as Mathematical Programming with Equilibrium Constrain (MPEC) [47]. 2. ii)

Logit-based GESTA

In the Logit-based GESTA, travelers’ perception of travel costs is assumed to follow Gumbel distribution. The variance/covariance of the path cost $C$ is not considered. Fisk [27], Janson [37] cast the Logit model to its dual form, a convex optimization problem that minimizes the following objective function:

[TABLE]

Then $\Phi^{+}_{\text{Logit}}$ is defined as:

[TABLE] 3. iii)

Probit-based GESTA

In the Probit-based GESTA, travelers’ perception errors follow Normal distribution [22], as part of the probability distribution of path cost $C$. There does not exist a explicit form on $\Phi^{+}_{\text{Probit}}$. Any pair of $(p,f)$ satisfying the Probit route choice model are in $\Phi^{+}_{\text{Probit}}$. The details of Probit choice model can be found in Daganzo et al. [22], Sheffi [72].

Formulation 22 is also known as the bi-level formulation of ODE with many existing studies [58, 28, 94, 92]. Since the solution uniqueness for the bi-level formulations varies by route choice models, we discuss them in Section 4. Other properties of Formulation 22 are discussed in the following remarks.

Remark 5**.**

Formulation 22 is non-convex.

Since $\Phi^{+}$ is clearly not a convex set regardless of the route choice models, Formulation 22 is not convex. A sensitivity-based algorithm by Josefsson and Patriksson [39] and a heuristic algorithm by Yang [92] are commonly used to solve for it. In addition, Nie and Zhang [59] relaxes the UE-based ODE to a one-level optimization problem, enhanced by Shen and Wynter [73] with a convex relaxation program on a one-level optimization problem.

Remark 6**.**

Logit-based SUE can be approximated by a specific Probit-based SUE.

In the Logit-based SUE, the only parameter for the route choice model is the dispersion factor $\Theta$. Given $\Theta$, a Probit model with a diagonal path cost variance matrix can be used to approximate the Logit model. Details regarding the transformation can be found in Greene [30].

3.7 Estimating the variance/covariance matrix of O-D demand

To estimate O-D demand variance and covariance matrix, we assume the link flow mean $x$ and path flow mean $f$ are provided and known. Consequently, the route choice probability $p$ and O-D demand mean $q$ are also known. We find an MLE to estimate the $\Sigma_{q}$. We first present the basic formulation to estimate $\Sigma_{q}$ given that $\Sigma_{x}$ follows the Wishart distribution [90]. Due to the high dimension of $\Sigma_{q}$, Lasso [78] regularized formulation is proposed to search for a sparse estimation of $\Sigma_{q}$ that makes trade off between variance and bias.

3.7.1 Basic formulation

First we define the empirical covariance matrix of the observed link flow to be $S_{x}^{o}=\frac{1}{n}\sum_{i=1}^{n}({\mathsf{x}}_{i}^{o}-\bar{{\mathsf{x}}}^{o})({\mathsf{x}}_{i}^{o}-\bar{{\mathsf{x}}}^{o})^{T}$, which is the maximum likelihood estimator of covariance matrix. $\bar{{\mathsf{x}}}^{o}$ is the averaged observed link flow, $\bar{{\mathsf{x}}}^{o}=\frac{1}{n}\sum_{i=1}^{n}{\mathsf{x}}_{i}^{o}$. Note $S_{x}^{o}$ is different from the sample covariance matrix $P_{x}^{o}=\frac{1}{n-1}\sum_{i=1}^{n}({\mathsf{x}}_{i}^{o}-\bar{{\mathsf{x}}}^{o})({\mathsf{x}}_{i}^{o}-\bar{{\mathsf{x}}}^{o})^{T}$. Since the link flow variance-covariance matrix follows the Wishart distribution, a maximum likelihood estimator can be built to solve for $\Sigma_{q}$.

Proposition 5**.**

Given the variance-covariate matrix for observed link flow $S_{x}^{o}$, the maximum likelihood estimator of $\Sigma_{q}$ is

[TABLE]

The first constraint in Formulation 27 is obtained from GESTA through Propositions 1 and 2. The convexity of Formulation 27 depends on the rank of $\Delta^{o}$. If $\Delta^{o}$ is fully ranked, Formulation 27 is non-convex. But if $\Delta^{o}$ is not fully ranked, then we can first find the optimal $\Sigma_{x}$ for the objective function, and then solve for $\Sigma_{x}=\Delta\Sigma_{f|q}\Delta^{T}+\Delta\tilde{p}\Sigma_{q}\tilde{p}^{T}\Delta^{T}$. Both steps are convex optimization problems. In addition, $\Sigma_{q}$ contains $\frac{1}{2}|K_{q}|(|K_{q}|-1)$ elements, which is usually in a higher dimension than the number of observed data, so the optimal estimator of $\Sigma_{q}$ may not be unique.

Next we introduce the regularization and relaxation of Formulation 27 to achieve convexity and uniqueness.

3.7.2 Sparse model selection

Since the number of entries in the O-D variance-covariance matrix is usually much greater than the size of data, a Lasso penalization is used to select the O-D variance-covariance matrix as in Formulation 28,

[TABLE]

$\lambda$ is a Lasso parameter to adjust the sparsity of $\hat{\Sigma}_{q}$. $\Sigma_{f|q}$ is constructed using $p$ and $\hat{q}$. We note Formulation 28 obtains a biased but robust estimator of OD variance/covariance matrix. Formulation 28 is hard to solve due to its non-convexity [8]. Although non-linear optimization methods can be employed to solve this formulation, none of them can guarantee computationally efficiency thus not suitable for large-scale networks. We would prefer to approximate it using a convex optimization problem with Lasso regularization. Inspired by [98], a second order approximation to the MLE of the covariance matrix is used in Formulation 29.

[TABLE]

where $\left\lVert A\right\rVert_{F}=\sqrt{\text{Tr}(A^{T}A)}$ and $\left\lVert A\right\rVert_{1}=\sum_{ij}|A_{ij}|$. The former one is known as Frobenius Norm, equivalent to the element-wise $L_{2}$ norm [38]. The latter one is the element-wise $L_{1}$ norm.

Proposition 6** (Convexity).**

The optimization problem 29 is convex.

Proof.

$\Sigma_{q}$ can only be positive semi-definite matrix, which forms a convex set. Then we plug $\Sigma_{x}=\Delta\Sigma_{f|q}\Delta^{T}+\Delta\tilde{p}\Sigma_{q}\tilde{p}^{T}\Delta^{T}$ into the objective function. The objective function with respect to $\Sigma_{q}$ is also convex, so the entire formulation is convex. ∎

Formulation 29 is a desired optimization problem to solve for a sparse O-D demand variance/covariance matrix $\Sigma_{q}$, because the convexity allows its computational efficiency. In principle, we can use proximal methods [61] to solve Formulation 29. Details of the solution algorithms can be found in B.

3.8 Incorporating day-to-day travel time data in probabilistic ODE

The day-to-day travel time/speed data can be added to the Formulation 22 to further enhance the ODE. Real-time traffic speed data vendotrs, such as INRIX and HERE, can provide traffic speed data covering major roads in most of U.S. cities. Some studies [4, 48, 40] regarded the observed travel time/speed as another objective to minimize, and thus enhance the Formulation 22 to become Formulation 30.

[TABLE]

where $t(\cdot)$ is the link performance function that maps the link flow to link costs (such as the well known BPR functions). $w_{1}$, $w_{2}$ and $w_{3}$ are weights assigned to each objective.

Similar to estimating the covariance matrix of link flow by Formulation 29, one can also estimate the covariance matrix of travel cost/time given an estimator for its mean. However, a bigger issue is that the mapping from traffic speed to traffic hourly volume on the road segment is not a one-to-one mapping. A link performance function, though uniquely maps travel cost/time to volume in both ways, can be very sensitive in determining volume given near free-flow travel time. When travel speed/time data based on probe vehicles is highly biased, the error can be amplified through the link performance function. Therefore, the ODE relying on travel speed data in the static network settings is practically challenging. We believe that applying travel speed/time data can be more useful when extending GESTA and probabilistic ODE to dynamic network settings where the traffic dynamics is captured using microscopic or mesoscopic flow models. We hope to address probabilistic dynamic ODE in a future research paper.

4 Some properties of the formulations

In this section, we discuss how to evaluate the accuracy and effectiveness of our estimated O-D demand mean and variance/covariance matrix. We analyze variance/covariance by its decomposition into three main sources, O-D demand variance, route choice variance and unknown (unexplained) error. The observability of the proposed probabilistic ODE framework is also discussed.

4.1 Goodness of fit

For classical ODE, the goodness of fit indicator measures how close the estimated O-D demand mean $\hat{q}$ can, if loaded into the network following a deterministic traffic assignment model (UE or SUE), reproduce the observed traffic conditions. Commonly used indicators are summarized in Antoniou et al. [2]. Similarly, for the proposed probabilistic ODE, the goodness of fit can be measured by how close the estimated probabilistic O-D demand $\hat{Q}$, if loaded into the network following GESTA, can reproduce the probability distribution of observed traffic flow. Define the estimated link flow on observed links from the probabilistic ODE $\hat{X}^{o}\sim{\cal N}(\hat{x}^{o},\hat{\Sigma}_{x}^{o})$ and estimated link flow on observed links directly from data ${\mathsf{X}}^{o}\sim{\cal N}(\bar{{\mathsf{x}}}^{o},P_{x}^{o})$, $P_{x}^{o}=\frac{1}{n-1}\sum_{i=1}^{n}({\mathsf{x}}_{i}^{o}-\bar{{\mathsf{x}}}^{o})({\mathsf{x}}_{i}^{o}-\bar{{\mathsf{x}}}^{o})^{T}$. The goodness of fit indicator can be computed by the Hellinger distance or Kullback-Leibler distance between $\hat{X^{o}}$ and ${\mathsf{X}}^{o}$.

4.2 Variance analysis

In Section 3, we do not consider $\Sigma_{e}$. With the real world data, the observed data cannot be fully explained by the ODE, and thus contain unexplained errors. After the probabilistic ODE process, we can analytically decompose the link flow variance to check how much variance can be explained by the ODE.

Proposition 7** (Link flow variance decomposition).**

The variance of link flow can be decomposed into three parts, O-D demand variance, route choice variance, and unknown errors.

[TABLE]

There are many ways to quantify the variance ratio. In this study, we determine the ratio based on matrix traces, which is widely adopted in the statistics literature. Trace-based variance ratio is closely related to the spectral analysis of recurrent link/path flow data. Details and examples can be found in Ma and Qian [49].

4.3 Observability

ODE is notoriously difficult because it is underdetermined. Studies on O-D observability problem specifically discusses the issue of solution non-uniqueness [74, 95]. In this subsection, we discuss the uniqueness property of the proposed probabilistic ODE.

As we discussed in Section 3, under no congestion, Formulation 20 is able to estimate O-D mean $q$ uniquely once prior information $q^{H}$ is introduced [6, 93]. Studies also suggested estimate $\hat{q}$ by taking the pseudo-inverse of $\Delta^{o}$ that encodes a singular value decomposition (SVD) process in its formulation [60]. For a congested network with the Formulation 22, its solution $\hat{f}$ may not be unique. The observability of Formulation 22 varies by the constraints $\Phi^{+}$, dependent on the specific route choice model adopted under GESTA.

1. i)

UE-based GESTA

Generally path flow under the UE condition is not unique given O-D demand $q$ [75]. When UE is used as the constraint for Formulation 22, it cannot guarantee the optimal path flow estimator $\hat{f}$ to be unique. However, we can find an extreme point solution from the feasible domain, and then this solution to Formulation 22 is unique [80]. The extreme point solution can be obtained through column generation. The optimal estimator $\hat{f}$ is also locally stable and the upper level object function is strongly convex [80, 63]. Thus, if we use history O-D demand as the initial point and the history O-D is near the true O-D demand, then the optimization process is likely to find the optimal solution without trapping into a local minimum [92]. 2. ii)

Logit-based GESTA

Since the Logit model is strictly convex on $f$, the optimal solution to Formulation 22 is unique [63]. 3. iii)

Probit-based GESTA

Since $\hat{\Sigma}_{x}$ is given in Formulation 22, the probability distribution of path costs $C$ is uniquely determined. As a result, the solution to the route choice probability $p$ is unique under the Probit model. If the history O-D information is used, the optimal estimator of O-D demand mean $\hat{q}$ is unique. Since both $p$ and $\hat{q}$ are unique, the optimal estimator of path flow $\hat{f}$ is also unique.

As for Formulation 29 to estimate O-D demand variance/covariance matrix, the optimal solution is non-unique since Lasso regularization is not strictly convex [79]. However, practically, Lasso regularization can largely shrink the solution domain towards being unique.

Proposition 8**.**

The optimal solution $(\hat{q},\hat{\Sigma}_{q},\hat{f},\hat{x})$ to the IGLS framework consisting of both formulations 22, 29 may not be unique for any given observed link flow data set ${\mathsf{x}}^{o}$.

The major reason for the non-uniqueness is that the number of rows (namely the number of links that are covered with sensors) in $\Delta$ or $\Delta^{o}$ is much smaller than the number of its columns (namely the number of paths) in a general large-scale network. Consequently, $\Delta^{o}Pq$ is unique, but $q$ is not unique in Formulation 18. $\Delta^{o}f$ is unique, but $f$ is not unique in Formulation 21. Similarly, $\Sigma_{x}$ is unique, but $\Sigma_{q}$ is not unique in Formulation 27. Generally, the entire IGLS framework estimates both the mean and variance/covariance matrix, and thus has to search a much larger domain space than a deterministic ODE. Therefore, its observability is worse off. This will be further demonstrated in the numerical experiments.

Though Proposition 8 declares a challenge for estimating probabilistic O-D demand, we argue that by the proposed IGLS framework, the O-D mean estimator is no worse than a best possible estimator by an error that reduces with respect to the sample size, and thus no worse than the O-D demand estimator using deterministic ODE methods.

Proposition 9**.**

Suppose observations of link flow on observed links that are i.i.d drawn from the probability distribution of $X$ on each day, and they are used to estimate the mean and variance/covariance matrix of link flow. $\hat{\Sigma}_{x}\succ 0$. For an arbitrary route choice probability vector $p\geq 0$ (or equivalently an underlying route choice model), the statistical risk of the estimated O-D mean $\hat{q}$ from Formulation 18 (or Formulation 22 without history O-D information) is of ${\cal O}\left(\frac{1}{n}\right)$ where $n$ is the sample size (namely the number of days with observations).

Proof sketch.

Note that the observed link set $A^{o}$ does not change from day to day. First we define a risk function to measure the performance of a specific estimator. The risk is low when the estimator provides an accurate solution given any observed link flow data from a probability distribution of $X$, whereas the risk is high when the estimator is either inaccurate or not robust to the observed link flow. We then rewrite the ODE formulation and bound the risk. We show that the risk of the estimator for O-D demand mean is of ${\cal O}(\frac{1}{n})$ regardless of the quality of the estimated link variance/covariance matrix $\Sigma_{x}$. A detailed proof is provided in C. ∎

Proposition 9 is one of the major features that distinguish this research from other existing probabilistic ODE methods. Though probabilistic ODE works with a much larger solution space than the deterministic O-D ODE, Proposition 9 guarantees that using the proposed IGLS framework, the estimator for O-D demand mean is no worse than the best possible estimator by ${\cal O}(\frac{1}{n})$, and thus the mean estimated by deterministic ODE. Provided with a large data sample, the proposed probabilistic ODE will not “get lost” due to enlarged searching solution domain regardless of the variance/covariance matrix. In other words, the estimator for O-D variance/covariance matrix can be seen as additional information to be inferred using day-to-day traffic data, in addition to the mean estimator. The variance/covariance matrix does not impair the performance of estimated O-D mean vector. This is one critical feature that distinguishes our research from Shao et al. [68], where the formulation solving for both mean and covariance matrix simultaneously may not necessarily guarantee a robust estimator for the O-D demand mean.

5 Solution algorithms

In this section we present the solution algorithm for the proposed IGLS framework. The goal is to compute the estimators for O-D mean and variance/covariance matrix $(\hat{q},\hat{\Sigma}_{q})$. The proposed formulations are path based. The number of paths with positive flow increases exponentially when the network grows. For small networks, path enumeration is possible. When the networks are large, we can simply enumerate $K$ shortest paths [97, 26] for each O-D pair and then search for the solution in the prescribed path set. In addition, the proposed IGLS framework can also fit the column generation method [89, 66]. At each iteration, one or several additional paths that possess minimal path cost at the time of iteration can be generated and added to the prescribed path set.

For the sub-problem of estimating O-D demand mean vector $q$, two heuristic algorithms can be used to directly solve the bi-level formulation [92, 39]. A single-level convex relaxation to the formulation can also be adopted [73]. In addition, two algorithms, Iterative Shrinkage-Thresholding Algorithm (ISTA) [56] and Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [57] solve for the sub-problem of estimating O-D demand variance-covariance matrix $\Sigma_{q}$.

The solution algorithm is summarized as follows,

[TABLE]

6 Numerical experiments

We first examine the probabilistic ODE method on two small networks. Results are presented, discussed, and compared among different route choice models. The sensitivity of the data quantity and historical O-D information are also tested and analysed. The impact of penalty term on Lasso regulrazation is analysed. We also compare the efficiency of different gradient-based methods to solve Formulation 29. In addition, the proposed method is also applied on two large-scale real-world networks to examine its efficiency and scalability.

The true O-D demand in the real world is notoriously difficult to obtain. We adopt two methods to validate the proposed probabilistic ODE. In the first method, we construct probabilistic O-D demands then regard it as the “true” O-D demand. We run GESTA with Probit as the route choice model to obtain the probability vector $p$, again as the “true” choice probabilities. Then we randomly sample a set of O-D demand from its probability distribution, further randomly sample the route choice for each of the trips to obtain link/path flow, and finally add a perturbation of error (as the unknown error) to the link flow. The perturbation of error from $-\varepsilon$ to $\varepsilon$ is generated as follows: the perturbed value is $\xi_{p}=\xi(1+\texttt{rand})\varepsilon$ when the actual value is $\xi$, where rand is a sample uniformly distributed between $[-1,1]$. We do this random sampling for a sequence of many trials, each of which is seen as the observation on one day. A subset of the link flow is used as the observations. The performance of the probabilistic ODE is assessed by comparing the estimated O-D demand to the “true” O-D demand [2]. In the second method, we estimate the probabilistic O-D demand using real-world day-to-day traffic flow count data using this IGLS framework, and check if the estimated O-D demand, along with GESTA, can reproduce the actual traffic flow observed for a set of days.

To measure the error of estimated O-D mean, we use Percentage Root Mean Square Error (PRMSE). Kullback-Leibler distance (KL distance) is used to measure the error of estimated O-D probability distribution.

[TABLE]

In the numerical experiments, we test a few different settings for the sub-problem of estimating O-D demand mean, as well as for the sub-problem of estimating variance/covariance. “w/o EC” implies that Formulation 21 is adopted without an equilibrium constraint. We use “Logit” and “Probit” to denote Logit-based GESTA and Probit-based GESTA, respectively. When estimating O-D demand variance, “w/o Lasso” represents the setting without the Lasso regularization and “w/ Lasso” with the Lasso regularization.

6.1 A small three-link network

We first work with a toy network with three links, three paths and two O-D pairs as shown in Figure 3. The Bureau of Public Roads (BPR) link travel time function is adopted,

[TABLE]

where $t^{0}_{a}$ is the free-flow travel time on link $a\in A$, $\beta=4$, $\alpha=0.15$ are constant parameters. $cap_{a}$ denotes the capacity of link $a$. Link settings are $t^{0}_{1}=t^{0}_{2}=10,t^{0}_{3}=5$, $cap_{a}=360,\forall a\in A$. OD pairs $(1\to 3)$ and $(2\to 3)$ are considered, demand means are $q_{1\to 3}=700,q_{2\to 3}=500$. The variance of the O-D demand is set to be $25\%$ of the O-D demand mean. Probit-based GESTA is used as the “true” underlying statistical traffic assignment model.

6.1.1 Estimation results

Suppose we observe $500$ days’ traffic counts on Link $1$ and Link $3$. Those observations are used to estimate the probabilistic O-D demand by different modeling settings. We tested three demand patterns where the “true” synthesized correlation of the demand between the two O-D pairs $\rho$ is $-0.5,0,0.5$ respectively. For each correlation, the observed data is synthesized and presented in Figure 4. The ODE would not know the true underlying Probit route choice model. Instead, ODE would “speculate” a specific route choice model, Logit, Probit or no route choice model, to examine their respective performance. Furthermore, the historical O-D information is unknown. Given traffic observations of 500 days and a speculated route choice model, we apply the proposed probabilistic ODE method to estimate the probability distribution of O-D demand, all results are presented in Table 2. Note that when we apply the ODE with the Logit model, we identify the dispersion factor $\Theta=1$ such that it produces the best estimation performance of all possible values for $\Theta$. In fact, the value of $\Theta$ only marginally impacts the ODE performance, so the choice of its value does not affect our findings as much.

First of all, we obtain not only the O-D mean vector but also the O-D variance/covariance (namely the correlation in this example) by applying the probabilistic ODE to make the best use of all 500 days’ data. As can be seen in Figure 4, the daily average flow counts are the same for different $\rho$. If we apply the deterministic ODE, then the 500 days’ data are taken the average to estimate the O-D demand mean. Day-to-day traffic data are not fully used, and the correlations of O-D demand among O-D pairs are overlooked. The probabilistic ODE provides more insights of O-D demand that would be needed for both transportation planning and operation.

When comparing different settings of the probabilistic ODE, all those settings yield acceptable estimations. The RMPSEs of the probabilistic ODE, under different true demand patterns, route choice models, and whether equilibrium constraints are considered in PFE, are all below $4\%$. In general, the probabilistic ODE with different settings accurately estimates not only the mean but also the variance and covariance of O-D demand.

Since the “true” link flow is generated based on the Probit model, it is no surprise that the ODE with the Probit model yields the best performance among all route choice models, consistently for all $\rho$ values. The ODE with the Logit model also provides a good estimate, close enough to the ODE with the Probit model, as long as the dispersion factor $\Theta$ is set properly. However, the Logit based ODE is biased since it only approximates, but cannot fully capture the true demand correlation between the two O-D pairs. If no equilibrium constraints are used, the accuracy of estimate O-D demand can go both ways. This is also no surprise as a result of a much larger domain space for the path flow compared to Logit-based or Probit-based GESTA. For this small network, it outperforms the Logit model when $\rho=-0.5$ but is less accurate when $\rho=0.5$.

Lasso regularization leads to a more accurate estimation when the true demand is independent between the two O-D pairs. The sparse O-D variance/covariance can be better interpreted with more insights, and may furthermore allow causal inference/analysis of trips made among traffic analysis zones. Lasso regularization will, however, make the estimator biased, a disadvantage of using sparse variance/covariance matrix. This is why the ODE performance with Lasso when the correlation is 0.5 or -0.5 is worse than the ODE without Lasso. If all the variance and covariance are substantial, using Lasso leads to a substantial bias for the estimator. Note that variance and bias is always a trade-off to make in the probabilistic ODE. When the weight parameter $\lambda$ for Lasso is carefully chosen (for example by cross validation), the bias may be small in exchange for a much more reliable estimator (namely a much smaller variance) comparing to the settings without Lasso regularization.

6.1.2 Variance decomposition

After obtaining the estimated probabilistic O-D demand, we can conduct variance analysis for the flow on each link. As an example to demonstrate the variance decomposition, we consider the estimated probabilistic O-D demand using observations drawn from $\rho=0.5$ and the Probit model as the route choice model (listed in Table 2). The decomposition of link flow variance is presented in Figure 5. Most of the day-to-day flow variance on Link $2$ comes from the route choice variation. In other words, change in demand levels does not affect the link flow on Link 2 as much. This is because Link $1$ has lower cost than Links 2 and 3 combined. Demand from node $1$ will prefer using Link $1$. As a result, the flow variance ratio attributed to O-D demand on Link $2$, $p_{2}^{2}q_{1\to 3}$, is low (where $p_{2}$ is the probability of choosing Link $2$ for demand from 1 to 3), while the flow variance ratio attributed to route choices, $p_{1}p_{2}q_{1\to 3}$, is much higher (where $p_{1}$ is the probability of choosing Link $1$ for demand from 1 to 3). Link 1 is not used by the demand from node 2 to 3, and thus this demand and its route choice have an indirect impact on the flow variation on Link 1. Comparing to Link 1, Link 3 is directly affected by the demand of both O-D pairs. Consequently, Link $3$ has the highest day-to-day flow variance of all links.

6.1.3 Sensitivity analysis

So far, we have not used any historical O-D demand information in the probabilistic ODE. Historic O-D information may have substantial impact to the performance of the probabilistic ODE. We now examine how sensitive the ODE results are with respect to the quality of historical O-D demand mean $q^{H}$. We change the historic O-D demand from $70\%$ to $130\%$ of the “true” synthetic O-D demand. Other settings are the same as in previous experiments. For three ODE settings (no route choice, Probit and Logit), the KL distance from the estimated O-D demand to the “true” O-D demand are shown in Figure 6 with respect to the quality of historical O-D demand data, represented by how far the historical O-D demand mean provided to the ODE is away from the “true” O-D demand mean. The variance/covariance matrix of the historical O-D demand is set to the identity matrix.

Shown in Figure 6, the probabilistic ODE without equilibrium constraints is the most sensitive to biased history O-D demand mean, as a result of a larger domain space. Its resultant KL distance is almost twice as much as the ODE with Probit-based GESTA, given the same inaccurate historical O-D demand mean. Provided with the accurate historical O-D demand mean, the ODE of both with Probit-based GESTA and without equilibrium constraints can accurately estimate the O-D demand mean. However, this is not the case for ODE with Logit-based GESTA. Generally, the ODE with Logit-based GESTA results a less KL distance than the ODE with Probit-based GESTA when the provided historical O-D demand is less than the “true” demand, and it results a great KL distance when the provided historical O-D demand is greater than the “true” demand. Clearly, the ODE with Logit-based GESTA embeds a prior bias on the demand mean, and tends to estimate more demand than the provided historical demand. Overall, the ODE with Probit-based GESTA looks like the most robust estimator to the historical information in this case, possibly because it happens to use the “true” route choice model, namely the Probit model.

Next, we examine the impact of data sample size to the ODE results. In this experiment, we again do not use the historical O-D information. The ODE without equilibrium constraints does not seem to improve as the sample size increases. This implies that 100 days of observation, in this case, does not necessarily guarantee reliable O-D demand estimate. However, if the GESTA with a route choice model is adopted, then increasing sample size can improve the ODE results. In this case, when the data sample size is more than 20, the ODE results are reasonably good. And when the size reaches 70, the ODE can provide the solution very close to the “true” probability distribution of the O-D demand.

6.1.4 Efficiency of solution algorithms: ISTA versus FISTA

We now examine the performance of solution algorithms, ISTA and FISTA, when solving Formulation 29. Changes in the objective function value against the number of iterations using both algorithms are presented in Figure 8. Clearly, FISTA is more efficient in solving the minimization problem than ISTA. In the following experiments, we use FISTA as the sole solution algorithm.

6.2 A second small network

A second small toy network is used to demonstrate the effects of Lasso penalty term $\lambda$ on the estimation results. This toy network contains $6$ links, $7$ nodes and $5$ O-D pairs, as shown in Figure 9. All 5 O-D pairs share the same link, Link 1. The free-flow travel time $t_{a}^{0}$ is $10$ for each link $a$. The capacity for Link $1$ is 1,800, and the capacity of Links 2, 3, 4, 5 and 6 are randomly drawn from $250$ to $750$. The standard BPR link performance function is adopted. The “true” O-D demand is synthesized by randomly drawing from $300$ to $700$ for each O-D pair. The “true” O-D variance/covariance matrix is set to contain $6$ zero entries out of the total $25$ entries. We synthesize the “true” variance of the O-D demand using its mean, and the “true” correlation factor is randomly drawn from $-0.5$ to $0.5$.

Suppose we have a full coverage on the network, namely, the flow counts on all the links are observed for in all $1000$ days. We estimate the probabilistic O-D demand using the proposed IGLS framework without using historical O-D information. Since we have a full coverage for all links, the estimated mean of the O-D demand is fairly accurate, whereas the estimated variance/covariance matrix is largely dependent on the LASSO penalization term $\lambda$. The relation between values of $15$ variance/covariance entries (due to the symmetry of this matrix) and $\lambda$ is presented in Figure 10, also known as coefficient paths for LASSO. Each path represents the value of one entry in the O-D variance/covariance matrix under different LASSO penalty $\lambda$. As can be seen from Figure 10, when the Lasso penalty is high, three variance/covariance entries are the most significant, implying these three pairs are the most correlated. The results are consistent with the “true” O-D demand. On the other hand, if Lasso penalty $\lambda$ is too small, most of the entries are selected and non-zero, which is inconsistent with the “true” demand that has 6 zero entries. Only when the Lasso penalty $\lambda\in[0.05,0.15]$, the variance/covariance matrix can be estimated accurately.

6.3 A large-scale network: California SR-41 corridor network

The proposed framework is now applied to a real-world network to demonstrate its computational efficiency. The SR-41 corridor network is located in the City of Fresno, California. This network consists of one major freeway and two parallel arterial roads connected with local streets. The network is presented in Figure 11, containing 2,413 links and 7,110 O-D pairs. Its O-D demand mean was calibrated by Liu et al. [45], Zhang et al. [99].

We assume that the “true” O-D demand variance is the same as its mean (similar to a Poisson distribution), and that $10\%$ of O-D pairs (randomly chosen) are mutually correlated with a correlation randomly drawn from $-0.5$ to $0.5$. We randomly choose $50\%$ of the links on the network to be observed for 1,000 days. Again, a standard BPR function is used for all links. We use Logit-based GESTA as the underlying statistical traffic assignment model, paired with Lasso regularization and historical O-D demand mean. The historical O-D demand variance/covariance matrix is set to the identity matrix. As for the historical O-D demand mean, we uniformly sample a perturbation value from -20% to 20%, independently for each O-D pair. It is also used as the initial values for the probabilistic ODE process. The path set is generated by running $3$-shortest paths algorithm for each O-D pair before the ODE process. In all 17,835 paths are considered. The path set for the entire network is assumed to be pre-determined and fixed during this estimation process. The proposed method is developed under MATLAB 2014a and runs on a regular desktop computer (Inter(R) Core i5-4460 3.20 GHz $\times$2, RAM 8 GB). As a result, the average computation time for one iteration of updating both mean and variance/covariance matrix is $301.82$ seconds. The memory usage over first $10$ IGLS iterations is presented in Figure 12. In general, the sub-problem of estimating O-D demand mean consumes less memory than the sub-problem of estimating the variance/covariance matrix. Peak memory usage is around $4.5$GB for this network settings. The memory usage is closely related to the sparsity level of the O-D variance/covariance matrix.

We perform $99$ iterations for the entire IGLS framework. Under each IGLS iteration, we perform $9$ iterations for each of the sub-problems. The convergence of both O-D demand mean and variance/covariance matrix is presented in Figure 13, at the level of iterations for sub-problems. As can be seen, both sub-problems can be solved very efficiently. The entire process of 900 iterations takes $486$ minutes, but the estimate is reasonably good within approximately $300$ minutes. In addition, we plot the estimated path flow mean against “true” flow mean, as well as estimated O-D demand variance/covariance against “true” variance/covariance in Figure 14. Figure 15 plots and the estimated link flow mean against “true” flow mean, as well as the estimated link flow variance, against “true” variance/covariance of the marginal distributions of link flow. The proposed probabilistic ODE seems computationally plausible on a sizable network and is able to achieve reasonably accurate results, approaching the synthesized “true” day-to-day demand mean and covariance.

One interesting result is that the Lasso regularization on the demand variance/covariance matrix unavoidably leads to a biased estimation, as can be seen from both Figure 14 and Figure 15. Most of the variance/covariance (for both O-D demand and link flow) are either estimated as zeros, or substantially underestimated, as a result of Lasso shrinking. However, without the Lasso regularization, most of those variance/covariance entries in the matrix that are substantial can go way off the chart. Again, a proper Lasso penalty ensures a good trade-off between bias and variance of the estimation. In practice, trial-and-error may be needed to identify a proper Lasso penalty.

We also conduct another experiment for this network setting with positively correlated O-D demand, namely 10% of O-D pairs (randomly chosen) are mutually correlated with a correlation randomly drawn from [math] to $0.5$. All other settings are the same as before. This experiment aims at examining the robustness of the proposed probabilistic OD estimator. The estimation results are presented in Figure 16. Clearly the proposed method accurately estimates the probabilistic O-D demand in terms of both the mean and variance-covariance matrix.

6.4 A second large-scale network: Washington D.C. Downtown Area

Previous experiments are conducted in a simulated environment where observations data were synthesized. In this subsection, we apply our probabilistic ODE to a real world network in Washington D.C. by using the actual day-to-day traffic count data.

This network is generally a grid network that consists of $984$ road junctions, 2,585 road segments and 4,900 O-D pairs, The overview of the network is shown in Figure 17. Red dots on the map represent those active fixed-location sensors. There are in all $51$ sensors in the region, while only $10$ of those sensors are working under healthy conditions and collecting data continuously from 2008. We obtain the traffic counts data from August 2008 through December 2015. Around 2,000 data samples were observed for each sensor. We box plot the aggregated traffic counts of all sensors during the morning peak over the years (four of them are presented in Figure 18). The box plot shows that the mean and variance of each sensor do not change as much over the $7$ years. Thus, we decide to use all the days for each sensor to estimate the probabilistic O-D demand representing the demand over a course of 7 years.

We again use Logit-based GESTA as the underlying statistical traffic assignment model, paired with Lasso regularization and historical O-D demand mean. The historical O-D demand variance/covariance matrix is set to the identity matrix. The historical O-D demand mean is obtained from the planning model of Year 2013 developed by Metropolitan Washington Council of Governments, which is also used as the initial demand mean for the solution process. The initial demand variance/covariance matrix is randomly generated.

The convergence of both sub-problems are presented in Figure 19. Overall, the solution algorithm performs well on both sub-problems. It takes $40.18$ minutes to complete all 900 iterations using the same programming environment and aforementioned computer. The estimated O-D demand, under the probabilistic ODE framework, is able to reproduce the day-to-day observations on those observed links, as shown in Figure 20. Again, the variance of the marginal distributions of most observed links is underestimated (or estimated as zeros) as expected, due to Lasso regularization. In addition, the estimation seems robust to a few outliers identified in Figure 18. Overall, the results of the probabilistic ODE are compelling and satisfactory. However, we speculate that the initial variance/covariance matrix can be critical to the final estimation results. Some prior knowledge about the variance/covariance of demand among O-D pairs can be obtained from traditional planning models, which may help improve the estimation results for real-world networks.

7 Conclusions

This paper develops a novel theoretical framework for estimating the mean and variance/covariance matrix of O-D demand considering the day-to-day variation induced by travelers’ independent route choices. The essential idea is to see the traffic data on each day as one data point and to use data points collected years along to estimate the probability distribution of O-D demand, as well as the probability distributions of links, paths and their generalized costs. As opposed to a real-valued estimation of flow and costs from traditional ODE, the probabilistic ODE estimates their probability distributions that are central to reliable network design, operation and planning. The probabilistic ODE framework is large-scale data friendly in the sense that it can make the best use of large-scale day-to-day traffic data to support complex decision making.

The framework estimates O-D mean vector and variance/covariance matrix iteratively, also known as iterative generalized least squares (IGLS) in statistics. IGLS holds great potential to converge faster than a formulation that estimates both simultaneously. It also decomposes a complex estimation problem into two sub-problems, which are relatively easier to solve. In the sub-problem of estimating O-D demand mean, we illustrate how to incorporate day-to-day traffic flow observations into the formulations and explain how the data size and historical O-D information effect the estimation results. In the sub-problem of estimating the O-D demand variance-covariance matrix, a convex optimization formulation is presented to approximate the solution. Lasso regularization is employed to obtain sparse covariance matrix for better interpretation and computational efficiency. We also discuss the observability of the probabilistic ODE problem. The non-uniqueness property of the probabilistic ODE under the IGLS framework is examined. Though probabilistic ODE works with a much larger solution space than the deterministic O-D ODE, we show that its estimator for O-D demand mean is no worse than the best possible estimator by an error that reduces with the increase in sample size.

The probabilistic ODE is examined on two small networks and two real-world large-scale networks. The solution converges quickly under the IGLS framework. In all those experiments, the results of the probabilistic ODE are compelling, satisfactory and computationally plausible. We also conduct the sensitivity analysis of estimation performance with respect to data sample size and historic O-D information. Increase the sample size and quality of historical O-D information can effectively approach the “true” probability distribution of O-D demand and path/link flow. Lasso regularization on the covariance matrix estimation leans to underestimate variance and covariance. A proper Lasso penalty ensures a good trade-off between bias and variance of the estimation. In practice, trial-and-error may be needed to identify a proper Lasso penalty.

In the near future, we plan to address a few computational issues before it can be widely deployed for practitioners. We will intensively test this probabilistic ODE method in other large-scale networks with a better data coverage than the D.C. network tested in this paper. Various modeling settings need to be tested, such as different route choice models, with and without Lasso regularization, with and without historical O-D, with and without traffic speed data. In addition, we speculate that the initial variance/covariance matrix can be critical to the estimator. Some prior knowledge about the variance/covariance of demand among O-D pairs can be obtained from traditional planning models. We plan to test how this prior knowledge can help improve the estimation results. In addition, this paper assumes the traffic observations on a set of days are i.i.d, but the set can be flexible. We can fit two probabilistic O-D distributions using solely workdays and weekends. We can even construct an unsupervised learning mechanism to cluster the traffic observations and then fit probabilistic O-D distribution for each cluster.

Furthermore, we plan to extend this research to estimate the probability distributions of time-varying O-D demand where mesoscopic traffic flow dynamics can be incorporated into the network modeling instead of naive BPR functions.

Acknowledgement

This research is funded in part by Traffic 21 Institute and Carnegie Mellon University’s Mobility21, a National University Transportation Center for Mobility sponsored by the US Department of Transportation. The contents of this report reflect the views of the authors, who are responsible for the facts and the accuracy of the information presented herein. The U.S. Government assumes no liability for the contents or use thereof. We would also like to thank anonymous reviewers for their valuable suggestions.

Appendix A Estimating the link flow on observed links

Identifying the estimator for the link flow on observed links can be cast into a quadratic optimization problem as shown in Equation 38.

[TABLE]

where

[TABLE]

First the above optimization problem is strongly convex, since $\Sigma_{x}=\Delta^{T}\Sigma_{f|q}\Delta+\Delta P\Sigma_{q}\Delta^{T}P^{T}+\Sigma_{e}$ and $\Sigma_{e}$, $\Sigma_{f|q}$ semi-positive definitive and $\Sigma_{e}$ is positive definitive. $\Sigma_{x}\succ 0$ always holds, and so does $\Sigma_{x}^{-1}\succ 0$

Proposition 3 implies that $\hat{x^{o}}$ can be derived in a closed form. This is proven by using the KKT condition of the minimization problem above.

Proof.

[TABLE]

Since ${\mathsf{x}}_{i}^{o}\succ 0$, $x^{o}=\frac{1}{n}x_{i}>0$. Then, $\lambda_{a}=0$, all four KKT conditions are satisfied. Since the minimization problem is strongly convex, therefore $\hat{x^{o}}=\frac{1}{n}{\mathsf{x}}_{i}^{o}$ is the optimal solution to the optimization problem. ∎

Appendix B Estimating the covariance matrix with Lasso regularization

We discuss the procedure for estimating the covariance matrix with Lasso regularization. We first use the proximal method Iterative Shrinkage-Thresholding Algorithm (ISTA) [56, 5]. ISTA can be modified to the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [57]. FISTA’s convergence rate is of $O(\frac{1}{n})$, compared to $O(\frac{1}{n^{2}})$ for ISTA.

Rewrite Formulation 29.

[TABLE]

where

[TABLE]

where $\mathcal{C}$ is a constant independent of $\Sigma_{q}$.

Using results presented in Petersen et al. [65], the derivative of $f(\Sigma_{q})$,

[TABLE]

Then define the soft-thresholding operator $S_{\lambda}(\beta)$:

[TABLE]

ISTA’s updating rule is:

[TABLE]

FISTA’s updating rule is:

[TABLE]

When applying the updating rules, proper step sizes need to be chosen to ensure $\Sigma_{q}\succeq 0$, backtracking methods can be used to select the step sizes [9]. Convergence analysis for both ISTA and FISTA can be found in Nesterov [56, 57].

Appendix C Proof of Proposition 9

Denote ${\cal X}$ as the sample space where we observe the data from, $x$ is one possible observation from the sample space. The ODE method $d(x)$ takes the observations $x$ as input and outputs the O-D estimators. The loss function $L(\cdot,\cdot)$ measures the discrepency between the estimated value and the true value. The risk is defined as the expectation of the loss between the true O-D demand mean $q$ and the estimated O-D demand mean by all possible observations. The intuition to define risk is that, we want to minimize the loss for all possible observations, rather than one specific observation (the empirical loss for specific numerical experiments). ODE methods that minimize the risk will have robust performance, even provided with noisy and limited data inputs. The risk function is defined by,

[TABLE]

The loss function can be any non-negative and strictly convex function. In this study, we use the quadratic loss since it is the most commonly used. However, the following proof would generally work for other norm operators.

[TABLE]

We first consider Formulation 18 without historical O-D information,

[TABLE]

where $\hat{\Sigma}_{x}^{-1}$ is the estimated inverse link variance/covariance matrix. $\hat{\Sigma}_{x}^{-1}$ can be estimated during the IGLS iteration, we rewrite the formulation in the form of $L^{2}$ norm,

[TABLE]

Formulation 61 is a standard non-negative least square problem. Many efficient algorithms [42] are practically ready to solve the formulation. Here we propose a fairly good estimation of O-D demand mean $q$ in Lemma 1 rather than directly solving for Formulation 61. We will later see that the proposed estimator $\hat{q}$ approximates the true O-D mean $q$ when the data size $n$ is sufficiently large.

Lemma 1**.**

When the data size $n$ is sufficiently large, the O-D demand mean can be estimated by,

[TABLE]

where $v^{o}=\hat{\Sigma}_{x}^{-\frac{1}{2}}x^{o}$, $\tilde{D}=\hat{\Sigma}_{x}^{-\frac{1}{2}}\Delta^{o}p$ and $\tilde{D}^{+}$ is the Moore–Penrose pseudoinverse of matrix $\tilde{D}$. Detailed proof can be found in Nie et al. [60].

If $|K|\leq|A^{o}|$, then $\tilde{D}^{+}=\tilde{D}^{-1}$. Otherwise, $\tilde{D}^{+}=(\tilde{D}^{T}\tilde{D})^{-1}\tilde{D}^{T}$. Based on Lemma 1, we prove Proposition 9.

Proof.

Our first target is to bound the risk when the data sample size increases. The risk for a given estimation method $d(x)$ presented in 62 is,

[TABLE]

In Equation 67, $\tilde{p}^{T}{\Delta^{o}}^{T}\hat{\Sigma}_{x}^{-1}\Delta^{o}\tilde{p}$ is invertible when $\Delta^{o}$ is fully ranked. Since $\left\lVert p^{T}{\Delta^{o}}^{T}\right\rVert_{2}^{2}$ is independent of the data size $n$, we can see it as a constant. As for $\left\lVert\left(\tilde{p}^{T}{\Delta^{o}}^{T}\hat{\Sigma}_{x}^{-1}\Delta^{o}\tilde{p}\right)^{-1}\right\rVert_{2}^{2}$,when the sample size increases, $\hat{\Sigma}_{x}$ approximates $\Sigma_{x}$. As long as the observed data $x$ is bounded, $\hat{\Sigma}_{x}$ can be bounded, independent of the sample size $n$. For $v^{o}$, in the sub-problem of estimating the O-D mean vector, we have,

[TABLE]

Again $\hat{\Sigma}_{x}^{-\frac{1}{2}}\Sigma_{x}\hat{\Sigma}_{x}^{-\frac{1}{2}}$ can be bounded. When $n\to\infty$, by Law of large number (LLN), we have,

[TABLE]

Also after assuming $\left\lVert\left(\tilde{p}^{T}{\Delta^{o}}^{T}\hat{\Sigma}_{x}^{-1}\Delta^{o}\tilde{p}\right)^{-1}\right\rVert_{2}^{2}\left\lVert\tilde{p}^{T}{\Delta^{o}}^{T}\right\rVert_{2}^{2}\left\lVert\hat{\Sigma}_{x}^{-\frac{1}{2}}\Sigma_{x}\hat{\Sigma}_{x}^{-\frac{1}{2}}\right\rVert_{2}^{2}\leq M$, we have

[TABLE]

This implies that as long as the estimation of $\Sigma_{x}$ is bounded, any estimator $d(x)$ can achieve the same level of accuracy provided with a sufficiently large $n$.

For Formulation 22, suppose we use heuristic methods to solve the bi-level formulation as in Yang [92]. Each iteration in solving the upper level problem is equivalent to solving Formulation 18 with certain route choice probability $p$. Note the bound applies for all route choice probability $p$. Therefore the statistical risk of the estimated O-D mean is still of ${\cal O}\left(\frac{1}{n}\right)$. ∎