Statistical inference of travelers' route choice preferences
with system-level data

add1]Pablo Guarda

add1,add2]Sean Qian

[add1]Department of Civil and Environmental Engineering, Carnegie Mellon University
[add2]Heinz College, Carnegie Mellon University

Traditional network models encapsulate travel behavior among all origin-destination pairs based on a simplified and generic travelers' utility function. Typically, the utility function consists of travel time solely and its coefficients are equated to estimates obtained from discrete choice models and stated preference data. While this modeling strategy is reasonable, the inherent sampling bias in individual-level experimental data may be further amplified over network flow aggregation, leading to inaccurate flow estimates. In addition, individual-level data must be collected from surveys or travel diaries, which may be labor intensive, costly and limited to a small time period. To address these limitations, this study extends classical bi-level formulations to estimate travelers’ utility functions with multiple attributes using system-level data. This data tends to be less subject to sampling bias than individual-level data. It is cheaper to collect and it has became increasingly diverse and available. To leverage system-level data, we formulate a methodology grounded on non-linear least squares to statistically infer travelers' utility function in the network context using traffic counts, traffic speeds, the number of traffic incidents and sociodemographic information obtained from the US Census, among other attributes. The analysis of the mathematical properties of the optimization problem and of its pseudo-convexity motivate the use of normalized gradient descent, an algorithm developed in the machine learning community that is suitable for pseudo-convex programs. More importantly, we develop a hypothesis test framework to examine statistical properties of coefficients attached to utility terms and to perform attributes selection. Experiments on synthetic data show that the coefficients of the travelers’ utility function can be consistently recovered and that hypothesis tests are a reliable statistic to identify which attributes are determinants of travelers’ route choices. Besides, a series of Monte-Carlo experiments showed that statistical inference is robust to noise in the Origin-Destination matrix and in the traffic count measurements, and to various levels of sensor coverage. The methodology is also deployed at a large scale using real-world multi-source data in Fresno, CA collected before and during the COVID-19 outbreak.

network models, stochastic user equilibrium, route choices, travel behavior, utility function, multinomial logit model, hypothesis testing, pseudo-convexity, normalized gradient descent, traffic count data

§ INTRODUCTION

System-level data, as opposed to individual-level data, measures characteristics of aggregated flow in transportation networks, and it is a valuable source of information for studying travel demand. An example is the use of traffic count data for estimating origin-destination (O-D) matrices <cit.>. Interestingly, few attention has been given to use of system-level data to estimate travelers' utility functions in route choice models. These models are key to depict travel demand in transport planning applications and their parameters are usually estimated with data collected from individual-level surveys and experiments. Modeling route choices at the system level using  individual-level data, however, may be inaccurate and expensive due to the inherent sampling bias and high collection cost of experimental data. The increase in availability and diversity of system-level data offers opportunities to overcome those limitations and to understand the impact of a broader set of factors that may influence travelers' route choice decisions. Some of the new sources of system-level data include crashes, weather and pavement conditions, land use characteristics, socio-demographics information from Census, travel time reliability and other relevant attributes. Clearly those factors may determine how a traveler makes his/her route choice, but surveying those for each individual traveler would be infeasible. Thus, this paper proposes a statistical method to better understand why and how a traveler make route choices using a diverse set of system-level data only.

Researchers have already leveraged system level data within traditional network models to estimate travelers' route choice preferences  <cit.>. However, their primary focus has been on the estimation of a utility function dependent on travel time only. Furthermore, most existing algorithms has been deployed to network of relatively small size, which raises questions on their potential to scale up to large transportation networks and to contribute to real world applications. A main difficulty to solve this estimation problem comes from the fact that network models are not data-driven. Discrete choice models, the gold standard to analyze individual-level data in travel behavior studies, may be considered a promising alternative to solve this problem given their data-driven nature. However, these models are not able to capture the endogeneity of travel times arising from the interdependence between travelers' decisions and from the effect of traffic congestion in transportation networks.

With the goal of estimating the coefficients of travelers' utility functions with multiple attributes and using system-level data, this paper enhances traditional network models with optimization algorithms developed in the machine learning literature. Mathematical properties of the optimization problem and its pseudo-convexity motivates the use of normalized gradient descent, a first order method that is suitable for pseudo-convex optimization. Our solution algorithm is tailored to perform a robust estimation of the coefficients of the travelers' utility function and to also reproduce the traffic conditions observed at the system level. To ease the identification of the determinants of travelers' route choices with system level data, we formulate hypothesis tests on the utility function coefficients.

The paper is structured as follows. We first survey the existing literature and identify the main contributions of our work. Then, we present an example on a toy network to illustrate the optimization problem arising from the estimation of the travelers' utility function coefficients using system level data. The following three sections describe the formulation of our methodology on a general transportation network, analyze the mathematical properties of the optimization problem and propose a solution algorithm to solve the problem. Subsequently, we describe the framework to perform statistical inference on the utility function coefficients and we present the results of numerical experiments conducted in networks of small and medium size. Next, we show the estimation results obtained in a large scale network and using real world system level data. Finally, we describe our main conclusions and suggests avenues for further research. The mathematical notation used for the remainder of the paper is included in <ref>.

§ LITERATURE REVIEW AND RESEARCH GAPS

The computation of network equilibrium requires to specify a travelers' route choice model. The choice of route choice model defines the class of network equilibria <cit.>. The simplest and most traditional class of equilibria is known as deterministic user equilibrium (). At no traveler has incentive to change to alternative routes and travelers' are assumed to pick, in a deterministic manner, the route that maximize their utility. The multinomial model () remains as a gold standard in travel behavior research to depict travelers' decision making <cit.>. A well-known application of the model in network modeling is for the computation of Stochastic User Equilibrium with assignment (-). In contrast to other probabilistic models of travelers' choices, the model exhibits a better compromise between behavioral realism and mathematical tractability. Similar to the model, the -enjoys of this mathematical convenience in the network modeling context. For instance, under mild conditions, -has solution uniqueness in link and path flow space whereas only has solution uniqueness in link flow space. Besides, the existence of a closed form expression for the choice probabilities in the multinomial model can be leveraged for the computation of -, e.g the Dial algorithm <cit.> or to solve the problem <cit.>.

A stream of research that brought our attention concerns the estimation of the travelers' utility function from traffic count data and which induces a traffic assignment consistent with -<cit.>. The solution of this problem and of network equilibrium in general requires to know a priori the coefficients of the travelers' utility function of the route choice model. Hence, a standard practice is to set the values of these coefficients equal to estimates obtained in previous travel behavior studies. Nevertheless, there are multiple advantages to estimate these coefficients. First, it avoids to search for external estimates which may be cumbersome and it may provide inaccurate coefficients. Second, it can improve the generalization performance on the estimated O-D matrix when, for instance, the coefficients fitted from existing count data become close to the population coefficients. The following sections describe the relevant literature that has studied the aforementioned problem.

§.§ Overview of the Logit Utility Estimation (LUE) problem

Seminal work in the literature focused on the problem of estimating the coefficient $\hat{\theta} \in \sR$ of a utility function dependent solely on travel time and where both the O-D matrix and travel costs among links/paths are assumed exogenous. To our knowledge, the problem was first studied by <cit.>, who estimated the coefficient $\hat{\theta}$ using link and path flows obtained from a traffic assignment consistent with the <cit.> method. The estimator $\hat{\theta}$ was assumed to be Gaussian distributed and it was obtained via maximum likelihood estimation (). A main limitation of this work was the requirement of knowing the traffic counts at every link/path in the network. The problem of estimating $\hat{\theta}$ from a subset of traffic count in the network was then addressed by <cit.>. Another major limitation of this work was assuming full knowledge of the average cost of traversing paths between O-D pairs, which is plausible in settings where costs are equated to travel times but arguably unrealistic in real world scenarios where travelers make route choices based on multiple attributes such as monetary cost and travel time reliability.

<cit.> focused on a more general setting where link flows observations were consistent with but not necessarily with a assignment. In contrast to <cit.>, the distribution of the estimator $\hat{\theta}$ was derived by imposing distributional properties on the link flows and by relying on the statistical properties of the estimator. This work is the first that conducts hypothesis testing on the coefficients of the travelers' utility function. However, similar to previous research, it assumes a utility function dependent on travel time only and exogenous link costs, which limits the application to uncongested networks or to settings where link costs at equilibrium are known. Under the assumption that link flow data is consistent with -, <cit.> extended prior work by allowing for endogenous travel costs. The authors studied the impact of estimating $\hat{\theta}$ without accounting for the endogeneity of the travel costs and observed that estimates of $\hat{\theta}$ become biased and that statistical inference is less consistent. Besides, the estimation of $\hat{\theta}$ was shown to perform better with least squares than maximum likelihood.

§.§ Overview of the Origin-Destination and Logit Utility Estimation (ODLUE) problem

