Ruiwei Chen
email: ruiwei.chen@entpe.fr

Cécile Becarie
email: cecile.becarie@entpe.fr

Ludovic Leclercq
email: ludovic.leclercq@entpe.fr

Learning link marginals from dynamic simulation to calculate sustainable system optimum

This work focuses on dynamic sustainable system optimum (SSO). When carrying out dynamic SSO assignment, instead of minimizing their own travel cost of each user, we aim to define an SSO status to minimize the on-road traffic emissions of the whole system. The method is based on an system-optimization program. This work proposes a simple, intuitive, simulation-based methodological framework for solving the SSO problem. The main contribution of this paper is to propose a statistically reliable learning process of link marginals that permits to derive link performance function from any DTA simulator. Then the obtained functions can be used to compute network equilibrium including SSO. A test case is carried out on a toy network. Numerical results show that with the help of the proposed method, we can reduce 5.79 % of carbon monoxide (CO) emission on the whole system under SSO equilibrium, when compared with UE.

Introduction

Exhaust emissions from road transport are one of the most important sources of air pollution at urban areas. With given traffic demand and infrastructure, various researches are focused on reducing on-road traffic emission via (i) emission pricing, (ii) traffic signal control and (iii) eco-cost-based traffic assignment models, etc. (Chen and Yang, 2012, [START_REF] Ma | Emission modeling and pricing in dynamic traffic networks[END_REF][START_REF] Osorio | Urban transportation emissions mitigation: Coupling high-resolution vehicular emissions and traffic models for traffic signal optimization[END_REF], Aziz and Ukkusuri, 2012[START_REF] Wang | Dynamic traffic assignment: A review of the methodological advances for environmentally sustainable road transportation applications[END_REF]. This paper focuses on the third topic. We propose a methodological framework to determine the dynamic traffic assignment (DTA) aiming at minimizing the total traffic emissions of the whole system. The obtained equilibrium status is called sustainable system optimum (SSO), in reference to the system optimum (SO) equilibrium where travel time is minimized [START_REF] Mo Ghali | A model for the dynamic system optimum traffic assignment problem[END_REF], Peeta and Mahmassani, 1995b[START_REF] Shen | On path marginal cost analysis and its relation to dynamic system-optimal traffic assignment[END_REF][START_REF] Zhen | System-optimal dynamic traffic assignment with and without queue spillback: Its path-based formulation and solution via approximate path marginal cost[END_REF][START_REF] Lu | Dynamic origin-destination demand flow estimation under congested traffic conditions[END_REF]. While calculation methods for SO have received lots of attention in the literature, very few studies focus on the environmental counter part. In recent years, some hEART 2020 Symposium Original abstract submittal pioneering researches propose eco-system optimal DTA models that extend cell-based DTA model (Aziz and Ukkusuri, 2012) and agent-based DTA model [START_REF] Lu | Eco-system optimal time-dependent flow assignment in a congested network[END_REF], with emission considerations. Aziz and Ukkusuri (2012) compute CO emission in each cell during a time interval based on a nonlinear regression model [START_REF] Lin | Impacts of traffic heterogeneity on roadside air pollution concentration[END_REF]Ge, 2006, Zhang et al., 2010), combined with computed average space-mean speed in each cell. In their mathematical program, the SO objective function takes (i) CO emission and (ii) travel time into account.

In classical SO problem, the path marginal cost (PMC) represents the increase of travel cost to the whole system due to an additional vehicle on that path. At SO state, users are on routes with equal and minimum PMC so that no user can shift to any other paths that have less marginal travel cost (Sheffi, 1985). In SSO problem, we consider path marginal emission (PME) and try to minimize it. Based on this concept, [START_REF] Lu | Eco-system optimal time-dependent flow assignment in a congested network[END_REF] propose an integer linear program based on a bottleneck model to compute link travel speed. They extended path marginal delay formulation proposed by [START_REF] Shen | On path marginal cost analysis and its relation to dynamic system-optimal traffic assignment[END_REF][START_REF] Zhen | System-optimal dynamic traffic assignment with and without queue spillback: Its path-based formulation and solution via approximate path marginal cost[END_REF], [START_REF] Lu | Dynamic origin-destination demand flow estimation under congested traffic conditions[END_REF], in order to model path marginal emission. Their numerical results show that when users choose paths under eco-system SO criteria, we can reduce significantly emissions and total travel time of the whole system. However, as mentioned in [START_REF] Lu | Eco-system optimal time-dependent flow assignment in a congested network[END_REF], there are approximation errors when the bottleneck capacities are highly dynamic, especially when the capacies changes due to the traffic signal and/or the queue spills back from the downstream link. Basically, this approach is only functional for network with slow varying traffic conditions and homogeneous vehicle distribution. In dense and large urban network, many signalised intersections and local heterogeneity prevent from using this approach to determine PME.

The main contribution of this work is that we extend the method of Peeta and Mahmassani (1995b) in order to compute time-dependent link marginal emission (LME) and path marginal emission (PME), based on link performance functions learned from DTA simulation results. The core idea is to learn link performance functions (speed v.s. density) from different network loading and then to derive SSO. Briefly, the methodological framework of this work is as follows:

(1) The dynamic SSO problem is formulated as a minimization program with total emission as objective function. COPERT model [START_REF] Ntziachristos | Copert: a european road transport emission inventory model[END_REF][START_REF] Eea | EMEP/EEA air pollutant emission inventory guidebook -2016[END_REF][START_REF] Lejri | Accounting for traffic speed dynamics when calculating copert and phem pollutant emissions at the urban scale[END_REF] is used for computing on-road traffic emission.

(2) We learn the link performance function, i.e., speed v.s. density, from dynamic simulation with multiple demand scenarios. This process permits us to estimate the link response to most loading scenarios and defines the link cost function for equilibrium calculation.

(3) Based on the link-cost function, we can compute time-dependent LME and PME. We can then solve the optimization problem using the classical fixed-point problem formulation and the MSA algorithm (Peeta andMahmassani, 1995a, Sbayti et al., 2007). In fact, for the same Origine-Destination (OD) pair and for each iteration of the MSA loop, users can switch to optimal routes with minimum PME in order to converge towards SSO. The solution approach consists of simulation-based heuristic hEART 2020 Symposium Original abstract submittal and iterative algorithm in order to obtain equal PME of all OD pairs of the whole system.

(4) Finally, a case study is carried out on a toy network during 30 minutes to evaluate the efficiency of the method.

The remainder of this paper is organized as follows. Section 2 overviews the general SSO problem, methodological framework, the calculation of link emission, link marginal emission, path emission and parth marginal emission. Section 3 present the case study and numerical results.

Method

This section presents mathematical formulation, methodological framework and solution algorithm for modeling and solving SSO problem. Table 1 (on the last page) shows the main notations used in this paper.

Mathematical program

The SSO refers to the status where no user of the network can reduce the whole system emission by changing its routes. The path emission is obtained by summing up link emission. The latter can be obtained as the product of (i) link traffic flow and (ii) emission factor computed by COPERT model. Similar to the classical SO problem (Sheffi, 1985, Peeta andMahmassani, 1995b), the mathematical program of SSO problem reads:

min E total = τ o d p E τ odp = τ a E τ a , ∀o, d, p, τ, a (1a) subject to p f τ o,d,p = q τ o,d ∀o, d, p, τ (1b) f τ o,d,p ≥ 0 ∀o, d, p, τ (1c) x τ a = τ o d p (f τ o,d,p × δ τ a odp ) ∀o, d, p, τ, a (1d)
In Equation ( 1a), E total is used to evaluate the total emission of the whole system. E τ odp and E τ a represent path emission and link emission, respectively. Equation (1b) and Equation (1c) are contraints that the traffic flows of the DTA problem should respect. Equation (1d) explains the relation between the link traffic flow and path traffic flow, with the indicator function δ τ a odp which: (i) equals to 1 if link a belongs to the path p departing at time τ from o to d, and (ii) equals to 0 otherwise.

The total duration of the simulation period H is divided to discret assignment periods T that is indexed by τ : τ ∈ [τ 0 , τ 0 + ∆T, τ 0 + 2∆T, . . . , H]. We assume that the demand hEART 2020 Symposium Original abstract submittal is constant over an assignment period for all OD pairs. Furthermore, traffic condition including TT are averaged over the period to transform the DTA problem into a quasistatic one. During T i , the program (1) can be reformulated as Equations ( 2) based on marginal emissions in order to transform the optimization problem as a fixed-point problem [START_REF] Sbayti | Efficient implementation of method of successive averages in simulation-based dynamic traffic assignment models for large-scale network applications[END_REF][START_REF] Lu | Eco-system optimal time-dependent flow assignment in a congested network[END_REF].

π p,τ (P M E p,τ -P M E * p,τ ) = 0; ∀o, d, p ∈ P od,τ , τ ∈ T (2a) P M E p,τ -P M E * p,τ ≥ 0; ∀o, d, p ∈ P od,τ , τ ∈ T (2b) π p,τ ≥ 0; ∀o, d, p ∈ P od,τ , τ ∈ T (2c)
We firstly present in subsection 2.2 the optimization process embedded with MSA method for searching the least PME paths. Then we present in detail the calculation of LME and PME in subsection 2.3 and subsection 2.4.

Equilibrium calculation for one assignment period

The optimization process embeds two layers [START_REF] Sbayti | Efficient implementation of method of successive averages in simulation-based dynamic traffic assignment models for large-scale network applications[END_REF]: (i) the outer layer for updating the candidate least PME paths set of all OD pairs and (ii) inner layer to update the traffic assignment on cadidate paths using classic MSA method (Peeta and Mahmassani, 1995b) based on the calculation of PME. Figure 1 illustrates the simulation-based two-layer optimization process. In this work, we use Symuvia as a trip-based DTA simulator for computing variables that are used to solve the program 1. Symuvia is developed by the LICIT laboratory at Gustave Eiffel University. It is a microscopic simulator based on a Lagrangian discretization of the LWR model [START_REF] Leclercq | The lagrangian coordinates and what it means for first order traffic flow models[END_REF]. The convergence check of the outerloop (Step 4) is based on the total gap on PME during the corresponding time interval. Regarding the innerloop convergence check, we rely on the total PME gap of all paths of the path set. The convergence check of the outerloop (Step 4 in Figure 1) is based on the total gap on PME during the corresponding time interval. Regarding the innerloop convergence check, we rely on the total PME gap of all paths of the path set. We explain in detail the calculations of LME and PME in the following subsection 2.3.

Emission and marginal emission of links and paths