In the following decades, the research interest shifted to addressing the joint estimation of the O-D matrix and $\hat{\theta}$ from traffic count data following -, which we coined as O-D and Logit Utility estimation () problem[We denominated this problem as O-D Logit Utility estimation () because the O-D matrix is estimated in addition top the travelers' utility function coefficients as it is done in ]. Interestingly, the extension of the problem to settings where the utility function was dependent on multiple attributes besides travel time did not gain the same attention. By this time, the literature had made significant progress to estimate O-D matrices with traffic count data consistent with -and these advances could be directly leveraged to solve the problem. A rich set of solution methods had been developed in the literature, which include sensitivity analysis <cit.> and the alternating optimization of the upper and lower levels of the classic bilevel formulation <cit.>.

To our best knowledge, <cit.> conducted the first study that studied the problem and that presented an application using system level data collected from a real transportation network. Traffic count and travel time data was collected from a small network at the Purdue University campus and estimates of $\hat{\theta}$ were obtained for different time periods of the day. The authors proposed a two stage calibration method that leveraged the closed form of the route choice probabilities to solve for $\hat{\theta}$ using the secant method. The O-D matrix was estimated via a least square minimization. Key limitations of this work were to not account for the congestion effects in the transportation network and to assume that traffic counts and travel time data were available for the full set of links. This assumption could be mild in small networks where traffic counts can be collected at every link but implausible in larger networks where this data is typically available for a small proportion of the links.

<cit.> overcame some of the aforementioned limitations by solving a bilevel optimization problem that resembled the traditional mathematical program with equilibrium constraints () that had been used for solving the problem. A main contribution of this work was to estimate both $\hat{\theta}$ and the O-D matrix while also accounting for the impact of traffic congestion. A Sequential Quadratic Programming () algorithm was used to find the parameters of interest and to minimize the gap between observed and predicted traffic counts. The optimization problem also included constraints to ensure that the -equilibrium conditions were satisfied over parameters' updates. The authors also derived for the first time the analytical expression of the gradients of the outer level objective respect to the parameters associated to the O-D matrix and to the travelers' utility function. <cit.> proposed a similar framework but where a instead of a non linear least square () problem was solved at the upper level of the bilevel formulation. The authors performed on both synthetic and real world data collected from the Tuen Mun Corridor network in Hong Kong. Later, <cit.> implemented a similar approach using real world data gathered from a small network in Seattle, WA. They also performed experiments on synthetic data to study the robustness of the parameter estimates to noise in the ground truth O-D matrix and the parameter $\hat{\theta}$ used to generate synthetic traffic count data.

§.§ Extensions of the O-D Logit Utility Estimation (ODLUE) problem

A number of studies have estimated additional parameters on top of the O-D matrix and of the travelers' utility function coefficients.  <cit.> and <cit.> estimated the shape parameters of the link performance function using traffic counts and travel time measurements consistent with -.  To leverage the travel time data, the outer level objective of their problem incorporated the squared difference between predicted and observed travel times. Results on synthetic data suggested that all parameters could be consistently recovered. However, few attention was paid to the identifiability of the parameters and on overfitting due to the additional degrees of freedom introduced in the model. In fact, is known to be an underdetermined problem which is subject to identifiability issues even when a full set of traffic counts is available in a network <cit.>. Naturally, this issue is worsen if additional parameters are introduced in the model. Besides the problem of identifiability (also known as observability), a high model complexity may cause overfitting and thus, a poor generalization of the model on unseen data even if a good fit is observed on training data. <cit.> also argued that in an uncongested network, the optimization problem is convex to claim solution uniqueness. However, no attention was paid to the non-convexity of the feasible set in the optimization problem, i.e. which makes the problem non-convex, and also to the fact that convexity is not a sufficient condition for solution uniqueness.

Other line of research extended the application of the beyond the traffic assignment stage of the the classic four step model <cit.>. Under the assumption of an exogenous demand matrix, <cit.> simultaneously estimated the parameters of a category index trips generation model, the parameters associated to the distribution step and the mode specific constants, and the travel time and cost coefficients of the utility function of a multinomial logit model of mode choice. The authors tested their method with real data collected from two Italian cities and found that the magnitude and sign of the estimated parameters were consistent with their expectations. However, the treatment of congestion in the road network in their work to capture the endogeneity of travel times is missing. In similar research, <cit.> estimated the utility function parameters of a route choice model and the mode choice parameters using a nested logit model via Nonlinear Generalized Least Squares (). A bilevel optimization problem was formulated to account for traffic congestion and it was solved via an alternating optimization of the upper and inner level problems. Similar to most prior work, no statistical tests were conducted to check if those parameters were statistically significant or to assess if the increase in goodness of fit compensated the additional degrees of freedom introduced in the full model.

More recently, <cit.> used a computational graph to estimate the parameters of a model that incorporated the steps of trip generation, spatial distribution and path flow-based traffic assignment of the 4-step model. Data from household travel surveys, mobile phones and on traffic counts were leveraged to estimate these parameters and the model was deployed on a network with more than 2500 nodes and 5000 links. This research illustrated the potential of computational graphs to increase the scalability of models that have been developed in the literature and to also integrate heterogeneous data sources that inform about transport decisions. Nevertheless, a main drawback of this work is to not account for traffic congestion and to focus on a travelers' utility function dependent on travel time only.

§.§ Contributions of this research

This paper enhances existing formulations of the problem with the goal to statistically infer the coefficients of travelers' utility functions with multiple attributes using system-level data. A bilevel optimization program is formulated to both minimize the gap between estimated and observed traffic counts, and to account for the endogenous effect of traffic congestion on travelers' route choices. More importantly, a framework of hypothesis tests is proposed to examine statistical properties of the coefficient estimates. Under the pseudo-convexity of the optimization problem, we implement a variant of gradient descent suitable for pseudo-convex optimization. First order optimization methods seem promising to speed up computation and thus, to make model training scalable to large transportation networks.

Below are the four main contributions of this paper to existing literature:

* It presents a methodology to statistically infer the coefficients of travelers' utility functions consisting of multiple attributes and using system-level data. It does not only extend the problem, but it also enables learning realistic travel behavior from diverse datasets at the system level.

* It conducts for the first time a mathematical analysis of the non-convexity of the problem respect to the utility function coefficients estimated from traffic count data.

* For the solution of the problem, it shows that the integration of first order and second order optimization methods outperform existing approaches in the literature.

* It presents a statistical framework to perform hypothesis testing and attributes selection on the coefficients of multi-attribute utility functions

* It implements the methodology in a large-scale transportation network and with real system-level data.

All our analyses are replicable, open-sourced and ready to be used for the transportation community (Section <ref>).

§.§ Illustrative example

Consider a network with two parallel links and thus, with two alternative paths only (Figure <ref>). The costs functions associated to each link/path are $t_1(x_1) = t^{0}_{1}(1+\alpha(\frac{x_1}{\gamma_1})^\beta)$ and $t_2(x_2) = t^{0}_{2}(1+\alpha(\frac{x_2}{\gamma_2})^\beta)$, where $\gamma_1,\gamma_2 \in \sR$ are the link capacities and $t^0_1=t^0_2$ are the links' free flow travel times. The parameters of the link performance functions are $\alpha = 0.15$ and $\beta = 4$. Assume that the travelers' utility function is dependent of the travel time $t$ and the monetary cost $c$ of traversing each link of a path only. In addition, suppose that the coefficients $\theta_t, \theta_c \in \sR$ are linearly weighting each attribute of the utility function and that they are common among travelers. Thus, the deterministic component of the utility attained to each link/path can be expressed as $v_1 (x_1) = \theta_t t_1(x_1) + \theta_c c_1$ and $v_2 (x_2) = \theta_t t_2(x_2) + \theta_c c_2$.

at (0,1) [circle, draw] (N1) $1$ ;
at (2,1) [circle, draw] (N2) $2$ ;

[draw, ->] (N1) edge[bend left = 60] node[below left = 0.4cm and 0.9cm]$a_1$  (N2) ;
[draw, ->] (N1) edge[bend right= 60] node[below]$a_2$ (N2);

at (1,1.6) [rectangle, draw, fill = white] (Toll) $\phantom{aaaa}$;

Illustrative example

Consider $\bar{q} \in \sR_+$ individuals traveling between origin destination 1-2 and making route choices consistent with a model. $\bar{q}$ is small enough such that links' travel times are approximately equal to their free flow travel times, i.e. $t_1(x_1) \approx t_1^0, \ t_2(x_2) \approx t_2^0$. Suppose link $a_1$ has a toll fee equal to 1 USD, i.e $c_1 = 1, c_2 = 0$ and a sensor recording traffic counts. Then, path choice probabilities can be modeled with a sigmoid function and link flows can be obtained as follows:

\begin{equation}
\label{eq:link-flow-solution-illustrative-example}
x_1 = \frac{\exp(\theta_t t_1+\theta_c c_1)}{\exp{(\theta_t t_1+\theta_c c_1)}+\exp{(\theta_t t_2+\theta_c c_2)}} \ \bar{q} =  \sigma(\theta_c) \bar{q},\quad x_2 = \bar{q} -x_1 = (1-\sigma(\theta_c))\bar{q} = \sigma(-\theta_c)\bar{q}
\end{equation}

Suppose the goal is to estimate the coefficients $\theta_t$ and $\theta_c$ of the travelers' utility function with the traffic count measurement $\bar{x}_1$. From the link flow solutions (Eq. <ref>), we note first that $\theta_t$ is not identifiable because travel times are the same in the two alternative paths. In contrast, because the monetary costs are different between paths, $\theta_c$ is identifiable and equal to the solution $\theta^{\star}_c$ of the following nonlinear least square minimization problem:

\begin{align}
\label{eq:objective-function-illustrative-example}
\theta^{\star}_c
= \arg \min_{\theta_c } \ (x_1(\theta_c) - \bar{x}_1)^2
%= \arg \min_{\theta_t } \left[ \left(\frac{\bar{q}}{1+\exp{(-\theta_t)}} - \bar{x}_1\right)^2 + \left(\frac{\bar{q}}{1+\exp{(\theta_t)}} - \bar{x}_2\right)^2\right]
%&= \arg \min_{\theta_t } \left[ \left(\sigma(\theta_t) - \frac{\bar{x}_1}{\bar{q}}\right)^2 + \left(\sigma(\minus\theta_t) -_{} \frac{\bar{x}_2}{\bar{q}}\right)^2\right]
= \arg \min_{\theta_t } \left(\sigma(\theta_c) - \frac{\bar{x}_1}{\bar{q}}\right)^2
\rightarrow
\theta^{\star}_c = \ln(\bar{q}/\bar{x}_1-1)
\end{align}

Consider two scenarios for the solution of this problem (Eq. <ref>). In the first scenario $\bar{x}_1 = 0.5 \bar{q}$ meaning that $\bar{x}_1$ is deterministic and that $\theta_c^\star = 0$. In the second scenario, $\bar{x}_1 = x_1 + \epsilon$ , $\epsilon \sim \mathcal{N}(0,1)$, i.e. $\bar{x}_1$ is not deterministic and $\bar{x}_1 \sim \mathcal{N}(0.5\bar{q},1)$. Note that a single realization of $\bar{x}_1$ in the second scenario could lead us to conclude with the same probability that the traveler's utility and the monetary cost are positively or negatively associated. However, in expectation and similar to the first scenario, this association does not exist.

This example illustrates important features of the problem studied in this paper. First, it shows that the non-linear least squares formulation of the problem is suitable to estimate the coefficients of the travelers' utility function with traffic count data. Second, the identification issue associated to $\theta_t$ illustrates the special considerations that must be made to specify the travelers' utility function and to prevent a naive estimation of coefficients. Third, the estimation of $\theta_c$ in the non-deterministic scenario illustrates the relevance of statistical inference tools to assess if a feature is a determinant of the travelers' route choices in the presence of randomness of the data generating process.

For larger networks, the analytical derivation of the optimal solution for travelers' utility function parameters is intractable and thus methods that can iteratively approach to the optimal solution are needed. Besides, due to the class of non-linearity of the link flows respect to the utility function coefficients, the objective function cannot be guaranteed to be convex, which imposes additional difficulties to search for an optimal solution. Another challenge is related to the impact of traffic congestion on travelers' choices and which is not addressed in the illustrative example.

§ ESTIMATION OF TRAVELERS' UTILITY FUNCTION COEFFICIENTS WITH SYSTEM LEVEL DATA

This section discusses in detail our methodology to estimate the travelers' utility function coefficients using system level data. We first introduce our main assumptions and then we describe the components of the bilevel optimization program formulated to estimate the utility function coefficients.

§.§ Assumptions

Network traffic flow follows stochastic user equilibrium with logit assignment

With the motivation of bridging transportation network analysis and the study of travel behavior with discrete choice models, we focused on the problem of estimating the travelers' utility function from link flow measurements consistent with stochastic user equilibrium under assignment (-). On one hand, -is one of the many alternative representations of user equilibrium used in the networking modeling community, hence it may be considered a strong assumption. On the other hand, -is in line with the state of the practice in travel behavior research where the model remains as one of the gold standards to depict individual route choices. The good compromise between behavioral realism and the mathematical tractability of is a key property for the estimation in discrete choice models and it will be also leveraged in the formulation and solution of our problem.

[Exogenous and deterministic O-D]

The origin-destination (O-D) demand matrix is deterministic and exogenous

O-D demand estimation () is a complex problem that has been largely studied in the network modeling literature. As an underdetermined problem, there is no guarantee on solution uniqueness for the estimated O-D matrix and this makes difficult to verify the consistency of the solution with a ground truth O-D matrix. The literature has extended to estimate the travelers' utility function coefficients on top of the O-D matrix. The problem is harder to solve due to the additional degrees of freedom introduced by an endogenous O-D matrix and therefore, as , it also suffers of solution non-uniqueness.

The assumption of an exogenous and deterministic O-D matrix allows us to focus on the estimation of the utility function coefficients, a first step to understand the mathematical properties of our problem. It can always be easily extended with any method to solve for . In the context of the existing literature, our methodology can be seen as a case where a reasonably accurate reference O-D matrix is available and assumed as the ground truth O-D matrix. As we will show later, we can still obtain satisfactory results if the exogenous O-D demand is noisy but not too far from the ground truth O-D.  Note also that our methodology can be extended to an iterative process where and alternate if the goal is to simultaneously estimate both the travelers' utility function coefficients and the O-D matrix.

[Linearity and homogeneity of utility function] The utility function is a linear weight of attributes and coefficients and the coefficients are common/homogeneous for all individuals traveling in the transportation network.

The choice of a linear-in-parameters utility function and with homogeneous coefficients is a standard assumption for the estimation of discrete choice models in travel behavior studies. The preference of linear over non-linear specifications of the utility function is motivated by parsimony arguments and on the convenience of dealing with a concave likelihood function in the estimation of multinomial logit models () models. In our case, this assumption significantly facilitates the proofs of theoretical guarantees in our optimization problem. Note that this assumption could be relaxed by allowing a variation of the utility function coefficients at the O-D pair level or by letting the coefficients to be drawn from some arbitrary probability density function. Both modeling strategies can be seen as analogies to the systematic taste variation used in models and to the mixing distributions used in Mixed Logit models <cit.>.

[Monotonic increase of and exogeneity of parameters of link performance functions] The link performance functions are monotonically increasing respect to their link flows and their parameters are exogenous

The exogeneity of the parameters of the link performance is a standard assumption in transportation network analysis. Following the state of the practice in prior literature, we adopt the Bureau of Public Roads (BPR) link travel time function. For each link, the travel time monotonically increases respect to the link flow. For the sake of simplicity, we also assume that there is no interaction between links and hence, the travel time of a link is dependent on the traffic flow in that link only. A direct consequence of this assumption is the uniqueness of the link flow solution of -. Some work in the literature has proposed methods to estimate the parameters of the link performance function, namely, relaxing the exogeneity assumption. Although it is feasible to estimate these parameters on top of the travelers' utility function coefficients within our modeling framework, the introduction of these additional degrees of freedom may cause identifiability issues or it could make the statistical inference of the utility function coefficients less reliable.

[Constrained path set]

The set of paths between each O-D pair has a fixed size and it is dynamically updated over iterations with column generation methods

Most prior work on the problem assumes that the consideration set in an O-D pair is equal to the set all reasonable paths in that pair <cit.>. Under this assumption, each path choice probability can be decomposed into a weight of exponentials of their link utilities. On one hand, this property avoids the enumeration of all paths in the network and it can significantly speed up the computation of -that it is required to estimate the travelers' utility function coefficients. On the other hand, it implicitly assumes that travelers' consider all paths to travel between an O-D pair and hence, it may induce overlapping among paths. Path overlapping is known to affect the computation of path choice probabilities in multinomial logit models () and it has been identified as a drawback of algorithms that rely on the set of all reasonable paths <cit.>.

Navigation apps are widely used nowadays and they recommend a constrained set of options to travel between an O-D pair. This makes us believe that in modern real world applications the assumption of a constrained path set between O-D is more behaviorally plausible than considering all reasonable paths. Besides, a constrained path set is also expected to induce less path overlapping, which is detrimental route choice models. A challenge of working with a constrained path set is the impossibility of performing a link level decomposition of the path choice probabilities and hence, the potential increase of computational cost in the network loading stage of -. In addition, while path sets can be reasonable approximated with existing information, e.g. Google Maps recommendation, they may be inaccurate or not available for all O-D pairs. To improve our prior about the true composition of the path set, the path sets are dynamically updated via column generation methods <cit.>. The size of the consideration sets is treated as a hyperparameter and it is constrained to small values to reduce computational burden.

To capture the correlation between paths in the consideration set due to shared link segments, path utilities are corrected with the same principle used in the path size logit (PSL) model <cit.>. Path utilities incorporate the logarithm of a factor that increase with the amount of overlapping among paths within the same consideration set. A path with no overlapping links needs no utility adjustment since the PSL factor has a size of one. Thus, the PSL correction reduces the chances of a violation of the Independence of Irrelevant Alternatives (IIA) assumption of the model. Previous literature suggests to generate paths that do not overlap with existing paths during the column generation phase <cit.>. However, this approach requires to solve a combinatorial problem on all possible subset of paths and it implicitly assumes that travelers' consideration set have paths with few or no overlap.

§.§ Stochastic user equilibrium with logit assignment (-)

In a congested network, travel times are endogenous and dependent of the traffic flows. At equilibria, both traffic flows and travel times are expected to reach an stationary point where travelers have no incentives to switch to alternative paths. Depending on the underlying behavioral representation used to model travelers' route choices, different types of network equilibrium are induced. In Deterministic User Equilibrium (), travelers' choices are deterministic, meaning that the modeler fully knows the specification of the traveler' utility function. In contrast, in stochastic user equilibrium (), travelers' choices are assumed to be probabilistic and modelers' may be ignorant of a set of unobservable components of the traveler' utility function.  In the network context, at , no user believes he can improve his travel cost by unilaterally changing routes <cit.>. Under with assignment (-), travelers' are assumed to make choices consistent with a multinomial model ().

§.§.§ Multinomial Logit () route choice model

The model remains as the gold standard to model route choices due to its good compromise between behavioral realism and mathematical tractability (McFadden, 1973). Consider a traveler $l \in L$ choosing a path $i$ within her consideration set $J_{l}$ and based on the latent utility $U_{jl}$ of each alternative path $j \in J_{l}$. Suppose that the modeler knows the observable component $V_{jl}$ of the travelers' latent utility but ignores a component $e_{jl}$ which can be safely assumed to be stochastic. Therefore, if any attribute relevant to route choice decision is unobserved, travelers' choices could look probabilistic but not deterministic from the modelers' perspective. In particular, if $e_{jl} \overset{\text{i.i.d}}{\sim} \ EV(0,\mu)$ the probability $p_i$ that a traveler $l$ chooses a path $i \in J_{l}$ will have the following closed form:

\begin{align}
\label{eq:mnl-choice-probabilities}
\displaystyle p_i = \dfrac{\displaystyle \exp(\mu V_i)}{\displaystyle \sum_{j \in A_{l}} \exp(\mu V_j) }
\end{align}

where $\mu > 0$ is a scale parameter of the Extreme Value (EV) Type 1 distribution and which is set to 1 for convenience and identification purposes. $\mu$ is inversely proportional to the variance of the random component, and thus, it is expected to be lower as smaller is the unobserved component of the latent utility. At a extreme case where the latent utility function is fully observed, $\mu \to \infty$ and the choice probabilities reduce to an indicator function taking the value 1 if the utility of an alternative is the highest within a choice set and 0 otherwise, i.e. the case.

§.§.§ SUE-logit with multi-attribute utility function

The original formulation of the -problem entails that route choices are made based on travel cost/time only <cit.>. Utility, however, may depend on additional attributes such as travel time reliability, waiting time or monetary costs (e.g. toll fees). To make a more explicit bridge between the travelers' utility function specified in route choice models and to extend the analysis for a multi-attribute case, we reformulated the -problem as follows:

,   ∑_a ∈A∫_0^x_a v_a(u,) du - ⟨, ln⟩

where $\vx \in \sR_{\geq 0}^{|A|}, \vf \in \sR_{\geq 0}^{|H|}, \vq \in \sR_{+}^{|W|}, \mIq \in \sR_{\{0,1\}}^{|W| \times |H|}, \ \mIx \in \sR_{\{0,1\}}^{|A| \times |H|}$. Besides, $v_a (u, \vtheta) = \theta_t t_a(u) + \sum_{k \in K_{\mZ}} \theta_k \cdot Z_{ak}$ is the utility associated to link $a$ and at a traffic flow level $u$, $Z_{ak}$ is the value of the exogenous attribute $k \in K_{\mZ}$ at link $a$, $K_{\mZ}$ is the set of exogenous attributes at each link, $\vtheta_Z \in \sR^{|K_Z|}$ and $\vtheta = \begin{bmatrix} \theta_t & \vtheta_Z\end{bmatrix} \in \sR^{|K|}$ are the vector of coefficients associated to the exogenous attributes and to all attributes, respectively. Exogenous attributes may include the number of traffic lights, streets intersections or the level of income of the area where a link is located. From the specification of $v_a(\cdot)$, it follows that travel time ($t$) is the only endogenous attribute in the utility function and it is linearly weighting the preference coefficient $\theta_t$. The first order optimality condition of Problem <ref> gives the following path flow solution:

\begin{equation}
\label{eq:sue-logit-path-flows-solution}
%	https://tex.stackexchange.com/questions/38352/resize-scale-equation-in-beamer
f^\star_{h} = q_{w} \frac{\displaystyle \exp\left( \sum_{a \in A} {v_a^{\star}} \delta_{ah}\right)}{\displaystyle\sum_{j \in H_{w}}^{} \exp\left(\sum_{a \in A} {v_a^{\star}} \delta_{aj}\right)}
%		\to
%\mathbf{f} = (\mIqT \vq) \circ \left(\exp\left(\displaystyle \mIxT \mathbf{v} \right) \oslash \left(\displaystyle   \mIqT \mIq \exp(\mIxT \mathbf{v})\right)\right)$
%		\vf = (\mIqT \vq) \circ \vpf(\vvx)
\end{equation}

where $\delta_{ah} = \mathbb{I}(\textmd{link } a \in \textmd{ path } h)$ and $v^\star_{a}$ is the link utility of link $a \in A$ at -. As expected, the path flow vector $\vf^\star_{h}$ at -follows a distribution. From a micro-level travel behavior standpoint, $q_{w}$ individuals traveling in the O-D pair $w \in W$ are making route choices according to a model (Eq. <ref>), and hence the path flow $\vf^\star_{h}$ is the aggregate of these individual decisions (Eq. <ref>). A detailed derivation of the extension of the <cit.> formulation from a single to a multi-attribute case is included in <ref>.

In line with the rules of parameterization used to derive the ,  the specification of the utility function defined to compute -should account for the fact that the utility function coefficients are scaled by a factor $\mu \in \sR_{+}$ proportional to the variance of the unobservable component of the utility function i.e. $\vtheta = \mu \tilde{\vtheta}$ where $\tilde{\vtheta}$ is the unscaled vector of coefficients and which is not identifiable. A direct consequence in the single attribute case is that only the sign but not the magnitude of the coefficients should assumed to be known by the modeler. This is particularly relevant in cases where estimates from external travel behavior studies are used to set the utility function coefficients of the route choice model defined for the computation of -. Therefore, accounting for the scale factor of the parameters here is critical for a correct economic and behavioral interpretation of the coefficients of multi-attribute utility functions but that is oftentimes overlooked in the and literature.

§.§.§ Stochastic network loading

Suppose the vector of the travelers' preferences $\vtheta \in \sR^{|D|}$, the values of the matrix of exogenous attributes $\mZ$ and the travel times at -are known and that the goal is to find the resulting path/link flows in the transportation network. This is precisely the problem that stochastic network loading () aims to solve. Assumptions about the composition of the consideration set can significantly speed up the computation of . A well-known example is the Dial algorithm which, under the assumption of path sets containing all reasonable paths, allows to decompose path probabilities into link level weights that can be then used to recursively obtain the resulting links flows in the transportation network. Alternatively, and in line with Assumption <ref>, we assume path sets with a constrained size and that are dynamically updated via a column generation method. The pseudo-code of our method is shown in Algorithm <ref>, <ref>.

§.§.§ Solution methods for -

If travel times were insensitive to variation of link flows, namely, exogenous, a single computation of would suffice to obtain the link and path flows at -. In practice, travel times are endogenous variables and thus, the computation of -requires to perform multiple times.

The method of successive average () is one of the prominent algorithms used to compute -in the literature. At the initial iteration, is computed to generate a feasible link flow solution. For each of the following iterations, the convex combination of the link flow solutions at the current and the previous iteration is computed to obtain a new feasible solution. defines a step size $\lambda_i = 1/(1+i)$ to set the weights $\lambda_i$ and $1-\lambda_i$ that are used to compute the convex combination of solutions at iteration $i$. Because the feasible set of the -is convex, the convex combination necessarily gives a feasible point. The process is repeated until some convergence criterion has been achieved, such as the difference between the current and previous feasible link flow solution.

Frank-Wolfe (F-W) is a general purpose method for constrained convex optimization and that can outperform by allowing to make a smarter choice of the step size parameter $\lambda$. A key thing to notice is that the convex combination of two feasible link flows solutions obtained via provides a feasible descent direction for the objective function of traffic equilibria problem (Eq.  <ref>). Thus, the value of $\lambda$ used for the convex combination of solutions is chosen to best improve the -objective. This solution strategy for -was originally implemented by <cit.> with a utility function dependent on travel time only and it is also integrated in the disaggregate simplicial decomposition () algorithm developed by <cit.>. An interesting feature of is the introduction of a column generation phase that dynamically update the consideration sets among O-D pairs. Consideration sets can be augmented with new paths according to different criterion such as the choice probabilities or the level of dissimilarity of the set of candidate paths respect to the existing paths among consideration sets.

§.§ Nonlinear least squares

Non-linear least squares () is a standard estimation method used in the and literature. In contrast to maximum likelihood estimation <cit.>, there is no need to make distributional assumptions about the data generating process <cit.> of the traffic counts and there is a variety of specialized algorithms to minimize the optimization objective <cit.>. Note that when the response function is linear, reduces to ordinary least squares () and thus, the solution could be proved to be unique under mild conditions. Under a non-linear response function, there may exist multiple local minima and saddle points that could make the optimization to be sensitive to the starting points for optimization and to converge toward a local but not the global minima. Therefore, the performance of the optimization algorithms heavily depends on the class of non-linearity of the response function respect to the parameters of interest.

§.§.§ Problem formulation

Consider the regression equation:

\begin{equation}
\label{eq:NLLS-regression-equation}
\vy = m(\vbeta,\mX) + \vu
\end{equation}

where $m(\cdot)$ is the response function, $\mX \in \sR^{|N| \times |K|}$ is the matrix of values of a set of exogenous attributes, $\vy \in \sR^{|N|}$ is the vector of values of the dependent variables,  $\vbeta \in \sR^{|K|}$ is the true parameter vector and $\rvu  \in \sR^{|N|}$ is a vector of random perturbations of some arbitrary distribution.

The problem consists in finding the estimator $\hat{\vbeta}_{NLLS}$ that minimizes the residuals, namely, the deviation between the vector of observed measurements $\bar{\vy}$ and the predictions of the response function:

\begin{equation}
\label{eq:NLLS-argmin-problem}
\hat{\vbeta}^{\star}
%= \arg \min_{\hat{\vbeta}} \ell(\hat{\vbeta})
= \arg \min_{\hat{\vbeta}} \|\bar{\vy} - m(\hat{\vbeta},\mX)\|_2^2 % =  S_T(\hat{\rvbeta})
\end{equation}

The application of the formulation to our  problem is direct. Let's be $\hat{\vtheta} \in \sR^{|D|}$ the vector of estimated coefficients in the travelers' utility function,  $\vt \in \sR^{|A^o|}$ the vector of link travel times, $\mZ \in \sR^{|A| \times |K_Z|}$ the matrix of exogenous link attributes and $\vx(\hat{\vtheta}, \mZ, \vt)$ the vector with the link flow functions of any link $a \in A^o$. In particular, when all attributes of the travelers' utility function are exogenous and known, such that $\vt = \bar{\vt}$, $\vx(\hat{\vtheta}, \mZ, \bar{\vt})$ becomes a vector valued function that resembles the response function $m(\vbeta,\mX)$ in and which has the following closed form:

\begin{align}
\label{eq:NLLS-link-flow-equation}
\vx(\hat{\vtheta}, \mZ, \bar{\vt})
= \mIx \vf
&= \mIx \left((\mIqT \vq) \circ  \vpf(\hat{\vtheta}, Z, \bar{\vt})  \right)
\end{align}

where $\vpf(\vtheta, \mZ, \bar{\vt})$ is the vector of path choice probabilities:

\begin{align}
\label{eq:NLLS-path-choice-probabilities-equation}
\vpf(\hat{\vtheta}, \mZ, \bar{\vt}) &= \exp\left(\displaystyle \mIxT \vv_{x}(\hat{\vtheta}, \mZ, \bar{\vt})     \right) \oslash \left(\displaystyle   \mIqT \mIq \exp(\mIxT \vv_{x}(\hat{\vtheta}, \mZ, \bar{\vt})  )\right)
\end{align}

Note that $\oslash$ is an operator for element wise division and $\vv_{x}(\vtheta, \mZ, \bar{\vt})$ is the vector of link utilities:

\begin{equation}
\label{eq:NLLS-link-utilities}
\vv_x (\hat{\vtheta}, \mZ, \bar{\vt})
%= \vt^{\star} \theta_t  +  \mZ \vtheta_Z
= \bar{\vt} \hat{\theta}_t  +   \mZ \hat{\vtheta}_Z
\begin{bmatrix}
\bar{\vt} & \mZ
\end{bmatrix}
\begin{bmatrix}
\hat{\theta}_t
\\
\hat{\vtheta}_Z
\end{bmatrix}
\begin{bmatrix}
\bar{\vt} & \mZ
\end{bmatrix}
\hat{\vtheta}
\end{equation}

Finally, if $\bar{\vx} \in \sR^{|A^{o}|}$ is the vector of observed traffic counts $\bar{x}$, the estimator $\hat{\vtheta}^{\star}$ of $\vtheta \in \sR^{|D|}$ can be obtained as follows:

\begin{equation}
\label{eq:NLLS-argmin-learning-problem}
\hat{\vtheta}^{\star} = \arg \min_{\hat{\vtheta}} \|\vx(\hat{\vtheta}, \mZ, \bar{\vt}) -\bar{\vx}\|_{2}
%= \arg \min_{\hat{\vtheta}} \left\|\mIx \left( (\mIqT \vq) \circ  \vpf(\hat{\vtheta})  \right) -\bar{\vx} \right\|_{2}
\end{equation}

§.§.§ Identifiability of utility function coefficients

The use of the model in our problem imposes rules of identifiability that are similar to those used in discrete choice models. For the computation of the choice probabilities, differences in utility among the alternatives is what matters <cit.>. Therefore, if an attribute of the paths' utility function is the same within each path set, then the coefficient weighting that attribute will be, by construction, not identifiable. From a behavioral point of view, travelers would perceive no difference in utilities when facing a path choice decision and thus, their choices will be no informative about the strength of their preference for that attribute. From an optimization standpoint, there will be no gradient respect to the attribute's coefficient because the choice probabilities are invariant to changes in the value of that coefficient.

A second rule of identifiability refers to the inclusion of alternative specific constants, which may lead to an over-specification of the utility function. The practical implication in our problem is that the specific constant in one of the paths connecting each O-D pair should be fixed, e.g. setting a constant to zero. A third rule of identifiability is related to the minimum amount of traffic counts measurements required to estimate a certain number of coefficients in the utility function, which can be referred as empirical identification <cit.>. Similar to the rank condition in , the number of observations should higher or equal than the number of parameters but, in practice, a larger amount of data may be needed due to the existence of collinearity between observations.

To show the relevance of the identification rules, let's consider the network in Figure <ref> and suppose that there are link flow measurements available for the two links. Assume the network is uncongested, $t_1(x_1) \approx t_1^0, \ t_2(x_2) \approx t_2^0$, and that the observed link counts perfectly matched the true link flow solution at equilibria, i.e. $(\bar{x}_1,\bar{x}_2)  = (x^{\star}_1,x^{\star}_2), \ \bar{x}_1 + \bar{x}_2 =\bar{q} $.
From Eq. <ref>, Section <ref>, we note that the link flow solution of the inner level problem can be written in closed form as:

\begin{equation}
\label{eq:sue-solution-small-uncongested-network}
x^{\star}_1 = \frac{\bar{q}}{1+\exp{(\theta_t(t^0_1-t^0_2) + \theta_c(c_1-c_2))}} \ , \  x^{\star}_2 = \frac{\bar{q}}{1+\exp{(\theta_t(t^0_2-t^0_1) + \theta_c(c_2-c_1))}}
\end{equation}

which satisfies both the conservation constraints between travel demand and path flows and the non-negativity constraints of link flows (Eq. <ref>, Section <ref>) for $\forall \theta_t \in \sR, \forall\theta_c \in \sR$. From the two link count measurements, we can derive a system of two equations which are dependent on the two unknowns coefficients of the utility function:

\begin{align}
\label{eq:inverse-sue-solution-uncongested-small-network}
\theta_t(t_1(\bar{x}_1)-t_2(\bar{x}_2)) + \theta_c(c_1-c_2) &= \ln\left(\frac{\bar{q}}{\bar{x}_1}-1\right) = \ln\left(\frac{\bar{x}_1 + \bar{x}_2}{\bar{x}_1}-1\right) = (\ln\left(\bar{x}_2\right)-\ln\left(\bar{x}_1\right)) \nonumber \\
%x^{\star}_1 = \frac{\bar{q}}{\exp{\theta_c(c_1-c_2))}}  \to   x^{\star}_2 = \frac{\bar{q}}{\exp{\theta_c(c_2-c_1))}}
\theta_t(t_2(\bar{x}_1)-t_1(\bar{x}_2)) + \theta_c(c_2-c_1) &= \ln\left(\frac{\bar{q}}{\bar{x}_2}-1\right) = \ln\left(\frac{\bar{x}_1 + \bar{x}_2}{\bar{x}_2}-1\right) = -(\ln\left(\bar{x}_1\right)-\ln\left(\bar{x}_2\right))
%x^{\star}_1 = \frac{\bar{q}}{\exp{\theta_c(c_1-c_2))}}  \to   x^{\star}_2 = \frac{\bar{q}}{\exp{\theta_c(c_2-c_1))}}
\end{align}