Link level emission

The key step is to calculate E a in Equation (1a). It is calculated as E a = T T D τ a × e, which represents the product of (i) total traveled distance of all vehicles on a and (ii) emission factor e(v). On the one hand, based on the assumption that the traffic demand during [τ 0 , τ 0 + ∆T ] is constant, the traffic flow is obtained as

x τ a = k τ a × v t a (n t a )
. The total traveled distance on link a reads: according to COPERT model, the emission factor can be modeled as a function of average speed: e a (v τ a ). Therefore, the link emission E τ a is calculated as:

T T D τ a = L a x τ a ∆T = L a ∆T k τ a v τ a = ∆T n τ a v τ a .
E τ a = ∆T n τ a v τ a e a (v τ a ). (3)
The link average speed can be obtained via a learning process, in order to get a map from

n to v: v τ z (•) : n τ a → v τ a .
We focus on the spatially-averaged number of vehicls because it represent the link density (k a = na La ), yet it is easy to be obtained via DTA simulator:

n τ = T T T τ σ .
Therefore, at link level during a fixed time interval T , link emission can be hEART 2020 Symposium Original abstract submittal considered as a function of n t a : E t a = F a (n t a ). F t a depends only on (i) time-dependent n on link a and (ii) link characteristics of a (link length, speed limit, jam density, traffic light cycle, etc.). For simplicity, in the following parts of this subsection, we focus on link-level formulations on a during the time interval [τ, τ + ∆T ] and ignore the subscript a and suscript τ . To ensure that the solution of the mathematical program exits, the first and second derivative of F (n) should be non-negative, i.e., ∂F (n) ∂n ≥ 0, ∂ 2 F (n) ∂n 2 ≥ 0. Nonetheless, we cannot ensure the nonnegativity of ∂F (n) ∂n since the term n dv dn is negative as shwon in Equation (4).

∂F (n) ∂n = ∂nv(n)e(v(n)) ∂n = e(v)(v + n dv dn ) + nv ∂e(v) ∂n . ( 4
)
Furthermore, computationally speaking, we cannot solve the problem over the full time horizon H. Therefore, we calculate the optimal solution over T where path flow and traffic condition can be assumed constant (quasi-static approximation). In the latter case, the link marginal definition may drive the system to a solution where link speed are very small as it reduces the travel distance and thus emission. It appears that the term in Equation ( 3) is not the most adequate formulation for quasi-static simulations. We then assume that all vehicles have to travel the full link distance after some time. This assumption seems reasonable since we set the assignment period ∆T to 15 min. Then T T D a is calculated as L a × n a . We can then reformulate the link emission function depending on n, and at the same time, respects the definition of link traffic emission. The link emission is calculated by the product of (i) total traveled distance of all vehicles on link a and (ii) the emission factor, as shown in Equation (5a). The corresponding link marginal emission is then calculated by Equation (5b).

F (n) = L × n × e(v). ( 5a
) LM E = ∂ F (n) ∂n = L × (e(v(n)) + n ∂e(v) ∂n ) = L × (e(v(n)) + n de(v) dv dv(n) dn ) (5b)
There are two keys to ensure the nonnegativity of ∂ 2 F (n) ∂n 2 : (i) to choose a reasonable emission factor model e(v) with d 2 e (v) dv 2 ≤ 0, and (ii) to build link performance function v(n) with non-positive second derivative. On the one hand, based on COPERT model, e(v) is approximated as a polynomial regression of v. In urban areas, the link average speed is usually limited under 60 km h -1 , which ensures the condition of d 2 e (v) dv 2 ≤ 0 [START_REF] Ntziachristos | Copert: a european road transport emission inventory model[END_REF][START_REF] Eea | EMEP/EEA air pollutant emission inventory guidebook -2016[END_REF][START_REF] Lejri | Accounting for traffic speed dynamics when calculating copert and phem pollutant emissions at the urban scale[END_REF]. On the other hand, v(n) can be learned from data obtained from DTA simulators.

Path marginal emission

To compute PME, we take similar assumptions made in (Peeta and Mahmassani, 1995b), stating that the PME is the sum of time-dependent LME. Therefore, based on Equation (5b), the PME can be calculated as hEART 2020 Symposium

Original abstract submittal

P M E τ o,d,p = τ o d p LM E τ a (n τ a )δ τ a odp , ( 6
)
where δ τ a odp = 1 if the link a is on the path p of the OD pair od with departure time τ . Otherwise, δ τ a odp = 0.

Learning of link performance function v-n

The main objective is to define the v -n curve of all links based on data obtained from microscopic DTA simulator. We assume that v(n) represents the intrinsic characteristic of link a. It represents the response of the speed to the level of loading (n) at link level. Here, we assume that we can estimate a single valuated decreasing function for each link of the network from previews simulations (learning process). There are three main concepts for building the v -n curve, learned from DTA simulator results.

(1) We should make sure that the regression curve is monotone on n, and the second derivative should be non-positive:

d 2 v dn 2 ≤ 0.
(2) We need enough observations of (n, v) on links to obtain good quality fits. This requires a thorough and comprehensive design of simulation experiments deriving the learning. To this end, we introduce a full OD matrix with various demand profiles, in order to reach a relatively congested status on the network, yet without gridlock. The simulations are carried out with sufficient iterations to ensure that a maximum number of links have enough sample points (n, v). Furthermore, we classify different values of n into K classes and update the statistical mean values and standard deviations of T T D and T T T of each class. This helps us to get statistically reliable link average speed computed by T T D T T T . This also enables us to eliminate points with extreme values of T T D and T T T .

(3) We rely on quadratic regression to approximate the v -n curve based on all the mean values of n and v in each class. Figure 2

Final expression for the LME

We devide the values of samples of n into K classes on each link. The outputs of the method from Section 2.4 give sample points (n k , v k ) of each class k. We then carry out a quadratic regression between all sample points (n k , v k ) and the final learned function v(n) reads:

v(n) = a × n 2 + b × n + c, ( 7
)
where the values of a and b are bounded in order to respect d 2 v dn 2 ≤ 0. Therefore, at time t, the calculation of LME reads:

LM E t = ∂E(n) ∂n = ∂Lne(v) ∂n = L × (e(v t ) + n t de(v) dv | v=v t dv(n) dn | n=n t ). (8)
With the help of COPERT formula [START_REF] Ntziachristos | Copert: a european road transport emission inventory model[END_REF][START_REF] Eea | EMEP/EEA air pollutant emission inventory guidebook -2016[END_REF][START_REF] Lejri | Accounting for traffic speed dynamics when calculating copert and phem pollutant emissions at the urban scale[END_REF], the generalized formulus for computing the emission factor is considered as a 3 order polynomial, hEART 2020 Symposium Original abstract submittal

e(v) = p 3 × v 3 + p 2 × v 2 + p 1 × v + p 0 , ( 9
)
where p 0 , p 1 , p 2 and p 3 are scalar coefficients, obtained from empirical regression of COP-ERT model [START_REF] Lejri | Accounting for traffic speed dynamics when calculating copert and phem pollutant emissions at the urban scale[END_REF]. By substituting Equation ( 9) and Equation ( 7) into Equation ( 8). We then get Equation ( 10) for computing link level emission in function of n,

E(n) = L × (a 6 × n 6 + a 5 × n 5 + . . . + a 2 × n 2 + a 1 × n + a 0 ). ( 10
)

Simulation results

We carry out DTA simulation under SSO equilibrium condition on a toy network during 30 minutes. The reference case is DTA simulation under UE condition with the same networks and OD matrix. We then compare the total emissions of the whole network under UE and SSO condition.

The test network has 5 links with one OD pair and one traffic light, presented in Figure 3 (left). Figure 3 (right) presents the traffic demand rate (in veh h -1 ). We carry out two DTA simulations under UE and SSO equilibrium condition, respectively. The total simulation horizon H = 30 min. H is divided into 3 successive assignment periods of 10 min. The results are shown in table 2. Results show that when users choose their path based on SSO condition, the total emission of CO is 5.79 % less than the total emission when users are traveling based on UE condition.

Conclusions and perspectives

This work proposes a methodological framework to solve the dynamic sustainble system optimum (SSO) problem. The SSO is formulated as an optimization program in terms of minimizing the total emissions of the whole system. It is then transformed into a fixedpoint problem using classical Beckman transformation and solved with a MSA algorithm.

The key contribution of this paper are twofold. First, we provide new definitions for link marginals that now address the SSO equilibrium instead of the classical SO problem. Second, the calculation of link marginals with dynamic simulation has always been a challenge which is not well-documented in the literature. Here, we propose a new framework to learn the link performance function from multiple simulations on the same network with different loading scenarios. The implementation of the learning process is simple and can be coupled with other networks and/or DTA simulators.

To evaluate the method efficiency we carry out a case study on a simple network during 30 minutes. The proposed method solves the SSO problem and success to obtain a reduction of total emissions from UE status to SSO equilibrium. At SSO status, the total emissions of CO is 5.79 % less than the CO emissions when users choose their routes based on classic user equilibrium (UE). The authors are working on case studies on a real-world network.

Keywords

Figure 2 :

 2 Figure 2: Link v -n curve earned from simulation results. v min and v max are minimum and maximum link average speed. The red points are obtained from sample points (n k , v k ). The numerical number beside each point denotes the total number of points in each class k.

Figure 3 :

 3 Figure 3: Input for case study on a simple network. Left: small network with traffic light, right: traffic demand rate from node A to node B.

:

  Sustainable system optimum (SSO), dynamic traffic assignment (DTA), link marginal emission, link performance function, learning process, COPERT 5 Acknowledgment This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 646592 -MAGnUM project). References HM Abdul Aziz and Satish V Ukkusuri. Integration of environmental objectives in a system optimal dynamic traffic assignment model. Computer-Aided Civil and Infrastructure Engineering, 27(7):494-511, 2012. Linxi Chen and Hai Yang. Managing congestion and emissions in road networks with tolls and rebates. Transportation Research Part B: Methodological, 46(8):933-948, 2012. hEART 2020 Symposium Original abstract submittal