Note that in line with the first identification rule, Eq. <ref> can be solved for $\theta_c$ or $\theta_t$ if the difference of travel time or cost are zero, respectively (see Section <ref>). Besides, from the third rule of empirical identification, the two traffic count measurements would suffice to identify the two parameters. However, the two equations in  Eq. <ref> are identical, hence linearly dependent and thus, one equation is not providing additional information to identify another coefficient of the utility function.

The objective function of the optimization problem may be globally minimized in a point that is not attainable or where the gradient does not vanish. Suppose that the traffic counts are consistent with , such that  $(\bar{x}_1,\bar{x}_2) = (0,\bar{q})$ or $(\bar{x}_1,\bar{x}_2) = (\bar{q},0)$. With utility function dependent on travel time only, these link flow measurements are consistent with the limiting cases where $\theta_t \to -\infty$ or $\theta_t \to \infty$, respectively. While the identification rules are satisfied, the data generating process does not follow -. To our knowledge, testing for the -assumption remains an open question in the literature. If this assumption does not hold in practice, the optimization algorithm chosen to fit the travel time coefficient is expected to improve the objective function over iterations even when the global optima is not attainable.

Suppose that the monetary cost and the free flow travel times are equal in the two links, i.e. $c_1 = c_2, \ t^{0}_1 = t^{0}_2$. Under , by definition, the link travel times will be the same for any level level of demand $\bar{q}$. Note that the later will be true even if the free flow travel time of the links are different, provided that the demand level surpasses a certain threshold such that both links are utilized. Because the travelers are uniformly distributed among paths and the path utilities are the same, -will reproduce the same equilibria. In line with the first rule of identifiability, the travel time coefficient $\theta_t$ is not identifiable but regardless, any value of $\theta_t$ will perfectly reproduce the observed link flows. An optimization library will not necessarily warn about these identification problems and it may return an arbitrary estimate for $\theta_t$. A basic sanity check is to analyze if the objective function improves over iterations. Note, however, that this identifiability issue will not affect the goodness of fit but only statistical inference.