H

  Set of optimal paths (with minimum path marginal emissions) from node o to d departing at τ Total origin node o to destination node d departing at time τ E τ odp Emission due to traffic departing from node o at time τ to node d via path p f τ odp Path traffic flow from origin node o to destination node d departing at time τ , using path p number of vehicle on link a. na = T T Ta ∆T and na = La × ka K The total number of classes of values of n on link a πp,τ Number of users departing at τ using the path p

  On the other hand,

	1. Initialization (OD matrix road network)	Microscopic Simulation with SIMUVIA
		Link flow, density, travel time and LME
		2. Time-dependent least PME path calculation
		3. Update cost function wrt emission
		4. Outerloop convergence reached?	Yes	-Save outputs -Repeat the process time interval -Go to the next
	Yes	outerloop iteration k'+1 No: update path set P
		5. Initialization of Step 1		6. Update path assignment
		No: innerloop iteration k+1	simulation 7. Microscopic
		10. Innerloop convergence reached?		Link flow, density, travel time and LME
		9. Update cost function wrt emission		8. Searching for the least PME path
				MSA method
	Figure 1: Flow chat of optimization process
	hEART 2020 Symposium			Original abstract submittal

Table 2 :

 2 Total computed CO emissions on the network in Figure3. The difference is calculated with respect to emissions obtained when users travel under UE condition.

	Equilibrium type	Total CO emissions (in g)	Relative difference of ∆E = E SSO -E UE E UE
	UE	24672	-
	SSO	23244	-5.79 %

Table 1 :

 1 Basic notations in this work Index for OD pxair from node o to node d via path p with departure time τ Set of paths from node o to d departing at τ

	Notations Explanations
	o, d, p, τ
	t	Current time
	a	Index of a link in a road network
	P od,τ
	P * od,τ


Mahendra Paipuri
email: mahendra.paipuri@univ-eiffel.fr

Yanyan Xu
email: yanyanxu@berkeley.edu

Marta C Gonzalez
email: martag@berkeley.edu

Ludovic Leclercq
email: ludovic.leclercq@univ-eiffel.fr

Marta C González

Calibration of Multi-region MFD Models using Mobile Phone Data

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Since the empirical existence of Macroscopic Fundamental Diagram (MFD) was first reported by [START_REF] Geroliminis | Existence of urban-scale macroscopic fundamental diagrams: Some experimental findings[END_REF] for the city of Yokohama, Japan, several applications like traffic state estimation [START_REF] Knoop | Network transmission model: a dynamic traffic model at network level[END_REF][START_REF] Yildirimoglu | Approximating dynamic equilibrium conditions with macroscopic fundamental diagrams[END_REF][START_REF] Kavianipour | Traffic state estimation in heterogeneous networks with stochastic demand and supply: Mixed lagrangian-eulerian approach[END_REF], perimeter control [START_REF] Keyvan-Ekbatani | Exploiting the fundamental diagram of urban networks for feedback-based gating[END_REF][START_REF] Haddad | Coordinated distributed adaptive perimeter control for large-scale urban road networks[END_REF][START_REF] Ampountolas | Macroscopic modelling and robust control of bi-modal multi-region urban road networks[END_REF][START_REF] Mohajerpoor | H ∞ robust perimeter flow control in urban networks with partial information feedback[END_REF], congestion pricing [START_REF] Gu | Optimal distance-and time-dependent area-based pricing with the Network Fundamental Diagram[END_REF] and cruising for parking [START_REF] Cao | System dynamics of urban traffic based on its parking-related-states[END_REF]Menendez, 2015, Leclercq et al., 2017), etc. were proposed based on MFD approach. The three main elements for any MFD-based models are the underlying MFD, the trip lengths for each macro-path and the path flow distribution. The present work addresses the calibration of all the stated elements from mobile phone data records.

The estimation of MFD for urban networks is far from trivial. There exists two main types of data sources to estimate empirical MFD namely, Loop Detector Data (LDD) and Floating Car Data (FCD). [START_REF] Wang | An empirical analysis of macroscopic fundamental diagrams for sendai road networks[END_REF], [START_REF] Ampountolas | Real-time estimation of critical values of the macroscopic fundamental diagram for maximum network throughput[END_REF] used the LDD to estimate the empirical MFDs for the urban networks of Sendai, China and Chania, Greece, respectively. However, [START_REF] Buisson | Exploring the impact of homogeneity of traffic measurements on the existence of Macroscopic Fundamental Diagrams[END_REF] demonstrated that the slope of MFD depends on the distance of loop-detectors to the downstream traffic signals. More recently, [START_REF] Shim | Empirical evaluation of influential factors on bifurcation in macroscopic fundamental diagrams[END_REF] studied the bifurcations in empirical MFDs by estimating them using roadside detectors hEART 2020 Symposium Original submittal

and [START_REF] Alonso | Effects of traffic control regulation on network macroscopic fundamental diagram: A statistical analysis of real data[END_REF] analyzed the shape of empirical MFDs for the urban corridors using LDD. On the other hand, FCD is more attractive than its counterpart, as it provides vehicular trajectory data. Typically, FCD is provided by either Global Positioning System (GPS) or mobile phones. FCD was used to estimate empirical MFD [START_REF] Geroliminis | Existence of urban-scale macroscopic fundamental diagrams: Some experimental findings[END_REF][START_REF] Bazzani | Towards congestion detection in transportation networks using GPS data[END_REF][START_REF] Tsubota | Macroscopic Fundamental Diagram for Brisbane, Australia: Empirical findings on network partitioning and incident detection[END_REF], traffic monitoring [START_REF] Herrera | Evaluation of traffic data obtained via gps-enabled mobile phones: The mobile century field experiment[END_REF] or travel time estimation [START_REF] Jie | Monitoring travel times in an urban network using video, GPS and bluetooth[END_REF]. [START_REF] Beibei | Determining the Macroscopic Fundamental Diagram on the basis of mixed and incomplete traffic data[END_REF], [START_REF] Ambühl | Data fusion algorithm for Macroscopic Fundamental Diagram estimation[END_REF] estimated an empirical MFD by data fusion of both LCD and FCD. Most of the works that employ FCD use only the GPS data from taxis to estimate MFD, as the GPS data of private cars is not readily available. Hence, the penetration rate of taxis is an important factor to estimate the MFD accurately from FCD [START_REF] Du | Deriving Macroscopic Fundamental Diagrams from probe data: Issues and proposed solutions[END_REF].