§.§.§ Solution methods

Second order optimization methods remain as the gold standard to estimate the travelers' utility function from traffic counts in the and literature. One of the earliest examples in the literature is found in <cit.>, who use the secant method, a Quasi-Newton method for unidimensional optimization, to estimate the coefficients of a utility function dependent on travel time only. More recent examples are <cit.> and <cit.> who solved the outer level problem of the via Sequential Quadratic Programming (SQP), a variant of the Newton method for constrained optimization. Applications of second order methods are also found in travel behavior research where Broyden–Fletcher–Goldfarb–Shannon () remains as a preferred optimizer to find the maximum likelihood estimates in discrete choice models <cit.>.

Second order optimization methods can achieve a superquadratic or superlinear convergence in convex problems. They are appealing when datasets are of moderate size and thus, when the computational cost for the calculation or approximation of the Hessian matrix is reasonable. However, their convergence guarantees heavily relies on how well the curvature of the objective function informs about the direction of its steepest descent/ascent. Besides, in non-convex problems, their performance strongly depends on how close the initial points for optimization are to a local optima and on whether these points are located within a locally convex region of the objective function. Furthermore, the frequent sign changes of the curvature in non-convex functions may cause that second order optimization methods does not converge to a local minima or that they get stuck at saddle points or flat regions of the optimization landscape.

First order methods have become a gold standard for non-convex optimization due in part to the huge computational gains that are attained when avoiding Hessian matrix computation. The machine learning community has made continued efforts to develop multiple variants of these methods that are able to accelerate the optimization and without compromising computational costs and that can perform reasonably well in highly non-convex optimization landscapes <cit.>. Some optimizers that remains in state-of-the-art applications include stochastic gradient descent () <cit.> and the Adagrad <cit.> and Adam <cit.> optimizers.

In the context of quasiconvex optimization problems, <cit.> proposed a modified version of gradient descent () in which the gradient direction is normalized by its norm.  The existence of flat regions in quasiconvex problems cause that the gradients become small despite that they may be pointing to the right direction of improvement. In this setting, the normalization of the gradient is helpful to keep improving the objective function in flat regions as well as to avoid gradient explosion in sharp regions of the feasible space. The application of normalized gradient descent () in unidimensional unconstrained optimization problems generates parameter updates according to the sign of the first derivative of the objective function.

Interestingly, the and literature has not explored the use of first order optimization methods to estimate the utility function coefficients using system level data. Further sections of this paper will show some theoretical and empirical results to support the choice of over the standard second order optimization methods in previous work.

§.§ Bilevel optimization

§.§.§ Problem formulation

The sole application of the method to estimate the travelers' utility function coefficients from network level data requires that travel times are known and exogenous. In uncongested networks, this assumption may be enforced by setting the travel times to be equal to the free flow travel times. In congested networks, this assumption is not plausible and hence, it is necessary to solve both a and a -problems of a bilevel formulation in an iterative fashion. Below is the standard formulation of the problem solved in the literature and that we extend to account for travelers' utility functions with multiple coefficients:

\begin{equation}
\label{eq:bilevel-formulation-learning-problem}
\begin{aligned}
%\text{Outer problem : }
\min \limits_{\hat{\vtheta}} \  & \ell(\hat{\vtheta}) = \| \vx(\hat{\vtheta}) - \bar{\vx}\|_2^2 \\
\st & \ \vx(\hat{\vtheta}) \in \arg\max_{\vx} g(\hat{\vtheta}, \vx,\vf) \\
&%\vq \in
\begin{aligned}[t]
%\text{Inner problem : }
\max_{\vx,\vf} \ & g(\hat{\vtheta}, \vx,\vf)
%= \  \displaystyle \sum_{a \in A}\int_{0}^{x_a}  {v}^{(a)}_{x}(u, \vtheta) du -
= \displaystyle \   -\left\langle \vf, \ln \vf\right\rangle +  \sum_{a \in A} x_a \sum_{k \in K_{\mZ}} \theta_k \cdot Z_{ak} +   \sum_{a \in A}\int_{0}^{x_a}  \theta_t t_a(u) \ du \\
\st & \ \mIq\vf = \vq \\
& \ \mIx\vf  = \vx\\
& \ \vx,\vf \geq \vzero
\end{aligned}\\
%&&&\vx \in [-10, 10]^2,\; \mathbf{y} \in [-10, 10]^2
\end{aligned}
\end{equation}

where $\ell(\hat{\vtheta})$ and $g(\hat{\vtheta}, \mZ, \vx,\vf)$ are the objective functions at the upper and lower level and $\hat{\vtheta} \in \mathbb{R}^{|D|}$ is the vector of estimated utility function coefficients. The idea behind expressing the objective functions in terms of $\hat{\vtheta}$ can be understood as follows. At the inner level, any value of $\hat{\vtheta}$ will lead to a different assignment of path/link flows in the network, hence, the objective function $g(\hat{\vtheta}, \vx,\vf)$ depends on $\hat{\vtheta}$. Note also that, thanks to the -property (Eq. <ref>), $\ell(\hat{\vtheta})$ can be approximated as an analytical function of $\hat{\vtheta}$.

§.§.§ Solution methods

A brute-force strategy to solve the bilevel formulation in Problem <ref> would consist of performing a grid search on all feasible values of $\hat{\vtheta}$. Therefore, for each value of $\hat{\vtheta}$, -would be first computed at the inner level and then the objective function at the upper level would be evaluated using the closed form expression derived in Eq. <ref>, Section <ref>. Lastly, a candidate solution would be the value of $\theta$ that minimizes $\ell(\hat{\vtheta})$. Note however that this solution strategy becomes intractable in a higher dimensional space resulting from a multi-attribute travelers' utility function or from the joint estimation of the utility function coefficients and an O-D matrix.

A standard heuristic employed in the and literature to solve bilevel formulations is an alternating optimization of the inner and outer level problem. Over iterations, this heuristic is expected to converge toward a stationary point that, at the outer level, minimizes the gap between observed and predicting flows and, at the inner level, induces a link/path flow solution that satisfies traffic equilibria. The solution strategy in our problem boils down to performing gradient-based updates of the utility function coefficients at the outer level and to compute -at the inner level using the coefficients $\hat{\vtheta}$ previously updated at the outer level. The resulting travel times from the -computation are then feed to the outer level and the process is repeated until some criterion is convergence is met. There are many variants of solution methods in the literature but most can be summarized with the following steps:

* Initialization of estimated parameters (e.g. utility function coefficients or O-D matrix)

* Inner level problem: solve -and update travel times in the network

* Outer level problem: solve with new travel times and update estimated parameters

* Repeat 2 and 3 until fulfilling some convergence criteria

In our application and similar to the literature, the utility function coefficients in $\hat{\vtheta}$ are the only free parameters.

A key advantage of -in network modeling applications is the existence of a closed form and unique solution at path and link flow spaces. This property can significantly ease the solution of the typical bilevel optimization structure that arises in problems. During the alternating optimization, the closed form solution is valid but only for the set of travel times obtained in the inner level problem. Thus, it is convenient to perform small updates of $\hat{\vtheta}$ at the outer level problem such that the inner level solution remains valid.

§ MATHEMATICAL PROPERTIES OF THE OPTIMIZATION PROBLEM

This section studies the mathematical properties of our optimization problem. Note that our problem falls into the class and since is a subclass of , the analysis of mathematical properties of our problem is also relevant for the problem.

§.§ Exogenous case

To ease the analysis of the mathematical properties, it is convenient to start assuming that the utility function is dependent on exogenous attributes only. This assumption significantly facilitates the solution of the bi-level formulation because the inner level problem  at each iteration consists of a single pass of . This assumption is plausible in an uncongested network where travel times are close to the free flow travel times or when there is access to travel time measurements for every link in the network. Alternatively, this assumption may be enforced by either setting the travel time coefficient of the utility function or the parameter $\alpha$ of the cost performance to zero but this is arguably less realistic.

§.§.§ Non-convexity

The non-convexity of the and problems respect to the travelers' utility function coefficients has been discussed in prior literature but via counterexamples only <cit.>. To our knowledge, no formal proof of the non-convexity has been given for neither and . A key observation to prove the non-convexity of these problems is that their objective function $f$ is upper-bounded respect to the utility function coefficients. In a one-dimensional case where $\theta \in \mathbb{R}$, the upper-bounds of $f$ corresponds to the solutions of when $\theta \to \infty^+$ and  $\theta \to \infty^-$. To extend this intuition to a multi-dimensional case where $\vtheta \in \mathbb{R}^{|D|}$, it is useful to introduce a general definition of convexity <cit.>.

A function $f$ is convex iff $\forall y_1, y_2 \in S, \lambda \in [0,1]$:

f(\lambda \vy_1 + (1-\lambda)\vy_2) \leq \lambda f(\vy_1) + (1-\lambda)f(\vy_2)

The objective function of the problem is lower-bounded and upper-bounded

The objective function $\ell: \mathbb{R}^{|D|} \to \mathbb{R}$ of the problem can be expressed as:

\ell(\vtheta) = \|\vx(\theta) -\bar{x}\|_2^2
= \sum_{i \in N} \left(x_i(\vtheta)- \bar{x}_i\right)^2
%= \sum_{n \in N} \left(\sum_{i \in A_{rs}} q_i \dfrac{\displaystyle \exp(V_i(\vtheta))}{\displaystyle \sum_{j \in A_{rs}}  \exp(V_j(\vtheta)) }\delta_{in}
%- \bar{x}_n\right)^2

where $N$ is the set of observed traffic count measurements. Each traffic function $x_i, \forall i \in N$ can be bounded as:

\begin{align}
\label{eq:bound-traffic-count-function}
0 &\leq x_i(\theta) \leq \sum_{w \in W} q_w = Q
\end{align}

where $Q \in \sR_+$ is the total demand, namely, the sum of all cells in the O-D matrix $\mQ \in \sR^{V\times V}$. Then, lower and upper bounds of the objective function can be found as follows:

\begin{align}
\label{eq:bound-loss-function}
0 &\leq x^2_i(\theta) \leq Q^2 \nonumber \\
-2\bar{x}_i +\bar{x}_i^2 &\leq x^2_i(\theta) -2\bar{x}_i +\bar{x}_i^2\leq Q^2 -2\bar{x}_i+\bar{x}_i^2 \nonumber  \\
\bar{x}_i(\bar{x}_i-2)&\leq (x_i(\theta) -\bar{x}_i)^2\leq (Q -\bar{x}_i)^2 \nonumber  \\
\sum_{i \in N} \bar{x}_i(\bar{x}_i-2)  &\leq \sum_{i \in N} (x_i(\theta) -\bar{x}_i)^2 \leq \sum_{i \in N} (Q -\bar{x}_i)^2
\nonumber \\
(\bar{\vx}^\top-2\cdot\vone^\top)\bar{\vx}  &\leq \|\vx(\theta) -\bar{\vx}\|_2^2 \leq \|\vone^\top \vq -\bar{\vx}\|_2^2
\end{align}

which completes the proof

Let's now prove the non-convexity of the problem with the following proposition:

Suppose the coefficients of the travelers' utility function in the problem are identifiable. Then, the problem is not convex if its objective function is upper-bounded

Let's prove this by contradiction. Assume that the objective function  $\ell(\vtheta) = \|\vx(\vtheta) -\bar{\vx}\|_2^2, \ \forall \vtheta \in \mathbb{R}^{|D|}$ is convex. Then, let's use Definition <ref> of convexity at points $\vy_1 = \frac{\vtheta_1 - (1-\lambda)\vtheta_2}{\lambda}$ and $\vy_2 = \vtheta_2$, with $\vy_1, \vy_2,\vtheta_1,\vtheta_2 \in \mathbb{R}^{|D|}$ and $\lambda \in ]0,1]$:

\begin{align*}
\ell(\lambda \vy_1 + (1-\lambda)\vy_2) &\leq \lambda \ell(\vy_1) + (1-\lambda)\ell(\vy_2)\\
\ell\left(\lambda \left(\frac{\vtheta_1 - (1-\lambda)\vtheta_2}{\lambda}\right) + (1-\lambda)\vtheta_2\right) &\leq \lambda \ell\left(\frac{\vtheta_1 - (1-\lambda)\vtheta_2}{\lambda}\right) + (1-\lambda)\ell(\vtheta_2)\\
\ell(\vtheta_1)  &\leq \lambda \ell\left(\frac{\vtheta_1 - (1-\lambda)\vtheta_2}{\lambda}\right) + (1-\lambda)\ell(\vtheta_2)
\end{align*}

which implies that:

\begin{equation}
\label{eq:non-convexity-proof1}
\ell\left(\frac{\vtheta_1 - (1-\lambda)\vtheta_2}{\lambda}\right) \geq \frac{\ell(\vtheta_1)-(1-\lambda)\ell(\vtheta_2)}{\lambda}= \frac{\ell(\vtheta_1)-\ell(\vtheta_2)}{\lambda} + \ell(\vtheta_2)
\end{equation}

If $\vtheta$ is identifiable, $f$ is not a constant function and $\exists \lambda \in [0,1]: \ell(\vtheta_1) > \ell(\vtheta_2)$. From Eq. <ref>, when $\ell(\vtheta_1) > \ell(\vtheta_2)$, the RHS grows with no bound as $\lambda \to 0^{+}$. Hence, the function $f$ is not upper bounded, which generates a contradiction and it completes the proof.

Convexity  is not a necessary but a sufficient condition for global optimality. Therefore, before ruling out the global optimality of our problem, it is convenient to study a generalization of convexity known as a pseudo-convexity that imposes weaker conditions on the optimization problem and that provides sufficient conditions for global optimality.

§.§.§ Pseudo-convexity

Pseudo-convex functions are a subclass of the quasi-convex functions. Hence, the properties of pseudo-convex functions are also valid for quasi-convex functions, while the converse is not true <cit.>. Pseudo-convexity is a generalization of convexity that, as quasi-convexity, extends the notion of unimodality to higher dimensional spaces. In contrast to quasi-convex functions but similar to convex functions, the first order necessary optimality condition of pseudo-convex functions suffices for global optimality <cit.>. This property will be key later to prove the existence of a global minimizer in our optimization problem. A formal definition of pseudo-convexity is the following <cit.>:

$f$ is a pseudo-convex function on the feasible set $S$ iff $\forall x_1,x_2 \in S$:

\begin{equation}
f(x_2) < f(x_1) \implies \nabla f(x_1)^\top(x_2-x_1) < 0
\label{eq:pseudo-convexity-definition-1a}
\end{equation}