Another key ingredient of the MFD-based modeling framework is the set of macro-paths and their corresponding trip lengths. It is not trivial to estimate neither the macro-paths nor their trip lengths using LDD without any further equipment. On the other hand, FCD from taxis can be processed to obtain a distribution of trip lengths. However, FCD is generally sparse and it fails to capture the repetitive trips made by the frequent users/residents. The importance of estimation of accurate and reliable trip lengths in the context of MFD-based simulation was discussed in-detail in [START_REF] Batista | Estimation of regional trip length distributions for the calibration of the aggregated network traffic models[END_REF] and it is as critical as the estimation of the MFD per se for the accurate resolution of traffic dynamics. Similarly, it is difficult to observe the path flow distributions along each macro-path with the existing LDD or FCD data sources. Most often path flow coefficients are estimated using Dynamic Traffic Assignment (DTA) with User Equilibrium (UE), Stochastic User Equilibrium (SUE) or Bounded Rational User Equilibrium (BRUE) equilibration [START_REF] Batista | Estimation of regional trip length distributions for the calibration of the aggregated network traffic models[END_REF]. However, recent studies show that route choice discipline may not valid for macro-paths for large-scale networks [START_REF] Mariotte | Calibration and validation of multi-reservoir mfd models: A case study in lyon[END_REF]. At the same time, there have not been any works proposed to estimate the path flow coefficient from the empirical data. These research gaps can be appropriately addressed by using massive phone data to calibrate the above discussed models.

The present work proposes a method to estimate the empirical MFD of any urban network from the Location Based Service (LBS) data. This type of data is generated by the smartphone Apps, which share their location data actively to the App developer. The major drawback of this data is the diversity of users providing GPS locations using different modes of transportation. Hence, it is important to pre-process the data to extract information on cars and taxis. Another important factor in the case of LBS data is the sampling rate, which can vary from a few seconds to a few minutes depending on the type of App being used. This work proposes a method to enhance the data with poor sampling rates to obtain accurate travel distance. This work also introduces a framework to estimate the time dependent penetration rates. The secondary contribution of this work is the analysis of trip lengths inside the network. A static and dynamic analysis of the trip lengths are elaborated in-detail. The final part of the manuscript discusses the estimation of path flow distribution. The proposed framework is generic and it can be used to process any type of phone data that has records of GPS coordinates and timestamps.

Data Description and Pre-processing

The data contains the positioning of the smartphone devices either by GPS or WiFi Positioning System (WPS) for the city of Dallas, Texas in the United States for a period of 2 months from March, 2017 to April, 2017. In the present work, downtown Dallas and neighboring suburbs as shown in Fig. 1a are considered to calibrate the MFD models. Fig. 1b shows the link level presentation of the area under consideration, which extends across 160 km 2 and contains 18 386 nodes and 48 287 links. The length of the road network is 4 800 km, which includes all types of roads. The raw data contains an anonymized user ID, timestamp, longitude, latitude and uncertainty of the location. The data consists of 3.7 million users and there are around 4.5 billion hEART 2020 Symposium Original submittal

Trip Enrichment Method

The main idea behind the enrichment method is to find the shortest path between the sparse records of a given trip. Two important tools are used in this context namely, OSMnx [START_REF] Boeing | OSMnx: New methods for acquiring, constructing, analyzing, and visualizing complex street networks[END_REF], which is a Python package that is used to analyze the road networks and NetworkX [START_REF] Hagberg | Exploring network structure, dynamics, and function using NetworkX[END_REF], another Python package used to study the dynamics of road networks. The enrichment method is explained with an example in the following. Figure 2a shows three example traces of a sample trip with sparse records from the given data. All three trips have either 3 or 4 records for relatively long trips lengths of around 10 km. The result is a very poor resolution of the trajectory of the trip, which is shown in Fig. 2a. A simple Haversine distance between successive data points of each trip introduces a huge approximation in the traveled distances. At the same time, the traveled time estimation is unaffected due to the sparse records. The combination of these two phenomena can introduce a considerable scatter in the MFD and poor estimation of trip lengths. Hence, it is desirable to map the trajectory of the trip to the underlying road network as closely as possible.

In the current work, the trips are enriched using the spatial geometry of the network. For each trip, the distance between successive records is estimated. If the distance is bigger than a threshold, defined a priori, the location of the nearest nodes close to the GPS positions of those records is obtained. Once the locations of the nodes are obtained between the sparse data records, the shortest path between those two nodes is computed using NetworkX shortest path subroutine. The GPS locations of the nodes in the shortest path are added to the trip between those two sparse records. Fig. 2b presents the traces of considered trips after computing the shortest paths between the sparse records. It can be clearly noticed that the trajectory of each trip is matched to the network after the enrichment process.

Whenever the trace of the trip is enhanced between two sparse records, the timestamps are also interpolated to match the spatial data. The timestamps are interpolated based on the average speed between the two records in the original trace. Hence, this method transforms the sparse spatial and temporal data into dense data, thereby improving the overall accuracy and the representation of the MFD.

Computation of Macroscopic Variables

Partitioning of the Network

Partitioning of the considered network into homogeneous reservoirs is the first step to estimate the macroscopic traffic variables. As the primary objective of the current work is to propose a methodology to estimate the MFD and the trip lengths from mobile phone data, a simple partitioning scheme is assumed. However, it is to be noted that the present framework can be used with any of the partitioning schemes proposed in the literature. The considered network is divided into 5 reservoirs as shown in Fig. 3a. The rationale behind the partitioning scheme is to have one reservoir for the downtown region and divide the region around the downtown into similar sized reservoirs. It is worth noting that the boundaries of reservoirs are not placed along with the road network, but in-between the road networks.

Estimation of Penetration Rates

Since the available mobile phone data accounts only for the part of the total traffic, it is necessary to estimate the penetration rate (ρ) of the vehicles to compute a well-defined MFD and path flow distributions. In this work, penetration rates are estimated by fusing the OD matrix data, LDD and LBS data. The OD matrix data is obtained from North Central Texas Council of Governments for the morning peak period, i.e., 06:30 AM to 08:59 AM. Since this data does not have dynamic information of the OD matrix, LDD is used to transform this static OD matrix into a dynamic OD matrix. Fig. 3c presents the estimated normalized demand profile from the loop-detectors in the network. Now, using the OD matrix data for the given 2.5 h period and the normalized demand profile, it is possible to expand the partial OD matrix to a typical full day OD matrix. hEART 2020 Symposium Original submittal

Table 1: Total number of vehicles (× 10 4 ) for 24 h period between different macroscopic OD pairs.

Reservoir 1 2 3 4 5 1 7.9 3.4 2.6 0.5 0.7 2 4.8 7.2 1.8 0.2 0.4 3 7.0 4.4 6.5 0.4 0.2 4 1.7 0.9 0.7 0.7 0.1 5 4.0 3.0 0.7 0.2 2.2 macroscopic OD pairs, i.e., vehicle counts from one reservoir to another for 24 h period.

According to the generalized definitions of [START_REF] Edie | Discussion of traffic stream measurements and definitions[END_REF], the average network density (k) and flow (q) can be expressed as,

k = N i=1 T T i L n ∆T and q = N i=1 T D i L n ∆T , ( 1
)
where T T i and T D i are time traveled and distance covered on a link l, respectively, L n is the total network length, ∆T is the aggregation interval and N is the total number of links in the network under consideration. In the present work, an aggregation interval of 15 min is used. It is clear from eq. ( 1) that the density is computed using Total Traveled Time (TTT) and the flow is estimated by Total Traveled Distance (TTD) in the network by all vehicles within each aggregation interval. The density and flow of the network can be estimated using eq. ( 1) only if the trajectories of all the vehicles are known a priori. Often, that is not the case in reality and only trajectories of a fraction of the vehicles are available. Hence, it is necessary to correct the expressions in eq. ( 1) with penetration rates. According to [START_REF] Nagle | Accuracy of networkwide traffic states estimated from mobile probe data[END_REF], density and flow of vehicles in the network can be re-defined as,

k = P i=1 T T i /ρ L n ∆T and q = P i=1 T D i /ρ L n ∆T , ( 2
)
where ρ is the penetration rate and it is defined as the number of probe vehicles to the total number of vehicles in the network. In eq. ( 2), the sum is made over the total number of probe vehicle trajectories P . The final step to compute the macroscopic variables is to estimate the penetration rate. The current work proposes two different types of penetration rates namely, OD specific penetration rate (ρ od ) and origin specific penetration rate (ρ o ). It is possible to estimate the trip OD matrix for each day from the mobile phone data based on the departure time of each trip. Since the number of trips between the same OD can vary from day-to-day, trip OD matrix based on the mobile phone data is estimated for each day separately. Let N I od,p be the number of trips from the data between origin o and destination d starting within the aggregation interval I. Similarly, N I od,n is the total number of trips estimated by data fusion of loop detectors and OD matrix from NCTCOG as elaborated earlier within a given interval I. Now, OD specific penetration rate at the aggregation interval I can be defined as,

ρ I od = N I od,p N I od,n . ( 3
)
Similarly, let N I o,p be the total number of trips from the data originating from origin o to all the destinations for a given interval I. In the same way, N I o,n is the total number of trips from hEART 2020 Symposium Original submittal

origin o to all the destinations computed from the NCTCOG data. Using these two quantities, origin specific penetration rate can be expressed as,

ρ I o = N I o,p N I o,n ≡ r d=1 N I od,p r d=1 N I od,n , ( 4
)
where r is the total number of macroscopic reservoirs in the network.

As the OD matrix provides the data of trip production to trip attraction, penetration rates should be used based on the departure time of the trip. Consider a trip i starting within an aggregation interval I between the OD pair od has traveled time and traveled distance of T T i and T D i , respectively. The expanded travel distance (T D e i ) and the expanded travel time (T T e i ) for the whole network considering the two types of penetration rates for that trajectory can be expressed as,