or equivalently,

\begin{equation}
\label{eq:pseudo-convexity-definition-1b}
\nabla f(x_1)^\top(x_2-x_1) \geq 0 \implies  f(x_2) \geq f(x_1)
\end{equation}

The definition of pseudo-convexity can be assessed via the augmented hessian of the objective function <cit.>. Formally, $f$ is pseudoconvex if there exists a scalar $v, \  0 \leq v < \infty$ such that its augmented Hessian $H(x) + v \nabla f(x) \nabla f(x)^T$ is positive semidefinite for all $x \in X$ <cit.>. This proposition can be also evaluated in terms of the signs of the principal minors of the bordered Hessian <cit.>.

For the remaining analyses of mathematical properties, it will be also helpful to introduce Definition <ref> of coordinate-wise pseudo-convexity. A function $f$ is coordinate-wise pseudo-convex if, when all coefficients of the utility function except for one coefficient are fixed, it is pseudo-convex respect to the non-constant coefficient. Note that if the function is dependent on a single attribute, the definitions of coordinate-wise pseudo-convexity and pseudo-convexity are equivalent. Formally:

$\ell: \sR^{|D|} \to \sR$ is a coordinate-wise pseudo-convex function iff $\forall i \in D, \forall \vtheta^1,\vtheta^2 \in S_i$, such that $S_i = \{\vtheta^1,\vtheta^2 \in \mathbb{R}^{|D|}\ | \ \theta_j^1 = \theta_j^2, \forall j \neq i \}$:

\begin{equation}
\ell(\vtheta^2) < \ell(\vtheta^1) \implies \dfrac{\partial \ell(\vtheta^1)}{\partial \theta_i}(\theta_i^2-\theta_i^1) < 0
\label{eq:coordinate-wise-pseudo-convexity-definition-1a}
\end{equation}

or equivalently,

\begin{equation}
\label{eq:coordinate-wise-pseudo-convexity-definition-1b}
\dfrac{\partial \ell(\vtheta^1)}{\partial \theta_i}(\theta_i^2-\theta_i^1) \geq 0 \implies  \ell(\vtheta^2) \geq \ell(\vtheta^1)
\end{equation}

Proposition <ref>, <ref> shows that coordinate-wise monotonicity of the traffic count (response) functions implies coordinate-wise pseudo-convexity of the objective function of the problem under exogenous travel times. The definition of coordinate-wise monotonicity is the following:

The traffic count (response) functions are monotonic respect to each coefficient of the utility function if, when all coefficients of the utility function except for one coefficient are kept constant, the response functions are monotonic respect to the non-constant coefficient.

Propositions <ref> and <ref>, <ref> illustrate two cases where under mild assumption the coordinate-wise monotonicity of the traffic flow functions holds in a general transportation network. While it is easy to create counter examples where the coordinate-wise monotonicity of the traffic flow functions is violated, analysis on synthetic data suggests that the assumption holds in practice (see Section <ref>).

§.§.§ Existence and uniqueness of global minima

In absence of measurement error in the traffic counts, the proofs for existence and uniqueness of a global optimality are direct (Propositions <ref> and <ref>, <ref>). While the absence of measurement error can be enforced in experiments with synthetic data, it is implausible in real world problems where traffic count data is expected to be noisy. Fortunately, the pseudo-convexity of the objective function in our problem can be leveraged to extend the proofs of existence and uniqueness in presence of measurement error. A key property of pseudo-convex functions that also holds in convex functions is that every local minima is also global minima <cit.>. This relaxes the positive (semi)definite assumption of the Hessian of the objective function used to prove (strict) global optimality in convex problems. Following the proof in <cit.> for multi-variable pseudo-convex functions, Proposition <ref> shows that this property can be also applied to prove the existence of global minima in our problem.

Assume the objective function $f$ of the problem under exogenous travel times is coordinate-wise pseudo-convex. Then, if the gradient of $f$ vanishes at $\vtheta^{\star} \in \mathbb{R}^{|D|}$, $\vtheta^{\star}$ is a global minima.

By assumption, $\grad_{\vtheta} \ell (\vtheta^{\star}) = 0$, which implies that $\dfrac{\partial \ell(\vtheta^{\star})}{\partial \theta_i} = 0, \forall i \in D$. By applying Definition <ref> of coordinate-wise pseudo-convexity to each coordinate $i \in D$, $\ell(\vtheta^2) \geq \ell(\vtheta^{\star}), \ \forall \vtheta^2 \in \mathbb{R}^{|D|}$, which completes the proof.

Proposition <ref>, <ref> gives a sufficient condition for the existing of a vanishing gradient (derivative) in the unidimensional case. While this seems restrictive, it addresses the standard setting studied in prior literature with a utility function dependent on travel time only.

Proposition <ref> proves the uniqueness of the global optima in our problem by leveraging the the strict quasi-convexity of differentiable pseudo-convex functions proven by <cit.>. A formal definition of strict quasiconvexity is the following <cit.>:

f is said to be strictly quasi-convex if for every $x_1, x_2 \in S, x_1 \neq x_2, \lambda \in [0,1]$:

f(x_1) < f(x_2) \implies f(\lambda x_1 + (1-\lambda_1)x_2) < f(x_1)

or equivalently:

f(\lambda x_1 + (1-\lambda) x_2)< \max\left(f(x_1), f(x_2)\right)

Analogously to Definition <ref>, the definition of strict quasi-convexity can be also applied coordinate-wise:

$\ell: \sR^{|D|} \to \sR$ is said to be coordinate-wise strictly quasi-convex iff $\forall i \in D, \forall \vtheta^1,\vtheta^2 \in S_i$, such that $S_i = \{\vtheta^1,\vtheta^2 \in \mathbb{R}^{|D|}\ | \ \theta_j^1 = \theta_j^2, \forall j \neq i \}$:

\ell(\vtheta_1) < \ell(\vtheta_2) \implies \ell(\lambda \vtheta_1 + (1-\lambda_1)\vtheta_2) < \ell(\vtheta_1)

or equivalently:

\ell(\lambda \vtheta_1 + (1-\lambda) \vtheta_2)< \max\left(\ell(\vtheta_1), \ell(\vtheta_2)\right)

If the objective function of problem is (coordinate-wise) pseudo-convex and a global minima exist, then the global minima is unique.

Let's prove this by contradiction. Let's be $\ell: \mathbb{R}^{|D|} \to \mathbb{R}$ the objective function of the problem and $\vtheta \in \mathbb{R}^{|D|}$ the vector of coefficients of the utility function. Suppose that there are two global minima $\vtheta^{\star}_1,\vtheta^{\star}_2 \in \mathbb{R}^{|D|}: \ell(\vtheta^1) = \ell(\vtheta^{\star}_2) = \ell^{\star}$ and $\theta^{\star}_1 \neq \theta^{\star}_2$. By the multivariate mean value theorem (MMVT), for $\vtheta^{\star}_3 \in \mathbb{R}^{|D|}$ , $\exists \theta^{\star}_3 \in (1-\lambda)\vtheta^{\star}_1 + \lambda \vtheta^{\star}_2 : \grad \ell (\vtheta = \vtheta^{\star}_3) = 0, \ \lambda \in [0,1]$ and by Proposition <ref>, $\ell(\theta^{\star}_3) = f^{\star}$. By using MMVT recursively, $\forall \vtheta^{\star}_i \in [\theta^{\star}_1,\theta^{\star}_2], \ \ell(\theta^{\star}_i) = \ell^{\star}, \ \theta^{\star}_i = (1-\lambda)\vtheta^{\star}_1 + \lambda \vtheta^{\star}_2, \ \lambda \in [0,1]$ and hence, $ [\vtheta^{\star}_1,\vtheta^{\star}_2]$ forms a convex set. By Property 2, <cit.>, if $\ell$ is pseudo-convex, $\ell$ is also strictly quasi-convex. Similarly, if $\ell$ is coordinate-wise pseudo-convex, $\ell$ is also coordinate-wise strictly quasi-convex. By applying Definition <ref> of strict quasi-convexity on $\vtheta^1, \vtheta^2 \in \mathbb{R}^{|D|}$ gives $\ell^{\star} = \ell(\lambda \theta^1 + (1-\lambda) \theta^2)< \max\left(\ell(\theta^1), \ell(\theta^2)\right) = \ell^{\star}$, which leads to a contradiction and it completes the proof.

§.§ Endogenous case

In a congested network, travel times are endogenous and this requires to solve a bilevel formulation of the problem (see Section <ref>). Note that even if the inner and outer level problems of the bilevel formulation were convex, the analysis of global optimality is complex <cit.>. Because convex functions are a subclass of the pseudo-convex functions, the pseudo-convexity of the outer problem is expected to make the analysis even harder. Despite the above difficulties, this section proves two properties of the inner level and outer level problems of the bilevel formulation that can improve the convergence of an alternating optimization algorithm.

§.§.§ Inner level problem

The -problem with a utility function dependent on travel time only is known to be strictly convex in path (and link) flow space under link performance functions that are monotonically increasing and dependent on the traffic flow on the links only. To our knowledge, Proposition <ref>
extends, for the first time, the aforementioned property to multi-attribute utility functions.

Suppose that the travelers' utility function is linear-in-parameters and that the link cost functions are monotonically increasing and dependent on the traffic flow on the link only. Then, the -problem has a unique solution in path flow space.

Let's transform the objective function of Problem <ref>, Section <ref> into a minimization problem with the following argument:

\begin{align*}
g(\vx,\vf, \vtheta)
%= \sum_{a \in A}\int_{0}^{x_a} t_a(w)dw  - \frac{1}{\theta_t} \left\langle \vf, \ln \vf\right\rangle
&= -\sum_{a \in A}\int_{0}^{x_a} v_a(u,\vtheta)du - \left\langle \vf, \ln \vf\right\rangle
= -\sum_{a \in A}\int_{0}^{x_a} \left(\theta_t t_a(u) - \sum_{k \in K_{\mZ}} \theta_k \cdot Z_{ak} \right)du + \left\langle \vf, \ln \vf\right\rangle\\
&= g_1(\vx, \vtheta) +g_2(\vf)
\end{align*}

From Theorem 1 in <cit.>, it follows that $g_2(\vf) = -\left\langle \vf, \ln \vf\right\rangle$ is strictly convex on the feasible set $S$ defined by the constraint in Problem <ref>, Section <ref>. Note that the matrix $\mZ \in \mathbb{R}^{A \times k_{\mZ}}$ with link attributes values is exogenous and thus, $g_1(\vx, \vtheta)$ can be rewritten as:

g_1(\vx, \vtheta) = -\theta_t \sum_{a \in A}\int_{0}^{x_a} t_a(u)  du -  \sum_{a \in A} x_a \left(\sum_{k \in K_{\mZ}} \theta_k \cdot Z_{ak}\right)

To analyze the convexity of the first term $g_1(\vx, \vtheta)$, it suffices to analyze its second derivative respect to any path flow variable $f_h$:

\begin{align*}
\dfrac{\partial g_1(\vx,\vf, \vtheta)}{\partial  f_h}
&= -\theta_t \dfrac{\partial }{\partial  f_h} \left(\sum_{a \in A}\int_{0}^{x_a} t_a(u)  du\right)-  \sum_{a \in A} \dfrac{\partial x_a}{\partial  f_h} \left(\sum_{k \in K_{\mZ}} \theta_k \cdot Z_{ak}\right)
%&= \theta_t \dfrac{\partial }{\partial  f_h} \left(\sum_{a \in A}\int_{0}^{x_a} t_a(u)  du\right)+  \sum_{a \in A} \left( \left(\sum_{k \in K_{rs}} q_{rs} \delta_{ak}\right) \left(\sum_{k \in K_{\mZ}} \theta_k \cdot Z_{ak}\right)\right)
%&= \theta_t \dfrac{\partial }{\partial  f_h} \left(\sum_{a \in A} (T_a(x_a,\vtheta) -T_a(0,\vtheta) )\right)+  \sum_{a \in A} \dfrac{\partial x_a}{\partial  f_h} \left(\sum_{k \in K_{\mZ}} \theta_k \cdot Z_{ak}\right)\\
%&= \sum_{a \in A} v_a(x_a,\vtheta) \dfrac{\partial x_a}{\partial  f_h}
%&= \sum_{a \in A} \left(v_a(x_a,\vtheta) \sum_{k \in K_{rs}} q_{rs} \delta_{ak}\right)
\end{align*}

Then, replacing the link flow variables by using the balance constraint for path and link flows in the second term of $g_1$:

\dfrac{\partial^2 g_1}{\partial^2 f_h}
%= \sum_{a \in A} \left(\dfrac{\partial v_a(x_a,\vtheta)}{\partial  f_h}   \sum_{k \in K_{rs}} q_{rs} \delta_{ak}\right)
= -\theta_t \dfrac{\partial^2 }{\partial^2  f_h} \left(\sum_{a \in A}\int_{0}^{x_a} t_a(u)  du\right)-  \sum_{a \in A} \cancelto{0}{\dfrac{\partial }{\partial  f_h}\left(\sum_{k \in K_{rs}} q_{rs} \delta_{ak}\right)}\left(\sum_{k \in K_{\mZ}} \theta_k \cdot Z_{ak}\right)
= -\theta_t \dfrac{\partial^2 }{\partial^2  f_h} \left(\sum_{a \in A}\int_{0}^{x_a} t_a(u)  du\right)

Given the class of link cost functions used in this proposition (Assumption <ref>, Section <ref>) and the fact that $\theta_t < 0$, it follows that $g_1$ is convex on the path flows variables $f_h$ <cit.>. As a consequence, $g(\vx,\vf, \vtheta)$ is a sum of a convex and strictly convex function, and thus, it is strictly convex. Also, all constraints in Problem <ref>, Section <ref> are affine and thus they define a convex set. Thus, the solution of problem <ref> has a unique global minima, which completes the proof.

§.§.§ Outer level problem

A pseudo-convex function can be convex and concave on its feasible space. Despite of the above, its gradient always point out to the global descent direction <cit.>. Proposition <ref> shows that this property also holds in the outer level objective of the problem.

Assume that the objective function of the outer level of the problem is coordinate-wise pseudo-convex and that a global optimal solution exist. Then, the negative gradient of the objective function points out to the global descent direction

Let's be $\vtheta^{\star} \in \mathbb{R}^{|D|}$ the global optima of the problem and $\bar{\vtheta} \in \mathbb{R}^{|D|}$ a feasible point that is not a global optima. By the coordinate-wise pseudo-convexity of $\ell: \sR^{|D|} \to \sR$:

\begin{equation}
\label{eq:coordinate-wise-pseudoconvexity-congested-case}
\ell(\vtheta^{\star}) < \ell(\bar{\vtheta}) \implies \dfrac{\partial \ell(\bar{\vtheta})}{\partial \theta_i}(\theta_i^{\star}-\theta_i) < 0, \quad \forall i \in D
\end{equation}

Since $\ell(\bar{\vtheta}) < \ell(\vtheta^{\star})$, $\dfrac{\partial \ell(\bar{\vtheta})}{\partial \theta_i}(\theta_i^{\star}-\bar{\theta}_i) < 0, \forall i \in D$. Then, if $\bar{\theta}_i < \theta_i^{\star}$, $\dfrac{\partial \ell(\bar{\vtheta})}{\partial \theta_i} >0$, and if $\bar{\theta}_i > \theta_i^{\star}$, $\dfrac{\partial \ell(\bar{\vtheta})}{\partial \theta_i} <0$. Therefore, if $\bar{\vtheta} - \vtheta^{\star} < 0, \  \minus\grad_{\theta} \ell(\bar{\vtheta}) < 0$ and if $\bar{\vtheta} - \vtheta^{\star}  > 0,\ \minus\grad_{\theta} \ell(\bar{\vtheta}) > 0$, which proves that the negative gradient is always pointing to the global descent direction.

In a congested network, each value of $\vtheta \in \mathbb{R}^{|D|}$ generates a different traffic assignment at the inner level of the bilevel formulation. Therefore, link travel times will change between iterations of an alternating optimization algorithm. As a consequence, the objective function of the outer level problem will be (coordinate-wise) pseudo-convex respect to the travelers' utility function coefficients only in a small neighborhood where travel times remain approximately constant.

§ SOLUTION ALGORITHM AND CONVERGENCE GUARANTEES

§.§ Inner level optimization

The solution method for the inner level problem of the bilevel formulation is described in Algorithm <ref>, <ref>. The stochastic network loading () is performed according to Algorithm <ref>, <ref>. The column generation phase is adapted from <cit.> and it is performed once per each iteration of the alternating optimization algorithm (Step 1, Algorithm <ref>, <ref>). Note that in the original implementation of , the column generation phase is performed at each iteration of -but in the context of a bilevel optimization, this path set augmentation strategy would not select good paths if the current solution for $\vtheta$ is far from the global optima. The best convex combination of solutions at each iteration of -(Step c, Algorithm <ref>, <ref>) is searched over a uniform grid in $[0,1]$ and with an arbitrary granularity depending on the desired level of accuracy. To reduce computational cost, the size of the augmented path set is constrained to a fix quantity ($k_g$) and the path set augmentation is only performed in a proportion ($\rho_{W}$) of the O-D pairs at each iteration. The sample of O-D pairs is selected according to the level of demand. To capture the correlation between paths in the consideration set due to overlapped link segments, we correct path utilities using the path size logit (PSL) factor <cit.>. The utility term associated to the PSL factor is weighted by a coefficient $\beta \in \sR$ that can be estimated <cit.>. For the scope of this paper, the coefficient $\beta$ is treated as an hyperparameter with default value equal to one. The utility term associated to the path size correction is updated before generating new paths and after selecting paths depending if the composition of the path set change when performing the steps 1 and 3 of the inner level optimization (Algorithm <ref>, <ref>),

§.§ Outer level optimization

To study the convenience of using first or second order methods, the outer level objective function of the bilevel formulation is optimized with an hybrid algorithm (see Algorithm <ref>, <ref>). In a non-refined stage, the algorithm uses Normalized Gradient Descent () (Algorithm <ref>,<ref>). Subsequently, the best solution obtained from the non-refined stage is used a starting point for the refined stage, where second order optimization methods are used to obtain a more accurate solution (Algorithm <ref>, <ref>). Two specialized second order methods were initially considered. The first is Gauss-Newton (), one of the standard optimizer in the Econometrics literature on non-linear regression <cit.> and which can be seen as an unconstrained version of the sequential quadratic programming () method. The second is the Levenberg–Marquardt () algorithm which, by interpolating between and gradient descent, can increase the robustness of in flat regions of the optimization landscape <cit.>. Given that is a particular case of the method (Algorithm <ref>, <ref>), only is used in the refined stage.

§.§ Bilevel optimization

Algorithm <ref>, Section <ref> shows the pseudocode for the alternating optimization of the inner and outer level of our bilevel formulation (Problem <ref>, Section <ref>). The algorithm returns the vector of utility function coefficients $\hat{\vtheta} \in \mathbb{R}^{|D|}$ that minimizes the discrepancy between predicted and observed link flows. To guarantee that the predicted link flows follow -, the last iteration of the alternating algorithm only solves the inner level problem using the value of $\hat{\vtheta}$ obtained in the previous iteration of the outer level problem.

Iterations $I$, initial vector of utility function coefficients $\vtheta_0$, inputs of  $T, \mZ, \vt_f, \vq, \vlambda_g,\rho_W$, inputs of  $T_1, T_2,\eta,\bar{\vx}, \vx, \vpf, \texttt{refined-method}, \eta_1, \eta_2, b, \delta$

Step 0: Initialization:

\vtheta:= \vtheta_0

Step 1: Alternating optimization

$i = 1 \ldots I$

$\vx, \vpf, \bar{\vt} \gets$ ($\vtheta,T, \mZ, \vt_f, \vq, \vlambda_g,\rho_W$)

$i < I$
$\vtheta \gets$ ($\vx, \vpf,\bar{\vx},\vtheta, T_1, T_2, \texttt{refined-method}, \eta_1, \eta_2, b, \delta$)

$\vtheta^{\star} = \argmin_{\{\vtheta_0,\ldots,\vtheta_{I-1} \}} \ell_i(\vtheta_i):= \|\bar{\vx}-\vx(\vtheta_i)\|^2$

§.§ Convergence guarantees

Our bilevel formulation can be casted as a mathematical program with equilibrium constraints (). For this class of problems, the feasible region is typically non-convex <cit.> and this makes difficult to find theoretical guarantees even for convergence toward local optima <cit.>. Thus, we found convenient to analyze a case where all attributes of the utility function are exogenous (Section <ref>) and where we could leverage the pseudo-convexity of the optimization problem (Section <ref>).

§.§.§ Theoretical $\epsilon$-convergence under exogenous travel times

The analysis of mathematical properties in Section <ref> provided sufficient conditions for the existence and uniqueness of a global optima. Now it remains to prove that the global optima is attainable when running normalized gradient descent () on a finite number of steps. We will consider the problem with a utility function dependent of a single exogenous attribute and where only the no-refined stage of the outer level algorithm is conducted. While this seems restrictive, it addresses the standard setting studied in prior literature with a utility function dependent on travel time only.

<cit.> proved the convergence of to a global minimum for functions that are both strictly quasi-convex and locally-Lipschitz. These results are relevant for our problem, since pseudo-convex functions are also quasi-convex  <cit.>. Locally-Lipschitz functions have the property of having its first derivative bounded by an arbitrary positive constant and they are required to be Lipschitz in a small region around the optimum. Generally, the Lipschitz constant is relevant to study the theoretical convergence speed of first optimization methods. Formally:

A unidimensional function $f: \mathbb{R} \to \mathbb{R}$ is Lipschitz iff $\ \forall x,y \in \mathbb{R}, \ G>0$:

\begin{align}
\label{eq:lipschitzness-uncongested-network}
|f(x)-f(y)| &\leq G\|x-y\| %\\
%|f(x^{\star})-f(x^{\star}+\epsilon)| &\leq G\|x^{\star}-(x^{\star}+\epsilon)\| \\
\end{align}

Suppose that the travelers' utility function in the problem depends on a single exogenous attribute. Then, the objective function is Lipschitz

Let's define the the objective function of the problem a $\ell: \sR^{|D|} \to \sR$. Starting from the LHS of Eq. <ref> and using the reverse triangle inequality:

\begin{align}
\label{eq:proof-lipschitzness-uncongested-network-1}
|\ell(\theta^1_t)-\ell(\theta^2_t)| \nonumber
&= \Big|\|\vx(\theta^1_t) -\bar{x}\|_2^2 -\|\vx(\theta^2_t) -\bar{x}\|_2^2\Big|\\ \nonumber
&\leq \|(\vx(\theta^1_t) -\bar{x}) - (\vx(\theta^2_t) -\bar{x})\|_2^2\\
&= \|(\vx(\theta^1_t) - \vx(\theta^2_t)\|_2^2
\end{align}

We can now bound the range of the each traffic flow function $x_i, \forall i \in N$ and the norm of the difference between the vector of traffic functions in Eq. <ref> as follows:

\begin{align}
\label{eq:proof-lipschitzness-uncongested-network-2}
0 \leq x_i(\theta_t) \leq Q
\implies
-Q &\leq x_i(\theta^1_t)-x_i(\theta^2_t) \leq Q \nonumber \\
0 &\leq \left(x_i(\theta^1_t)-x_i(\theta^2_t)\right)^2 \leq Q^2 \nonumber \\
0 &\leq \|\vx(\theta^1_t)-\vx(\theta^2_t)\|_2^2 \leq N Q^2
\end{align}

where $Q \in \sR_+$ is the total demand, namely, the sum of all cells in the O-D matrix $\mQ \in \sR^{V\times V}$. Now define $|\theta^1_t-\theta^2_t| = \delta >0$ and set $G=NQ^2/\delta > 0$, and replace it into Eq. <ref>:

\begin{align}
0 &\leq \|\vx(\theta^1_t)-\vx(\theta^2_t)\|_2^2 \leq N Q^2 = G \delta  = G|\theta^1_t-\theta^2_t|
\end{align}

which proves that the objective function is $G$-Lipschitz

Now armed with Proposition <ref>, we can provide guarantees for the convergence of when $\theta \in \mathbb{R}$ and under exogenous travel times:

Assume that the objective function of the problem is (coordinate-wise) pseudo-convex
and that the travelers' utility function in the depends on a single exogenous attribute. If the optimization problem is solved with normalized gradient descent (), converges to an $\epsilon$-optimal solution in a finite number of iterations (poly(1/$\epsilon$))

By assumption, the objective function $\ell: \mathbb{R} \to \mathbb{R}$ is pseudo-convex and thus strictly quasi-convex. Under exogenous travel times, we can use Proposition <ref> to conclude that $\ell$ is locally Lipschitz. Thus, all conditions from Theorem 4.2, <cit.> are satisfied and hence the proof is complete.

§ STATISTICAL INFERENCE ON THE PARAMETERS ESTIMATED FROM SYSTEM LEVEL DATA

This section derives closed form expressions of the statistical tests used to analyze the coefficients of the travelers' utility function estimated with the methodology presented in Section <ref>. estimation is often preferred over  because the former relies on weaker distributional assumptions  <cit.>.

§.§ The estimator

The nonlinear regression formulation of the outer level objective of our problem defines each traffic count measurement $\rx_i \in \sR, \forall i \in N$ in the true data generating process as a scalar dependent random variable with conditional mean:

\begin{equation}
\label{eq:response-function-nlls-learning-problem}
\E(\rx_i | \mZ, \bar{\vt}) = x_i(\mZ, \bar{\vt}, \vtheta)
\end{equation}

where $x_i$ is the traffic count (response) function (Eq. <ref>, Section <ref>),  $\bar{\vt}$ is the set of link travel times and $\mZ$ is the matrix of exogenous attributes. Note that $\bar{\vt}$ is assumed exogenous in the problem solved at the outer level but it is iteratively updated via the alternating optimization with the inner level problem (Section <ref>). The regression equation of $\rx_i$ is defined as:

\begin{equation}
\label{eq:regression-equation-nlls-learning-problem}
\rx_i = \E(\rx_i | \mZ, \bar{\vt}) + \ru_i
\end{equation}

where $\ru_i$ is the random error coming from the true data generating process. The differentiation of the objective function of the problem (<ref>) leads to the following first order necessary optimality condition:

\begin{equation}
\label{eq:first-order-optimality-condition-NLLS}
\dfrac{\partial \ell(\vtheta)}{\partial \vtheta}
= \frac{2}{n} \sum_{i \in N} \dfrac{\partial x_i(\vtheta)}{\partial \vtheta} (\rx_i-x_i(\vtheta))
%\dfrac{\partial x(\vtheta)_i}{\partial \vtheta} \ru_i
\vzero
\implies
[D_{\vtheta} \ x(\vtheta)]^\top \ru = \vzero
\end{equation}

where $n$ is
the number of observations and  $D_{\vtheta} \ x(\vtheta)  \in \mathbb{R}^{n \times d}$ is a Jacobian matrix corresponding to the stacked gradient vectors for each observation $n \in N$. Also,  for the ease of notation and for the remainder of the paper $\vx(\vtheta) = \vx(\mZ, \vt, \vtheta)$. Note that Eq. <ref> defines a set of $d$ equations, namely, one equation for each each element of the vector of utility function coefficients $\vtheta \in \mathbb{R}^{|D|}$, and it restricts the residual to be orthogonal to $D_{\vtheta} \ x(\vtheta)$. A key difference of respect to is that the matrix of exogenous attributes of the regression equation is replaced by the Jacobian matrix $D_{\vtheta} \ x(\vtheta)$ (<ref>). Under this realization, the statistical properties of the problem can be derived analogously to the case.

§.§ Assumptions

The proof of consistency of the estimator relies on a series of assumptions that has been well-established in Econometrics literature <cit.>. An accessible reference on these assumptions is found in <cit.>. Below are the assumptions reformulated in the context of our problem:

[Model specification]

The response function is well-specified, i.e. $\forall \vtheta_0 \in \mathbb{R}^{|D|}, \ \rvx = \vx(\vtheta_0, \mZ, \bar{\vt}) + \vu$

[Orthogonality between errors and regressors]

In the data generating process, $\E[\vu|\mZ, \bar{\vt}] = \vzero$ and $\E[\vu \vu^{\top}|\mZ, \bar{\vt}] = \mOmega_0$, where $\mOmega$ is the covariance matrix of the errors terms

Assumption <ref> allows for heterocedastic errors in the problem but it requires knowing or estimating a functional form of the covariance matrix $\mOmega_0$ <cit.>.  To ease the analysis, we introduce an additional assumption:

[Spherical errors]

Errors are homocedastic and non-autocorrelated, that is, their covariance matrix is diagonal $\mOmega_0 = \sigma^2\mI$, for $\sigma \in \sR_{\geq 0}$

Assumption <ref> of spherical errors is satisfied if the errors are independent and identically distributed (iid) and with a constant variance

Each traffic flow function $x_i(\cdot)$ satisfies that: $\forall i \in N, \forall \vtheta_1, \vtheta_2 \in \mathbb{R}^{|D|}$, $x_i(\vtheta_1, \mZ, \bar{\vt})$ = $x_i(\vtheta_2, \mZ, \bar{\vt})$ iff $\vtheta_1 = \vtheta_2$

Assumption <ref> is directly satisfied when $\theta \in \sR$ and the traffic flow functions are (coordinate-wise) monotonic (see Definition <ref>, Section  <ref>).

[Rank of the limit distribution of the Hessian matrix approximation]

The matrix

\mA_0 = \textmd{plim} \frac{1}{N}  \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top  \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]

exists and is finite and nonsingular $\forall \vtheta_0 \in \mathbb{R}^{|D|}$.

Assumption <ref> formalizes the third rule of identifiability discussed in Section <ref>. Proposition <ref> shows that this assumption can be tested via the rank of the Jacobian of the traffic functions:

Suppose the Jacobian matrix of the traffic flow functions given by $D_{\vtheta} \phantom{'}\vx(\vtheta) \in \sR^{n \times d}$, with $n >d$, is full rank at $\vtheta_{\vzero}\in \mathbb{R}^{|D|}$. Then, the square matrix $ \left[D_{\vtheta = \vtheta_{\vzero}} \phantom{'} \vx(\vtheta)\right]^{\top} D_{\vtheta= \vtheta_{\vzero}} \phantom{'}
\vx(\vtheta) \in \sR^{d \times d}$ is positive definite.

Consider $\forall \vu \in \sR^{m} \neq \vzero$:

\begin{align*}
\vu^T \left(2 \left[D_{\vtheta} \phantom{'}\vx(\vtheta)\right]^{\top} D_{\vtheta} \phantom{'}\vx(\hat{\vtheta})\right) \vu
= 2  (D_{\vtheta= \vtheta_{\vzero}} \phantom{'}\vx(\vtheta) \vu)^\top \left(D_{\vtheta= \vtheta_{\vzero}} \phantom{'}\vx(\vtheta) \vu \right)
%	= 2  \langle \vu^{\top}, D_{\vtheta= \vtheta^{\star}} \phantom{'}\vx(\vtheta) \rangle \langle D_{\vtheta= \vtheta^{\star}} \phantom{'}\vx(\vtheta),\vu^{\top} \rangle
= \|D_{\vtheta= \vtheta_{\vzero}} \vx(\vtheta) \vu\|_2^2
\end{align*}

By assumption, the Jacobian matrix $D_{\vtheta = \vtheta_{\vzero}} \ \vx(\vtheta)$ is full rank at $\vtheta_{\vzero}$ and thus $D_{\vtheta = \vtheta_{\vzero}} \ \vx(\vtheta) \vu = \vzero$ iff $\vu = \vzero$. Since $\vu \neq 0$, $\|D_{\vtheta= \vtheta_{\vzero}} \vx(\vtheta) \vu\|_2^2 > 0$ which proves that the matrix $\left[D_{\vtheta = \vtheta^{\star}} \phantom{'} \vx(\vtheta)\right]^{\top} D_{\vtheta= \vtheta_{\vzero}} \phantom{'}
\vx(\vtheta)$ is positive definite and it completes the proof.

Thus, if the Jacobian of the traffic flow functions $D_{\vtheta} \ x(\vtheta)$ has full rank, by Proposition <ref>, the matrix specified in the argument of the probability limit defined in Assumption <ref> is positive definite and hence non-singular.

[Central limit theorem]

$N^{-1/2} \sum_{i = 1}^N [D_{\vtheta} \ x_i(\vtheta)]^\top \ru_i \xrightarrow[]{d} \mathcal{N}[0, \mB_0]$, where

\mB_0
= \textmd{plim} \frac{1}{N} \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top  \mOmega_0 \  \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]
= \textmd{plim} \frac{\sigma^2}{N} \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]

and $\mOmega_0 = \sigma^2 \mI$ by Assumption <ref> of homocedasticity.

When the errors $\vu$ are Gaussian and independently distributed, the linear combination of errors given by $[D_{\vtheta} \ x(\vtheta)]^\top \ru$ is a Multivariate Gaussian. As a result, Assumption <ref> holds even in small samples. Alternatively, if each term $[D_{\vtheta} \ x_i(\vtheta)]^\top \ru_i$ follows an arbitrary distribution, the central limit theorem () must be invoked to ensure the asymptotic normality of the terms in the sequence  $\{D_{\vtheta} \ x_i(\vtheta)\}_{i= 1 \in N}$.  Since $D_{\vtheta} \ x_i(\vtheta)$ is assumed exogenous for statistical inference, the distribution of each term $i$ of the previous sequence is determined by the distribution of $\ru_i$.

§.§ Statistical properties of the NLLS estimator

Suppose $\vtheta_0 \in \mathbb{R}^{|D|}$ satisfies the first order necessary optimality condition defined in Eq. <ref>. Then, under the assumptions stated in Section <ref>, we can establish the consistency of the estimator $\hat{\vtheta}$ at $\vtheta_0$ (see Proposition 5.6, <cit.>) and that

\begin{equation}
\label{eq:convergence-distribution-nlls-estimator}
\sqrt{N}(\hat{\vtheta}-\vtheta_{\vzero}) \xrightarrow[]{d} \mathcal{N}[0, \mA_{\vzero}^{-1}\mB_0\mA_{\vzero}^{-1}]
\end{equation}

\mA_{\vzero} = \textmd{plim} \frac{1}{N}  \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top  \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]
\quad
\mB_0
%= \textmd{plim} \frac{1}{N} \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top  \mOmega_0 \  \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]
= \textmd{plim} \frac{\sigma^2}{N} \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]
= \sigma^2\mA_{\vzero}

By the non-singularity of $\mA_{\vzero}$ (Assumption <ref>, Section <ref>):

\mA_{\vzero}^{-1}\mB_0\mA_{\vzero}^{-1}
= \mA_{\vzero}^{-1}\sigma^2 \mA_{\vzero}\mA_{\vzero}^{-1}
= \sigma^2 \mA_{\vzero}^{-1}

By rearranging terms in Eq. <ref>, the asymptotic distribution of $\hat{\vtheta}$ becomes:

\hat{\vtheta} \xrightarrow[]{a} \mathcal{N}\left[\vtheta_{\vzero}, \Var(\hat{\vtheta})\right]

\begin{align}
\label{eq:asymptotical-distribution-nlls-estimator}
\Var(\hat{\vtheta})
%&= \left(\frac{1}{N}  \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top  \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]\right)^{-1}
%\sigma^2 \left(\frac{1}{N}  \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top  \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]\right)
%\left(\frac{1}{N}  \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top  \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]\right)^{-1}
&= \sigma^2 \left[D_{\vtheta} \ x(\vtheta) \Big|_{\vtheta_{\vzero}}\right]^\top \left[D_{\vtheta} \ x(\vtheta)\Big|_{\vtheta_{\vzero}}\right]
\end{align}

is the asymptotic variance matrix of the estimator. Note that the distribution of the estimator resembles , except that the design matrix $\mX$ is replaced by $\tilde{\mX} = D_{\vtheta = \vtheta_{\vzero}} \ x(\vtheta)$ (<ref>). Also, similar to the case, the statistical behavior of the estimator and the variance of the error terms in large samples can be expressed in terms of $\rvu$ and approximated as <cit.>:

\begin{align}
\
\hat{\vtheta}
%		= (\tilde{\mX}^T \tilde{\mX})^{-1}\tilde{\mX}^T \tilde{\vy}
&= \vtheta + (\tilde{\mX}^T \tilde{\mX})^{-1}\tilde{\mX}^T \rvu \\
%		\to \hat{\rvtheta} \adistr \gN(\vtheta, \Var(\hat{\rvtheta}))
s^2 &= (n-p)^{-1} \rvu^{\top}(\mI-\tilde{\mX}(\tilde{\mX}^{\top}\tilde{\mX})^{-1}\tilde{\mX}^{\top})\rvu
\end{align}

Note that the variance $\sigma^2 \in \mathbb{R}_{\geq 0}$ of the error terms defined in Eq. <ref> is unknown in real applications but it can be consistently estimated as:

\begin{equation}
\label{eq:estimated-variance-nlls}
= \frac{\RSS(\hat{\vtheta})}{N-|K|}
= \frac{\|\bar{\vx} - x(\hat{\vtheta})\|_2^2}{N-|K|}
\end{equation}

where the numerator represents the residual sum of squares () of the model and $|K| \in \mathbb{Z}_{+}$ in the denominator represents the dimension of the $\NLLS$ estimator.

§.§ Hypothesis testing and confidence intervals

Based on the asymptotic normality of the estimator $\hat{\vtheta}$ proved in the previous section, we can derive confidence intervals and hypothesis tests of $\hat{\vtheta}$. The  $1-\alpha$% confidence interval of $\hat{\vtheta}$ can be found by using the percentile $\1-\alpha/2$ of a t-variate with $n-|K|$ degrees of freedom $t_{n-|K|, 1-\alpha/2}$ and the variance of the estimator as shown in the following formula <cit.>:

\begin{equation}
\texttt{CI}(\hat{\vtheta})_{1-\alpha/2} = \hat{\vtheta} \pm t_{n-|K|, 1-\alpha/2}\sqrt{\diag(\Var(\hat{\vtheta}|\mX))}
\end{equation}

The hypothesis $H_0: \theta_d = \theta_{H_0}$ can be contrasted at the $\alpha$% confidence level by computing the statistics $\bar{T}_{d, H_0}$:

\begin{equation}
\bar{T}_{d, H_0}
= \frac{\hat{\theta}_d - \theta_{H_0} }{\sqrt{\diag(\Var(\hat{\vtheta}))_d } }
\end{equation}

and $H_0$ is rejected when $|\bar{T}_{d, H_0}| > |t_{\1-\alpha/2}|$.

§.§ F-test

The F-test is another useful statistic for model comparison in non-linear regression and that allows to contrast hypotheses on multiple model parameters <cit.>. Let's define $\vtheta_{1}, \vtheta_{2} \in \mathbb{R}^{|D|}$ as the vector of coefficients in the restricted (or nested) and unrestricted models. Then, the F-statistic $\bar{F}_{1,2}$ is:

\begin{equation}
\label{eq:f-statistic-model-comparison}
\bar{F}_{1,2} = \frac{\dfrac{\RSS(\vtheta_{1})-\RSS(\vtheta_{2})}{\|\vtheta_{2}\|_0-\|\vtheta_{1}\|_0}}{\dfrac{\RSS(\vtheta_{2})}{N-\|\vtheta_{2}\|_0}}
\end{equation}

where $\|\vtheta_{1}\|_0$ and $\|\vtheta_{2}\|_0$ are the number of non-zero entries in the vectors of parameters of the non-restricted and restricted models. The null hypothesis is that the models are statistically equivalent and it is rejected when $\bar{F}_{1,2} > F_{\|\vtheta_{2}\|_0-\|\vtheta_{12}\|_0, N-\|\vtheta_{2}\|_0,\alpha} \ $, where $F_{\|\vtheta_{2}\|_0-\|\vtheta_{12}\|_0, N-\|\vtheta_{2}\|_0, \alpha}$ is the value of a central $F$-distribution with $\|\vtheta_{2}\|_0-\|\vtheta_{1}\|_0$ numerator degrees of freedom and  $N-\|\vtheta_{2}\|_0^{\phantom{1^2}}$ denominator degrees of freedom and at the $\alpha$-th upper percentile.

An important consideration for hypothesis testing concerns the identifiability of the parameters under the null hypothesis <cit.>. Assume the travelers' utility function is defined as $U = \theta_t t^\alpha$ where $\theta_t,\alpha \in \sR$  are the free parameters and $t$ is the travel time in a given path. If one tests the hypothesis $\theta_t = 0$, $\alpha$ becomes not identifiable. The latter violates Assumption <ref>, Section <ref> on the rank condition and which is required to prove the consistency and asymptotic normality of the estimator. Conversely, under a homogeneous linear-in-parameters utility function (Assumption <ref>, Section <ref>), the identification of the parameters on a restricted model can be guaranteed under mild conditions (Section <ref>).

§ NUMERICAL EXPERIMENTS

To study the performance of the proposed algorithm and some of the mathematical properties of the optimization problem reviewed in Section <ref>, we conduct experiments with synthetic data generated from four small networks. Subsequently, we study the statistical properties of the non-linear least squares estimators () of the utility function coefficients via a series of Monte Carlo experiments conducted in the Sioux Falls, SD network. For the experiments, we employ a validation framework where a hypothetical O-D matrix and a vector of the utility function coefficients are assumed as the ground truth and then used them to generate synthetic traffic count measurements consistent with -. In line with the literature, the utility function coefficients are estimated under the assumption of an exogenous O-D matrix. All the experiments in this section are conducted on a MacBook Pro with Intel Core i5 CPU 2.7 GHz × 2, 1867 MHz 8 GB RAM, 256 GB SSD. The computation time to run all the experiments in this section was approximately 18 hours.

§.§ Data generating process (DGP)

The to conduct the experiments assumes the modeler has perfect information about the link performance functions, the network topology, the paths sets and the ground truth coefficients and attributes of the travelers' utility function. For each network and replicate of the experiment, we generate a set of synthetic traffic counts under an ideal scenario where link flows perfectly matches -. To introduce randomness, we add i.i.d errors to each set of traffic counts measurements. This i.i.d sampling strategy (see Remark <ref>) ensures that Assumptions 	<ref> and <ref>, Section <ref> are satisfied and thus, the statistical inference for the estimator remains consistent. For the sake of convenience, we choose Gaussian errors but other choices of probability distribution are also reasonable, including the Poisson or Negative Binomial distributions. The Gaussian errors are generated with a standard deviation equals to 10% of the average value of the traffic counts. As a reference, the standard deviation of the error term used in the Monte Carlo experiments conducted by <cit.> was equal to the 3% of the average value of the dependent variable reported in his synthetic dataset.

§.§ Small networks

All networks except for the toy network shown in Figure <ref> have been analyzed in prior literature <cit.>. For the toy network, the number of trips between origin-destination pairs is assumed equal to $q_{1,4} = 50, \ q_{2,4} = 100, \ q_{3,4} = 150$ and $q_{i,j} = 0,\  \forall (i,j) \neq \{(1,4), (2,4), (3,4)\}$. For the remaining networks, the O-D matrices are equal to those provided by <cit.>. The travelers' consideration sets in the four networks include all acyclic paths connecting every O-D pair. The links' performance functions are BPR functions with parameters $\alpha = 0.15$ and $\beta = 4$ and the link capacities and free flow travel times are identical to those reported by <cit.>. The utility function is dependent on travel time only, and with a ground truth coefficient $\theta^{\star}_t = -1 < 0$. Note that prior research defines $\theta = -\theta_t >0$ as a positive dispersion coefficient that weights the link/path costs.

[>=Latex[scale=1], scale = 0.8, transform shape]

main node/.style=circle,draw,minimum size=1cm,inner sep=3pt,node distance=3cm

[main node] (C1) [align = center] $1$;
[main node, below = of C1] (C2)  [align = center]  $2$;
[main node, above right = 1.3cm and 3cm of C2] (C3) [align = center]  $3$;
[main node, right = of C3] (C4) [align = center]  $4$;

[draw, ->] (C1) edge node[above]$a_1$  (C3) ;

[draw, ->] (C2) edge node[above]$a_2$  (C3) ;

[draw, ->] (C3) edge[bend left= 60] node[above]$a_3$  (C4) ;
[draw, ->] (C3) edge[bend right= 60] node[above]$a_4$  (C4) ;

Toy network

[>=Latex[scale=1], scale = 0.8, transform shape]

main node/.style=circle,draw,minimum size=1cm,inner sep=3pt,node distance=3cm

[main node] (C1) [align = center] $1$;
[main node, right = of C1] (C2)  [align = center]  $2$;
[main node, below = of C1] (C4) [align = center]  $4$;
[main node, right = of C4] (C3) [align = center]  $3$;

[draw, ->] (C1) edge[bend left = 20] node[above]$a_2$  (C2) ;
[draw, ->] (C2) edge[bend left = 20] node[below]$a_1$  (C1) ;

[draw, ->] (C1) edge[bend right = 20] node[left]$a_{6}$  (C4) ;
[draw, ->] (C4) edge[bend right = 20] node[right]$a_{5}$  (C1) ;

[draw, ->] (C2) edge[bend right = 20] node[left]$a_{7}$  (C3) ;
[draw, ->] (C3) edge[bend right = 20] node[right]$a_{8}$  (C2) ;