T T e,od i = T T i ρ Ip od and T D e,od i = T D i ρ Ip od , ( 5a
)
T T e,o i = T T i ρ Ip o and T D e,o i = T D i ρ Ip o , ( 5b
)
where suffixes/prefixes od and o represent that the expansion is done by OD and origin specific penetration factors, respectively and I p is the aggregation interval corresponding to the departure time of trip i. Similarly, k and q can be expressed as,

k od = P i=1 T T e,od i L n ∆T and q od = N i=P T D e,od i L n ∆T , ( 6a
)
k o = P i=1 T T e,o i L n ∆T and q o = N i=P T D e,o i L n ∆T . ( 6b
)

Estimated MFDs

The average values of density and flow for each day are estimated using eqs. (6). Only data from the weekdays, which are 43 in total, is considered in the computation of the macroscopic variables to estimate a stable and reproducible MFD. One important factor to consider in this context is the phenomenon of hysteresis in the MFD for the urban networks [START_REF] Leclercq | Macroscopic traffic dynamics under fast-varying demand[END_REF]. It is normal to observe the clockwise hysteresis loops in the MFD due to network heterogeneity, demand pattern, driver's behavior, etc. However, the loading of the network from near empty state, which is observed during late night hours, is more stable and reproducible. Hence, a parabolic curve is fitted for the estimated MFD for each weekday for the data point ranging from midnight till the morning peak hour, which is 08:00 AM in the present case. Only the days that show similar MFD fit characteristics are chosen to compute the mean MFD for each reservoir. In the present case, it is observed that on an average 18 days show similar MFD fits. Fig. 4 shows the estimated flow MFDs for the reservoirs, where the macroscopic variables are computed by the OD specific penetration rate. Firstly, it can be noticed that all the MFD curves are relatively stable with very little scatter. Reservoir 1 experiences the highest flow rate among all the reservoirs in the considered region, which corresponds to the downtown Dallas area. Another important inference to be made from the MFD plots is the presence of clockwise hysteresis loops in the MFDs. Reservoirs 1 and 3 in Figs. 4a and4c, respectively exhibit clearer hysteresis during the morning peak hour compared to others. The hysteresis phenomenon is present in other reservoirs too, albeit the size of the loop is relatively smaller.

Trip Lengths Estimation -Dynamic Analysis

There had not been many works on the dynamic variation of trip lengths using empirical data due to lack of data. The trip lengths within and between the reservoirs tend to be dynamic, where the users tend to prefer the alternative paths during the congestion periods. [START_REF] Batista | Regional dynamic traffic assignment framework with time-dependent trip lengths[END_REF] proposed a framework to estimate the dynamic trip lengths explicitly and concluded that including the dynamic changes in the trip lengths improves the accuracy of the MFDbased simulations. Note that in the present section, all the macro-paths are represented as the reservoir sequence. For instance, a trip that starts in reservoir 1 and ends in 3 by transversing through reservoir 2, it is named as 1 → 2 → 3.

In the present work, the dynamic trip lengths between macroscopic OD pairs are estimated using the starting time of each trip. In order to do so, an aggregate time period of 60 min is considered. For a given OD pair, all the trips that are starting within a given aggregation period are selected and a mean trip length is estimated for that given period. Fig. 5 presents the mean evolution of the trip length for different macroscopic OD pairs. Consider Fig. 5a, where the evolution of internal macro-path 3 → 3 is presented. Two peaks, one at the morning peak hour and another at the evening peak hour, can be clearly noticed in the plot. This trend signifies that the users perhaps take the longer paths during peak hours to avoid the most used routes. Fig. 5b shows the dynamic trip lengths of the macro-paths 1 → 3 and 3 → 1. It is clear from the plots that the dynamic trip lengths evolution of the two macro-paths are nearly symmetric. This is due to the movement of the people from the suburban region to the downtown in the morning and vice-versa in the evening. It is difficult to predict these trends between different OD pairs a priori and appropriately calibrate the MFD models. This type of analysis has not been done before in the literature due to the lack of sufficient and reliable data. Therefore, this framework estimates the mean evolution of different macro-paths, which can be directly used in the MFD simulation framework.

The remaining work includes a deeper investigation of MFDs using two types of penetration rates presented. An analysis on speed MFDs will also be discussed. In addition, static analysis of the trip lengths, their distribution will be presented in-detail. Finally, the path flow coefficients are extracted from the data and corresponding empirical UE gaps are calibrated.

hEART 2020 Symposium

Original submittal

  Link level representation of the Dallas.

Figure 1 :

 1 Figure 1: Dallas network: map of the area and its link level description.

Figure 2 :

 2 Figure 2: Trip enhancement method: Original and enriched traces.

Figure 3 :

 3 Figure 3: Partitioning of Dallas network and normalized demand pattern for a typical day scenario.

)

  to 12:00 PM (e) MFD for reservoir 5.

Figure 4 :

 4 Figure 4: Flow MFD estimates using the OD specific penetration rate for computing the density and the flow.

Figure 5 :

 5 Figure 5: Mean dynamic trip lengths for different macro-paths.

Acknowledgment M. Paipuri and L. Leclercq acknowledge the funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 646592 -MAGnUM project). Y. Xu and M. González acknowledge the support by the Berkeley Deep Drive (BDD) consortium.