[draw, ->] (C3) edge[bend right = 20] node[above]$a_{3}$  (C4) ;
[draw, ->] (C4) edge[bend right = 20] node[above]$a_{4}$  (C3) ;

<cit.> network

[>=Latex[scale=1], scale = 0.7, transform shape]

main node/.style=circle,draw,minimum size=1cm,inner sep=3pt,node distance=3cm

[main node] (C1) [align = center] $1$;
[main node, right = of C1] (C2)  [align = center]  $2$;
[main node, right = of C2] (C3) [align = center]  $3$;
[main node, below = of C1] (C4) [align = center]  $4$;
[main node, right = of C4] (C5) [align = center]  $5$;
[main node, right = of C5] (C6) [align = center]  $6$;

[draw, ->] (C1) edge[bend left = 20] node[above]$a_1$  (C2) ;
[draw, ->] (C2) edge[bend left = 20] node[below]$a_2$  (C1) ;

[draw, ->] (C1) edge[bend left = 20] node[right]$a_{12}$  (C4) ;
[draw, ->] (C4) edge[bend left = 20] node[left]$a_{11}$  (C1) ;

[draw, ->] (C2) edge[bend left = 20] node[above]$a_3$  (C3) ;
[draw, ->] (C3) edge[bend left = 20] node[below]$a_4$  (C2) ;

[draw, ->] (C2) edge[bend right = 20] node[left]$a_{13}$  (C5) ;
[draw, ->] (C5) edge[bend right = 20] node[right]$a_{14}$  (C2) ;

[draw, ->] (C3) edge[bend right = 20] node[left]$a_5$  (C6) ;
[draw, ->] (C6) edge[bend right = 20] node[right]$a_6$  (C3) ;

[draw, ->] (C5) edge[bend left = 20] node[above]$a_{9}$  (C4) ;
[draw, ->] (C4) edge[bend left = 20] node[above]$a_{10}$  (C5) ;

[draw, ->] (C5) edge[bend left = 20] node[above]$a_8$  (C6) ;
[draw, ->] (C6) edge[bend left = 20] node[above]$a_7$  (C5) ;

<cit.> network

[>=Latex[scale=1], scale = 0.6, transform shape]

main node/.style=circle,draw,minimum size=1cm,inner sep=3pt,node distance=3cm

[main node] (C1) [align = center] $1$;
[main node, right = of C1] (C2)  [align = center]  $2$;
[main node, right = of C2] (C3) [align = center]  $3$;
[main node, below = of C1] (C4) [align = center]  $4$;
[main node, right = of C4] (C5) [align = center]  $5$;
[main node, right = of C5] (C6) [align = center]  $6$;
[main node, below = of C4] (C7) [align = center]  $7$;
[main node, right = of C7] (C8) [align = center]  $8$;
[main node, right = of C8] (C9) [align = center]  $9$;

[draw, ->] (C1) edge node[above]$a_1$  (C2) ;

[draw, ->] (C1) edge node[right]$a_2$  (C4) ;

[draw, ->] (C1) edge node[above]$a_3$  (C5) ;

[draw, ->] (C2) edge node[above]$a_4$  (C3) ;

[draw, ->] (C2) edge node[right]$a_5$  (C5) ;

[draw, ->] (C3) edge node[right]$a_6$  (C6) ;

[draw, ->] (C4) edge node[above]$a_7$  (C5) ;

[draw, ->] (C4) edge node[right]$a_8$  (C7) ;

[draw, ->] (C5) edge node[above]$a_9$  (C6) ;

[draw, ->] (C5) edge node[right]$a_{10}$  (C8) ;

[draw, ->] (C5) edge node[above]$a_{11}$  (C9) ;

[draw, ->] (C6) edge node[right]$a_{12}$  (C9) ;

[draw, ->] (C7) edge node[above]$a_{13}$  (C8) ;

[draw, ->] (C8) edge node[above]$a_{14}$  (C9) ;

<cit.> network

Topologies of small networks

§.§.§ Monotonocity of traffic flow functions

Figure <ref> shows the output of the traffic flow functions for four arbitrary links in each network and for $\theta_t \in [-15,15]$. The first two values of the tuples shown in the legend of each subfigure correspond to the origin and destination nodes of a link. The last element of the tuple is an index that distinguish parallel links connecting the same O-D pair, e.g $a_3, a_4$ in Figure <ref>(a). Note that most traffic flow functions are monotonic and some exhibit convex and concave regions that resemble sigmoidal functions. In the toy network, the traffic flow functions of links $(2,3)$ and $(1,3)$ are constant because their link flows are invariant to the value of the travel time coefficient. Therefore, all individuals traveling from node 1 or 2 are forced to traverse those links and thus, the traffic count measurements are not contributing to the identifiability of additional coefficients of the travelers' utility function.

Toy network

Wang network

Lo & Chan network

Yang network

Monotonicity of traffic flow functions in small networks

§.§.§ Pseudo-convexity of the objective function

Figure <ref> illustrates the coordinate-wise pseudo-convexity of the problem. The dashed vertical lines represents the true value of the travel time coefficient $\theta_t^{\star}$. Note that since noise was introduced in the traffic count measurements, the error is not necessarily minimized at $\theta_t^{\star}$. In the four networks, the objective function $\ell$ is monotonically decreasing or increasing for the ranges of values that are lower or higher than $\theta^{\star}$, respectively, meaning that the negative slope of the first derivatives are pointing toward the region where the global optima is located (Figure <ref>). As expected, the changes in curvature of $\ell$ within the feasible domain suggest that $\ell$ is not convex but (coordinate-wise) pseudo-convex.

Objective function

Sign of first derivative of objective function

Pseudo-convexity of objective functions in small networks

The proof of coordinate-wise pseudo-convexity of the problem presented in Section <ref> requires both coordinate-wise monotonocity of the traffic flow functions and a utility function with exogenous attributes only. Our empirical results suggest that the exogeneity of travel time and the monotonocity of traffic flow function are not necessary but only sufficient conditions for the pseudo-convexity of the objective function.

§.§.§ Convergence and consistency in parameters' recovery

To study the advantages of using instead of a second order method in the non-refined stage of the optimization, we interchange the use of the methods between stages. Figure <ref> presents the convergence results obtained by applying and in the no-refined and refined stages, respectively. Figure <ref> shows the converse case. To ensure that the global optima at $\theta_t = -1$ is attainable, no noise is introduced into the traffic count measurements. The learning rate for was set to $\eta =2$ which, for our unidimensional optimization problem on $\theta_t$, equates to solution updates with magnitude equal to 2. The starting point for optimization of the travel time coefficient is set to $\theta_t=-14$, which is the most challenging scenario of convergence analyzed by <cit.>.

No-refined stage: . Refined stage:

No-refined stage: . Refined stage:

Convergence and consistency in parameter recovery with synthetic data from small networks

We observe that the use of in the no-refined stage makes the convergence faster and stable in the four networks. In contrast, the use of in the no-refined stage results into unstable updates of $\theta_t$. This behavior in second order optimization methods is expected due to the constant change of curvature sign and the existence of flat regions within the range of pseudo-convex functions. In fact, for all networks the objective function is concave and flat at the initial value for optimization $\theta_t = -14$ and it is convex and sharp around the global optima at $\theta_t=-1$  (Figure	<ref>).

§.§.§ Hypothesis testing

Table <ref> shows the parameters estimates and t-tests obtained for the four small networks and under four scenarios; (i) only , (ii) only , (iii) in no-refined stage and in refined stage, (iv) in no-refined stage and in refined stage. Overall, this evidence confirms the advantages of the integration of first and second order optimization methods to perform statistical inference. The integration of and in scenarios (iii) and (iv) generate solutions closer to the global optima than the standalone application of each optimization method. The hypothesis tests correctly reject the null at a 1% confidence level but only when and are used in the no-refined and refined stages of the optimization, respectively.

Point estimates and t-tests for travel time coefficient

4c Network

1l Parameter (t-test)
Toy ${\!}$     <cit.>     <cit.>      <cit.>

0.000 (0.0)     $-$2.000$^{***}$ ($-$5.6)      $-$2.000$^{***}$ ($-$6.9)     $-$4.000$^{***}$ ($-$3.6)

$-$14.000 ($-$0.0)     $-$14.000 ($-$0.6)      $-$0.984$^{***}$ ($-$7.3)     $-$14.000 ($-$1.1)

$+$    $-$0.761$^{***}$ ($-$7.1)     $-$2.000$^{***}$ ($-$5.6)      $-$0.976$^{***}$ ($-$15.7)     $-$0.982$^{***}$ ($-$7.3)

$+$    0.000 (0.0)      $-$2.000$^{***}$ ($-$5.6)      $-$0.984$^{***}$ ($-$7.3)      $-$2.000$^{***}$ ($-$6.8)

Links (coverage)     4 (100%)      8 (100%)     14 (100%)      14 (100%)

Paths     6      24      44     28

O-D pairs     3      12      12     9

* Note: Significance levels: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01

§.§.§ Impact of error in reference OD matrix

To study the impact on convergence and inference of assuming an inaccurate reference O-D matrix, the experiment and testing network used by <cit.> were used as a baseline. The ground truth O-D matrix was defined as $\vq^{\textmd{true}} = (q_{1,6}, q_{1,8}, q_{1,9}, q_{2,6}, q_{2,8}, q_{2,9}, q_{4,6}, q_{4,8}, q_{4,9}) = (120,150,100,130,200,90,80,180,110)$ whereas the reference O-D matrix was $\vq^{\textmd{distorted}} = (100,130,120,120,170,140,110,170,105)$. In line with <cit.>, it is assumed that traffic counts are only observed in the following links of the network:  $(3,6), (5, 6), (5, 8), (5, 9), (7, 8)$.

Figure <ref> shows the convergence toward the optimal solution when the reference O-D matrix was assumed equal to (i) the true O-D matrix $\vq^{\textmd{true}}$ (blue curve) or to (ii) the distorted O-D matrix $\vq^{\textmd{distorted}}$ (red curve). Note that in the two scenarios, the algorithm converges close to the global optima ($\theta_t^{\textmd{true O-D}} = -1.058, \ \theta_t^{\textmd{distorted O-D}} = -1.054$) and the objective function ($\ell_\star^{\textmd{true O-D}} = 758.8, \ \ell_\star^{\textmd{distorted O-D}} = 12301.6$) is minimized over iterations. Note that when $\vq = \vq^{\textmd{true}}$, a zero objective function is not attainable because the introduction of random error in the traffic count measurements (see Section <ref>). Despite the bias in the reference O-D matrix, the null was correctly rejected under the two scenarios at a 95% level of confidence (true O-D: $T_{H_0:\theta_t = 0} = -10.3, \ p < 0.01$,  distorted O-D: $T_{H_0:\theta_t = 0} = -4.1, \ p < 0.05)$)

Convergence under true and distorted reference OD matrix

§.§ A medium-scale network: Sioux Falls

The Sioux Falls network is a standard testing bed used by transportation researchers in network modeling studies <cit.>. The network comprises 24 nodes and 76 links (Figure <ref>). We generate 1,584 paths corresponding to the three paths with the shortest distance to travel between each O-D pair. The O-D matrix is obtained from <cit.> and a heatmap with its values is presented in Figure <ref>. The setup to generate the synthetic data is the same as for the small networks except that the exogenous attributes for the monetary cost $c$and the number of intersections $s$ at each link/path are incorporated in the utility function. The attributes $c$ and $s$ are generated as continuous and discrete random variables, respectively, and the utility function coefficients are set to $\theta _t = -1, \ \theta_c = -6, \ \theta_s = -3$. If the units of the travel time and monetary cost attributes are minutes and US dollars, the value of time () become equal to 10 USD$ per hour. By Assumption <ref>, Section <ref>, the set of coefficients $\theta_t, \theta_c, \theta_s \in \sR$ are common among travelers and linearly weighting the attributes $t$, $c$, $s$ in the utility function. The true path set is assumed to be known and thus, the experiments do not perform the column generation phase of the inner level optimization algorithm (Step 1, Algorithm <ref>, <ref>)

Network topology

O-D matrix

Topology and O-D matrix of Sioux Falls network

§.§.§ Pseudo-convexity of the objective function

To have a detailed characterization of the coordinate-wise pseudo-convexity of the outer level optimization problem, Figure <ref> shows plots of the objective function $\ell$(top left), the first of derivative of $\ell$ (top right), the sign of the first derivative of $\ell$ (bottom left) and the sign of the second derivative of $\ell$ (bottom right) respect to a coefficient of the utility function. The vectorized expressions to compute the first and second derivatives of the objective function are presented in <ref> and <ref>. To avoid an excess of overlap of multiple curves when included in the same figure, we only analyze the curves associated to the travel time and cost coefficients. The true values of the coefficients are represented with vertical dashed lines in each figure and the curve associated to each coefficient was generated by fixing the value of the other coefficient to its true value. To satisfy the exogeneity assumption required for coordinate-wise pseudo-convexity in the objective function (Section <ref>), travel times are assumed known and equal to the travel times obtained when generating the synthetic traffic counts at -.

Similar to the results obtained for the small networks (Section <ref>), we observe that the objective function is locally convex and minimized at a point close to $\theta_t = -1, \theta_c = -6$. In addition, the sign of the first derivatives are always pointing toward the global optima, which illustrate this key property of pseudo-convex functions. In contrast to the results obtained in small networks, the sign of the second derivative changes more frequently. The latter may be associated with the higher complexity in larger networks which tends to increase the non-linearity of the traffic flow functions and hence, of the objective function, respect to the utility function coefficients.

Pseudoconvexity of the objective function in Sioux Falls network

§.§.§ Convergence with a multi-attribute utility function

The coordinate-wise pseudo-convexity of the outer level objective function provides theoretical guarantees for the convergence of toward global optima but only when the utility function depends on a single exogenous attribute (Section <ref>). Interestingly, our experimental results suggests the ground truth coefficients of a multi-attribute utility function can be also recovered in the endogenous case. The top and bottom plots in Figure <ref> show the values of the ratio between the travel time and cost coefficients, i.e. the value of time, and the value of the objective function over iterations, respectively. The non-refined and refined stages of the optimization performed 10 iterations of and  , respectively. The starting points for optimization of the utility function coefficients were set to zero. In the exogenous and endogenous cases, we observe that the value of time is perfectly recovered and that the value of the objective function is close to zero. In the exogenous case, travel times are assumed to be known and thus, each iteration reduces to minimize the outer level objective function only. In the endogenous case, the convergence is less stable due to the additional computation of -at each iteration of the bilevel optimization.

Convergence to ground truth value of time in Sioux Falls network

§.§.§ Impact of endogeneity of travel times on statistical inference

We conduct 100 replicates of two Monte Carlo experiments to study the impact of the endogeneity of travel times on the statistical inference of the utility function coefficients. The significance level for all hypotheses tests is set at $\alpha = 0.1$. Monte Carlo experiments are often used by the discrete choice modeling community to empirically study the properties of estimators in travel behavior models <cit.>. We also study the impact on statistical inference of three different setups of the optimization algorithm; (i) is used in the no-refined stage (ii) is used in the no-refined stage, (iii) and are used in the no-refined and refined stages, respectively. In scenarios (i) and (ii), no iterations of the optimization methods are performed in the refined stage.