# Continuum Dynamic Traffic Assignment Model for Autonomous Vehicles in a Polycentric Urban City with Environmental Consideration

Authors: Bosheng Rong; Hui Zhao; Shaohua Cui; Cuiping Zhang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper proposed a continuum dynamic model for autonomous vehicles in a polycentric urban city by considering the environment impact of traffic emission. The model assumes that homogeneous autonomous vehicles are continuously distributed over the urban areas which tend to choose a path to minimize their total travel cost from origin to destination. To describe the path choice behavior of travelers, we presented the continuum dynamic traffic assignment model which consists of a two-dimensional hyperbolic system of nonlinear conservation laws with source terms and an Eikonal-type equation. The elastic demand is considered using a function which associating each copy of flow with its total instantaneous travel cost. For the environmental impacts, here we consider the influence ofCO emission and include the cost of emission into the actual transportation cost. A solution algorithm for the model is designed as a cell-centered finite volume method for conservation law equations and a fast sweeping method for Eikonal-type equations on unstructured grids. Numerical examples are given to demonstrate the model and the proposed solution algorithm. Further, the results of the travel cost considering CO emissions and not considering CO emissions are compared.


---

## Body

## 1. Introduction

Air pollution has a serious negative effect in many cities, such as Beijing, Tianjin, and Shenzhen, which brings great adverse effects on people’s health. As a major cause of air pollution, the vehicle exhaust emission attracts public attention especially in metropolitans with large population and dense road network. Due to the industrial and traffic emission, people suffer the health related problems and invisibility caused by heavy haze in several metropolitans of China. Therefore, it is necessary to reduce vehicle exhaust emission, especially vehicle exhaust including carbon dioxide (CO2, CO), nitrogen oxides (NOx), volatile organic compounds (VOC), and particulate matter (PM) [1]. Motivated by this consideration, much research has incorporated vehicle exhaust emission with the traffic behavior analysis [2–4].For years, a variety of modeling methods for vehicle exhaust emissions have been proposed, including analysis models, numerical models, and statistical models to evaluate the environmental benefits of various traffic management and traffic control strategies [5, 6]. Most of methods related traffic emissions are calculated on the emission coefficient of each vehicle under actual traffic conditions (e.g., the density, average speed, and acceleration of traffic flows) [7]. This approach is low accuracy, but the speed of calculation emission has been developed. Therefore, these models, such as US federal’s MOBILE model [8], California’s EMFAC model [9], and European’s COPERT III model [10], are useful for estimating large-scale area emissions. However, all of these studies are focused on the traditional vehicles have to be operated by drivers, and few works have considered the situation with autonomous vehicles. To the knowledge of our extent, there is only one paper written by Fagnant and Kockelman [11] which discussed the travel and environmental implications of shared autonomous vehicles. Based on the fact, in this paper, we try to investigate modeling approach for the autonomous vehicles with environmental consideration.For autonomous vehicles, researchers proposed various models to describe the travel choice behaviors in urban area. Chen et al. [12] develop a mathematical model to characterize the routing behavior for autonomous vehicles on a transportation network with heterogeneous traffic stream consisting of both conventional vehicles and autonomous vehicles. Network equilibrium models are proposed for deploying time dependent managed lanes for autonomous vehicles. Chen et al. [13] present a mathematical model for formulate the network equilibrium model to capture the routing behaviors for autonomous vehicles. The model is assumed that the autonomous vehicles follow the user-optimum routing principle in some specific zones in urban networks. In addition, Levin and Boyles [14] developed a multiclass cell transmission model for modeling travel behavior for autonomous vehicles and also developed a car-following model that predicted a triangular fundamental diagram based on reduced reaction times. Motivated by the greater communications available for autonomous vehicles due to the frequency of lane reversals, they present a cell transmission model formulation for dynamic lane reversal with autonomous vehicles [15]. Besides, van den Berg and Verhoef  [16] also discussed the effects of autonomous vehicles using a dynamic equilibrium model of congestion to understand the impact of the increase in capacity and the decrease in the VOT and the heterogeneity of VOTs.Among vehicle exhaust emission, highCO emission has many negative effects on environment, health and economy. It is obviously one of the major concerns for transport planners and decision makers [17, 18]. In addition, according to the data from EPA in 2008, the CO emission of metropolitans has exceeded the standard value set by National Ambient Air Quality Standards (NAAQS) [18, 19]. Environment and vehicle operation characteristics affect CO emissions of vehicles. However, environmental conditions (temperature, humidity, snow, precipitation, etc.) cannot be controlled by autonomous vehicles operators, while operating characteristics (such as vehicle speed and parking times, and route choice behavior of road users) can be controlled [20]. Hallmark et al. [19] found that speed fluctuation and low speed will cause cars to generate more CO exhaust. Further, Stevanovic et al. [21] found that short time delay, minimum stopping times, and medium speed are the most appropriate conditions to reduce emissions and fuel consumption. Therefore, car exhaust CO should be considered during path selection to minimize the total travel cost including travel time and emissions.As previously reviewed, researchers used dynamic traffic assignment (DTA) model to capture the elements of traffic assignment of autonomous vehicles effectively [14]. As DTA models can better study traveler’s behavior and path selection behavior, they play an important role in traffic control and management [22]. As it have been proved that traffic dynamics are closely related to vehicle emissions [23, 24], DTA models can capture more realistic traffic flow characteristics, such as shock waves, expansion waves and queue spillback [25]. Therefore, there should be more attention on the implementation of the DTA model for solving transportation problems by considering environmental factors.According to the different assignment patterns, the DTA model can be divided into dynamic user optimization (DUO) and dynamic system optimization (DSO). In the DUO model, travelers’ route selection behavior is based on the assumption that each user is self-interested and will always move along the shortest time/cost path [26]. There are two distribution patterns in DUO, reactive dynamic user-optimal (RDUO) assignment and predictive dynamic user-optimal (PDUO) assignment [27]. RDUO assignment is that the user selects the path with the lowest cost based on instantaneous travel cost information gained through broadcasting [28, 29]. PDUO assignment assumes that users choose a path to minimize actual travel time/cost by their experience [27, 30]. The dynamic social optimal (D SO) model is to minimize the total travel time/cost of the system [31, 32]. In this paper, we assume that the user selects the path with the lowest cost by obtaining the instantaneous travel time/cost information, for example, from the broadcast.Generally speaking, there are two common traffic modeling methods to solve traffic equilibrium problems: one is the discrete model method and the other is the continuum model method [33]. In discrete modeling approach, where each road link in the network is modeled separately, it is assumed that the demand is concentrated in the hypothetical regional centroid, which is commonly adopted for the detailed planning and analysis of a transportation system [34–36]. Discrete DTA models are used in mathematical programs [32], variational inequality [29, 37], etc. The discrete model can effectively characterize the time-dependent property of traffic demand in transportation system. However, for a traffic network with a large number of links, intersections, confluences, and serious congestion, it is obviously very difficult to model it. Furthermore, some hypothetical questions raised during the initial or conceptual phase of such complex systems study do not require very detailed results.Because of the advantages, the continuum modeling method has recently been used in the modeling of highly dense transportation systems, such as cities with large population density, New York, Los Angeles, etc. With the continuum approach, a dense network is approximated as a continuum, where the users are free to choose their routes in a two-dimensional space [38]. In this modeling approach, the characteristics of a network, such as the flow intensity, demand, and travel cost, can be represented by smooth mathematical functions [39]. Therefore, compared with discrete or microscopic modeling methods, continuum modeling methods require less data to build the model, so the problem scale of large transportation networks is reduced, thereby saving computing time and memory.Actually, there is another class of models we used in the paper, which is elastic demand model. The elastic demand refers to the number of trips between each origin-destination (O-D) pair varies with the network conditions and considers user responses to changes in network costs in an aggregate way, different from rerouting [22]. Cantarella et al. [40] proposed the original differential variational inequality formula for dynamic network user equilibrium with elastic travel demand. Friesz et al. [41] found that elastic demand can be roughly assumed to be a monotonically decreasing function relative to the instantaneous minimum travel cost between each O-D pair. The elastic demand function can better simulate the user’s travel status. A number of discrete DTA models with multiple user classes and elastic demand are found in literature, but few of them are continuum DTA models.In this paper, we develop a continuum DTA model with elastic demand for autonomous vehicles in an urban city where the city has multiple compact central business districts (CBDs). In this model, theCO emission costs are considered in the travel cost. The model is characterized by continuous independent variables, and studies the path selection behavior of urban regional traffic flow on the basis of continuum model. A copy of flow represents the traffic flow movements heading to a CBD or a destination. For convenience, we assume that the traveler’s path selection decision is based on radio broadcast information or route guidance system information and each copy of flow can follow the RDUO principle when travelers try to choose a route with the lowest total travel cost between the departure point and the destination location. With the development of technology, each vehicle could be equipped with devices for detecting exhaust emissions, and recording the emissions at each time spot. Operators for autonomous vehicles can make strategies to control CO emissions. Finite volume method is deployed for the spatial discretization of flow conservation equations and a fast sweeping method for the spatial discretization of Eikonal-type equations on unstructured grids and the explicit total variation diminishing Tunge-Kutta time-stepping method for time discretization. Finally, one numerical example is designed to test the rationality of the model and the effectiveness of the algorithm.The rest of the article is organized as follows. Section2 discusses the formulation of the problem and the mathematical formula of the DTA model. The solution algorithm for the DTA model is presented in Section 3. In Section 4, we show some numerical examples to demonstrate the rationality of the model and the validity of the algorithm. Finally, Section 5 provides a short summary and our future work.


## 2. Problem Formulation

Consider a study urban area of arbitrary shape as shown in Figure1. Within this area, there are a number of CBDs, which are assumed to be sufficiently compact comparing with the whole study area. In the region, the road network is dense enough to be viewed as a 2D continuum. The study region and its boundary are denoted by Ω (inkm2) and Γ0 (km), respectively. To avoid singularity at the facilities, it is assumed that each facility is of finite size and enclosed by a boundary. Let Γcm (in km) be the boundary of the mth (m=1,2,∙∙∙,M) compact CBD and n be the unit normal vector going out of the region.Figure 1
The modeled city.It is assumed that heterogeneous autonomous vehicles in the area are characterized by their choice of trip destination. In the urban area, themth copy of flow represents the autonomous vehicles to themth CBD. We develop a is assumed that heterogeneous autonomous vehicles in the area are characterized by their choice of trip destination. In the urban area, the mth copy of flow represents the autonomous vehicles to themth CBD. We develop DTA model including the cost of vehicle CO emission for a polycentric urban city to investigate the characteristics of urban traffic flow and the path selection behavior of heterogeneous autonomous vehicles in a continuous network in order to achieve the lowest travel costs, in the case where people are environmentally conscious and intentionally want to reduce vehicle exhaust emissions (CO is mainly). We assume that each autonomous vehicle in the urban area knows exactly the current traffic conditions during a trip to take a faster route based on available travel time/cost information (which is instantaneous through broadcasting or predictive through experience). During the entire journey, the environmental pollution caused by the exhaust CO is considered in the travel cost. For each copy of flow, the time-varying demand function (or demand distribution) qm(x,,y,t) (in veh/km2/h) is assumed that it is a function related to the travel cost and represented by(1)qmx,,y,t=qmΦmx,y,t,twhere qmx,y,t (in veh/km2/h) is the demand of each copy of flow at location (x,y)∈Ω who chooses mth CBD at time t and Φmx,y,t (in $) is the total instantaneous travel cost from the each copy of flow demand location (x,y)∈Ω to patronize mth CBD at time t. Based on Equation (1), qmx,y,tdxdy is the traffic demand from the subarea at (x,y)∈Ω at time t to the mth CBD for each copy of flow. It is obvious that traffic demand qmx,y,t is a monotonically decreasing function of travel cost Φmx,y,t for mth copy of flow [42, 43]. As travel costs increase, traffic demand will decrease and vice versa.


### 2.1. Flow Conservation Equations

Our approach to modeling urban traffic flow is based on the assumption that the entire traffic flow can be compared with fluid particle flows that can write appropriate balance or conservation laws. For each copy of flow, the traffic density, flow vector, and traffic demand must satisfy the flow conservation condition inside the domain of the areaΩ. Further, the condition can be presented as(2)ρtmx,y,t+∇∙Fmx,y,t=qmx,y,twhere ρmx,y,t (in veh/km2/hr) is the traffic density of the mth copy of flow at location (x,y)∈Ω at time t, Fmx,y,t=(f1mx,y,t,f2mx,y,t) (in veh/km/h) is the flow vector of mth copy of flow who are heading to mth CBD at time t, and f1mx,y,t and f2mx,y,t are the corresponding flow fluxes in the x and y directions, respectively. The flow vector Fmx,y,t indicates the mth copy of flow movement direction at location (x,y)∈Ω heading to mth CBD at time t in the two-dimensional plane, and the flow intensity for mth copy of flow is determined by (3)Fmx,y,t=f1mx,y,t2+f2mx,y,t2which is the norm of the flow vector at (x,y) at time t. The average speed Umx,y,t (in km/h) can be formulated as(4)Umx,y,t=Ufmx,yexp⁡-βρ2x,y,tIt is a function of the sum of all of the densities at location (x,y)∈Ω and time t, where ρ2x,y,t=∑m=1Mρmx,y,t is the total density of area traffic flow, Ufmx,y (in km/h) is the free flow speed of the vehicles at location (x,y)∈Ω, and β is a scalar parameter related to street features such as road area ratio. The flow intensity Fmx,y,t can also be defined as(5)Fmx,y,t=Umx,y,tρmx,y,t


### 2.2. Path Choice Constraints

We assume that each copy of flow at location(x,y)∈Ω to the destination is based on instantaneous traffic information from a radio broadcasting service or route guidance system to select the path to achieve the minimum traffic cost at time t. The actual transportation cost per unit distance of mth copy of flow travel at location (x,y)∈Ω to destination at time t is denoted by Cm(x,y,t) (in $/km) and for each copy of flow (6)Cmx,y,t=k1Umx,y,t+πmρ~x,y,t+θΨmwhere k (in $) is the value of time and π (in h/km) is a positive constant, ρ~=ρ1x,y,t,ρ2x,y,t,⋯,ρMx,y,t, Ψm is CO emissions rate at location (x,y)∈Ω at time t for mth copy of flow [20], and(7)Ψm=-0.064+0.0056Umx,y,t+0.00026Umx,y,t-502In (6), Cmx,y,t is divided into three part, where k/Umx,y,t represents the cost associated with the travel time, kπm(ρ~(x,y,t) represents other associated costs, such as a preference for avoiding high-density regions, and θΨm represents the cost of CO emissions to the destination. θ (in $) is a conversion factor, converting CO emission rate into emission cost.The travel costΦmx,y,t for the mth copy of flow can be defined as(8)Φmx,y,t=Φmxc,yc,t+minp⁡∫pCmx,y,tdswhere p represents any path from origin (x,y)∈Ω to destination (xc,yc)∈Γcm and Φmxc,yc,t represents the cost of entering the mth CBD from the border of CBD at time t. Similar to Xia et al. [44], we obtain the reactive dynamic user-equilibrium condition as the equilibrated flow pattern in the urban system:(9)Cmx,y,tFmx,y,tFmx,y,t+∇Φmx,y,t=0From (9), the total instantaneous travel cost Φmx,y,t of urban traffic flow from origin to destination satisfies the following Eikonal-type equation:(10)∇Φmx,y,t=Cmx,y,tFrom (9)–(10), we can observe that Fmx,y,t/Fmx,y,t is a unit vector that is parallel to ds along the path and Fmx,y,t∕∕-∇Φmx,y,t, where “∕∕” means that two vectors are parallel.For any used pathp from origin (x,y)∈Ω to destination (xc,yc)∈Γcm for each copy of flow, the total instantaneous cost(11)C-pm=∫pCmx,y,tds=∫p∇Φmx,y,tds=∫p∇Φmx,y,tds=-∫p∇Φmx,y,t∙ds=-Φmxc,yc,t-Φmx,y,t=Φmx,y,t-Φmxc,yc,twhich is independent of the used paths. Therefore the travel costs at time t on all used paths are equal.For any unused routep~ between origin (x,y)∈Ω and destination (xc,yc)∈Γcm for each copy of flow, the total instantaneous cost is(12)C-p~m=∫p~Cmx,y,tds=∫p~∇Φmx,y,tds=∫p~∇Φmx,y,tds≥-∫p~∇Φmx,y,t∙ds=-Φmxc,yc,t-Φmx,y,t=Φmx,y,t-Φmxc,yc,tThe inequality in (12) can be further explained that the fact that for some segments along route p~, the normal vectors Fmx,y,t/Fmx,y,t and ds are not parallel. Therefore, for any unused route the total instantaneous travel cost is greater than or equal to the used routes. Hence this model will guarantee that each copy of flow will chooses a path in the continuum road network in a user-optimal manner by using instantaneous travel information that can represent current traffic flow conditions [45].A dynamic continuum model with elastic demand for a polycentric urban city can now be formulated as follows. Each copy of flow’s flow conservation equation is(13)fx=ρtmx,y,t+∇∙Fmx,y,t=qmx,y,t,∀x,y∈ΩFmx,y,t=-Umx,y,tρmx,y,t∇Φmx,y,tCmx,y,t,∀x,y∈ΩFmx,y,t∙n=0,∀x,y∈ΓΓcmρmx,y,0=ρ0mx,y,∀x,y∈Ωand is coupled with the Eikonal-type equation(14)∇Φmx,y,t=Cmx,y,t,∀x,y∈ΩΦmx,y,t=Φ0mx,y,t,∀x,y∈Γcmwhere Φ0mx,y,t represents the cost of entering the mth CBD from the border of CBD at time t and m=1,2,⋯,M.Although the dynamic continuum model (13) is reactive, travelers decide whether to travel before they set off. Therefore, for travelers who decide to travel, the next decision is the route choice, so traffic demand will not change in the middle.


### 2.3. Hyperbolicity of the Model

Let(15)ρ-=ρ1ρ2ρ3⋮ρM,f-1=ρ1U1v11ρ2U2v12ρ3U3v13⋮ρMUMv1M,f-1=ρ1U1v21ρ2U2v22ρ3U3v23⋮ρMUMv2M,q-=q1q2q3⋮qM,where v-m=(v-1m,v-2m) is a unit vector in the direction of motion for the mth copy of flow. Based on (4), Fm=ρmUmv-m.The conservation laws of (13) can be described as the following system:(16)ρ-t+∇∙F-=q-The linearized form of (16) is(17)ρ-t+Aρ-x+Bρ-y=q-where A and B are, respectively, calculated by(18)A=∂f-1∂ρ-,B=∂f-2∂ρ-From (3), the following formulate is set up:(19)∂Um∂ρ1=∂Um∂ρ2=∂Um∂ρ3=⋯=∂Um∂ρM=dUmdρThe composite Jacobian matrix of the system (16) is defined as Fn=Anx+Bny, where (nx,ny) is a nonzero unit vector. The characteristic polynomial of Fn is(20)J11c1⋯c1c1c2J22⋯c2c2⋮⋮⋱⋮⋮cMcM⋯cMJMMwhere J11=U1ω1+c1-λ, J22=U2ω2+c2-λ, and JMM=UMωM+cM-λ. In these equations, ωm=v1mnx+v2mny and cm=ρmωm(dUm/dρ), 1≤m≤M. The column transformation is conducted for the above characteristic polynomials, where the column i-1 is subtracted by column i(i=M,M-1⋯2). And then we will get the following results:(21)J11J12⋯c1c1c2J22⋯c2c2⋮⋮⋱⋮⋮cM-1cM-1⋯JM-1,M-1JM-1,McMcM⋯cMJMM≕PMλwhere J11=U1ω1+c1-λ, J12=λ-U1ω1, J22=U2ω2-λ, JM-1,M-1=UM-1ωM-1-λ, JM-1,M=λ-UM-1ωM-1, and JMM=UMωM-λ. See Zhang et al. [46]; we know that the characteristic polynomial PM(λ) has M unequal real roots. From above, the system (16) is strictly hyperbolic.


## 2.1. Flow Conservation Equations

Our approach to modeling urban traffic flow is based on the assumption that the entire traffic flow can be compared with fluid particle flows that can write appropriate balance or conservation laws. For each copy of flow, the traffic density, flow vector, and traffic demand must satisfy the flow conservation condition inside the domain of the areaΩ. Further, the condition can be presented as(2)ρtmx,y,t+∇∙Fmx,y,t=qmx,y,twhere ρmx,y,t (in veh/km2/hr) is the traffic density of the mth copy of flow at location (x,y)∈Ω at time t, Fmx,y,t=(f1mx,y,t,f2mx,y,t) (in veh/km/h) is the flow vector of mth copy of flow who are heading to mth CBD at time t, and f1mx,y,t and f2mx,y,t are the corresponding flow fluxes in the x and y directions, respectively. The flow vector Fmx,y,t indicates the mth copy of flow movement direction at location (x,y)∈Ω heading to mth CBD at time t in the two-dimensional plane, and the flow intensity for mth copy of flow is determined by (3)Fmx,y,t=f1mx,y,t2+f2mx,y,t2which is the norm of the flow vector at (x,y) at time t. The average speed Umx,y,t (in km/h) can be formulated as(4)Umx,y,t=Ufmx,yexp⁡-βρ2x,y,tIt is a function of the sum of all of the densities at location (x,y)∈Ω and time t, where ρ2x,y,t=∑m=1Mρmx,y,t is the total density of area traffic flow, Ufmx,y (in km/h) is the free flow speed of the vehicles at location (x,y)∈Ω, and β is a scalar parameter related to street features such as road area ratio. The flow intensity Fmx,y,t can also be defined as(5)Fmx,y,t=Umx,y,tρmx,y,t


## 2.2. Path Choice Constraints

We assume that each copy of flow at location(x,y)∈Ω to the destination is based on instantaneous traffic information from a radio broadcasting service or route guidance system to select the path to achieve the minimum traffic cost at time t. The actual transportation cost per unit distance of mth copy of flow travel at location (x,y)∈Ω to destination at time t is denoted by Cm(x,y,t) (in $/km) and for each copy of flow (6)Cmx,y,t=k1Umx,y,t+πmρ~x,y,t+θΨmwhere k (in $) is the value of time and π (in h/km) is a positive constant, ρ~=ρ1x,y,t,ρ2x,y,t,⋯,ρMx,y,t, Ψm is CO emissions rate at location (x,y)∈Ω at time t for mth copy of flow [20], and(7)Ψm=-0.064+0.0056Umx,y,t+0.00026Umx,y,t-502In (6), Cmx,y,t is divided into three part, where k/Umx,y,t represents the cost associated with the travel time, kπm(ρ~(x,y,t) represents other associated costs, such as a preference for avoiding high-density regions, and θΨm represents the cost of CO emissions to the destination. θ (in $) is a conversion factor, converting CO emission rate into emission cost.The travel costΦmx,y,t for the mth copy of flow can be defined as(8)Φmx,y,t=Φmxc,yc,t+minp⁡∫pCmx,y,tdswhere p represents any path from origin (x,y)∈Ω to destination (xc,yc)∈Γcm and Φmxc,yc,t represents the cost of entering the mth CBD from the border of CBD at time t. Similar to Xia et al. [44], we obtain the reactive dynamic user-equilibrium condition as the equilibrated flow pattern in the urban system:(9)Cmx,y,tFmx,y,tFmx,y,t+∇Φmx,y,t=0From (9), the total instantaneous travel cost Φmx,y,t of urban traffic flow from origin to destination satisfies the following Eikonal-type equation:(10)∇Φmx,y,t=Cmx,y,tFrom (9)–(10), we can observe that Fmx,y,t/Fmx,y,t is a unit vector that is parallel to ds along the path and Fmx,y,t∕∕-∇Φmx,y,t, where “∕∕” means that two vectors are parallel.For any used pathp from origin (x,y)∈Ω to destination (xc,yc)∈Γcm for each copy of flow, the total instantaneous cost(11)C-pm=∫pCmx,y,tds=∫p∇Φmx,y,tds=∫p∇Φmx,y,tds=-∫p∇Φmx,y,t∙ds=-Φmxc,yc,t-Φmx,y,t=Φmx,y,t-Φmxc,yc,twhich is independent of the used paths. Therefore the travel costs at time t on all used paths are equal.For any unused routep~ between origin (x,y)∈Ω and destination (xc,yc)∈Γcm for each copy of flow, the total instantaneous cost is(12)C-p~m=∫p~Cmx,y,tds=∫p~∇Φmx,y,tds=∫p~∇Φmx,y,tds≥-∫p~∇Φmx,y,t∙ds=-Φmxc,yc,t-Φmx,y,t=Φmx,y,t-Φmxc,yc,tThe inequality in (12) can be further explained that the fact that for some segments along route p~, the normal vectors Fmx,y,t/Fmx,y,t and ds are not parallel. Therefore, for any unused route the total instantaneous travel cost is greater than or equal to the used routes. Hence this model will guarantee that each copy of flow will chooses a path in the continuum road network in a user-optimal manner by using instantaneous travel information that can represent current traffic flow conditions [45].A dynamic continuum model with elastic demand for a polycentric urban city can now be formulated as follows. Each copy of flow’s flow conservation equation is(13)fx=ρtmx,y,t+∇∙Fmx,y,t=qmx,y,t,∀x,y∈ΩFmx,y,t=-Umx,y,tρmx,y,t∇Φmx,y,tCmx,y,t,∀x,y∈ΩFmx,y,t∙n=0,∀x,y∈ΓΓcmρmx,y,0=ρ0mx,y,∀x,y∈Ωand is coupled with the Eikonal-type equation(14)∇Φmx,y,t=Cmx,y,t,∀x,y∈ΩΦmx,y,t=Φ0mx,y,t,∀x,y∈Γcmwhere Φ0mx,y,t represents the cost of entering the mth CBD from the border of CBD at time t and m=1,2,⋯,M.Although the dynamic continuum model (13) is reactive, travelers decide whether to travel before they set off. Therefore, for travelers who decide to travel, the next decision is the route choice, so traffic demand will not change in the middle.


## 2.3. Hyperbolicity of the Model

Let(15)ρ-=ρ1ρ2ρ3⋮ρM,f-1=ρ1U1v11ρ2U2v12ρ3U3v13⋮ρMUMv1M,f-1=ρ1U1v21ρ2U2v22ρ3U3v23⋮ρMUMv2M,q-=q1q2q3⋮qM,where v-m=(v-1m,v-2m) is a unit vector in the direction of motion for the mth copy of flow. Based on (4), Fm=ρmUmv-m.The conservation laws of (13) can be described as the following system:(16)ρ-t+∇∙F-=q-The linearized form of (16) is(17)ρ-t+Aρ-x+Bρ-y=q-where A and B are, respectively, calculated by(18)A=∂f-1∂ρ-,B=∂f-2∂ρ-From (3), the following formulate is set up:(19)∂Um∂ρ1=∂Um∂ρ2=∂Um∂ρ3=⋯=∂Um∂ρM=dUmdρThe composite Jacobian matrix of the system (16) is defined as Fn=Anx+Bny, where (nx,ny) is a nonzero unit vector. The characteristic polynomial of Fn is(20)J11c1⋯c1c1c2J22⋯c2c2⋮⋮⋱⋮⋮cMcM⋯cMJMMwhere J11=U1ω1+c1-λ, J22=U2ω2+c2-λ, and JMM=UMωM+cM-λ. In these equations, ωm=v1mnx+v2mny and cm=ρmωm(dUm/dρ), 1≤m≤M. The column transformation is conducted for the above characteristic polynomials, where the column i-1 is subtracted by column i(i=M,M-1⋯2). And then we will get the following results:(21)J11J12⋯c1c1c2J22⋯c2c2⋮⋮⋱⋮⋮cM-1cM-1⋯JM-1,M-1JM-1,McMcM⋯cMJMM≕PMλwhere J11=U1ω1+c1-λ, J12=λ-U1ω1, J22=U2ω2-λ, JM-1,M-1=UM-1ωM-1-λ, JM-1,M=λ-UM-1ωM-1, and JMM=UMωM-λ. See Zhang et al. [46]; we know that the characteristic polynomial PM(λ) has M unequal real roots. From above, the system (16) is strictly hyperbolic.


## 3. Solution Procedure

To solve (13) and (14), we use a cell-centered finite volume method for (13) and a fast sweeping method for (14) on unstructured triangular meshes. We begin with a regular triangulation of Ω: Sh=UΓi∈ShΓi, where Sh denotes by a regular triangulation of the computational domain Ω. Let Nik,lik,Tik denote the kth vertex, edge, and neighboring triangle of Γi, respectively (see Figure 2). We define h≔max⁡Γi∈ShhTi as the size of the discretization, where hTi is the exterior diameter of the triangular cell Ti.Figure 2
A typical unstructured triangular grid.


### 3.1. Finite Volume Spatial Discretization

Here, the total instantaneous travel costΦmx,y,t is to be known at time t and define the cell-centered finite volume method: each cell Γi represents a control volume (CV); the unknowns are stored in the geometric center of the cells.Integrating (13) over Ti and using the Gauss theorem yields(22)∂ρmxi,yi,t∂t+1Ai∮∂TiFmx,y,t∙nids=qΦmxi,yi,t,tIn (22), (xi,yi), Ai, ∂Ti, and ni are the geometric center, area, boundary, and outward unit normal vector of the triangular cell Ti, respectively, and we use (23) to calculate the integral expression ∮∂TiFmx,y,t∙nids:(23)∮∂TiFmx,y,t∙nids=∑k=13F~m∙niklikwhere ∂Γi denotes the three edges of Ti and nik is the outer unit normal pointing to ∂Γi, F~ikm is the numerical flux through the interface lik between two neighboring cells, and lik is the length of the kth edge lik of Ti, k=1,2,3. More details on calculation F~ikm can be seen from Jiang et al. [47].


### 3.2. Calculation of∇Φm on Triangular Elements

Here, we assume traffic densityρmx,y,t is to be known at time t and the cost potential Φ is assumed to be approximated by the following linear function:(24)Φmx,y,tTi=∑k=13N^kΦkmand its gradient is well approximated by the constant vector(25)∇Φmx,y,tTi=∑k=13∇N^kΦkmwhere N^kk=13 are the Lagrangian interpolation bases and Φkmk=13 are the values of Φm at the three vertices of the triangular cell Ti, which are solved through (14). Further details of this procedure are described in Xia et al. [44] and Qian et al. [48].


### 3.3. Time Discretization

In this subsection, we use the explicit total variation diminishing Runge-Kurra time-stepping method for time discretization to discrete (22), and then (13) and (14) can be decoupled. Here, the flux vector Fm in (12) is relevant to the value of ∇Φm. Therefore, we must first obtain the values of (ρm)n and then obtain the values of ∇Φm at all of the grid points of the unstructured triangular meshes at the nth time step. The values of (ρm)n+1 at the (n+1)th can be obtained as follows:(26)ρmin+1=ρmin-∆tAi∑k=13F~m∙niknlik+∆tqminHere, time step ∆t needs to satisfy the Courant–Friedrichs–Lewy condition and m=1,2,⋯,M.


## 3.1. Finite Volume Spatial Discretization

Here, the total instantaneous travel costΦmx,y,t is to be known at time t and define the cell-centered finite volume method: each cell Γi represents a control volume (CV); the unknowns are stored in the geometric center of the cells.Integrating (13) over Ti and using the Gauss theorem yields(22)∂ρmxi,yi,t∂t+1Ai∮∂TiFmx,y,t∙nids=qΦmxi,yi,t,tIn (22), (xi,yi), Ai, ∂Ti, and ni are the geometric center, area, boundary, and outward unit normal vector of the triangular cell Ti, respectively, and we use (23) to calculate the integral expression ∮∂TiFmx,y,t∙nids:(23)∮∂TiFmx,y,t∙nids=∑k=13F~m∙niklikwhere ∂Γi denotes the three edges of Ti and nik is the outer unit normal pointing to ∂Γi, F~ikm is the numerical flux through the interface lik between two neighboring cells, and lik is the length of the kth edge lik of Ti, k=1,2,3. More details on calculation F~ikm can be seen from Jiang et al. [47].


## 3.2. Calculation of∇Φm on Triangular Elements

Here, we assume traffic densityρmx,y,t is to be known at time t and the cost potential Φ is assumed to be approximated by the following linear function:(24)Φmx,y,tTi=∑k=13N^kΦkmand its gradient is well approximated by the constant vector(25)∇Φmx,y,tTi=∑k=13∇N^kΦkmwhere N^kk=13 are the Lagrangian interpolation bases and Φkmk=13 are the values of Φm at the three vertices of the triangular cell Ti, which are solved through (14). Further details of this procedure are described in Xia et al. [44] and Qian et al. [48].


## 3.3. Time Discretization

In this subsection, we use the explicit total variation diminishing Runge-Kurra time-stepping method for time discretization to discrete (22), and then (13) and (14) can be decoupled. Here, the flux vector Fm in (12) is relevant to the value of ∇Φm. Therefore, we must first obtain the values of (ρm)n and then obtain the values of ∇Φm at all of the grid points of the unstructured triangular meshes at the nth time step. The values of (ρm)n+1 at the (n+1)th can be obtained as follows:(26)ρmin+1=ρmin-∆tAi∑k=13F~m∙niknlik+∆tqminHere, time step ∆t needs to satisfy the Courant–Friedrichs–Lewy condition and m=1,2,⋯,M.


## 4. Numerical Experiments

A numerical example based on the urban area shown in Figure1 is designed to illustrate the dynamic continuum model. This area is divided into finite triangular grids following Figure 2. The city has two compact CBDs, where CBD 1 is in (15.38km, 16.52km) and CBD 2 is in (32.81km, 21.42km) and their diameter is 1.5 km. Here, Copy 1 and Copy 2 are travelers to CBD 1 and CBD 2, respectively. The modeling period is 6:00 a.m.–11:00 a.m., i.e., t∈T=[0h,5h]. The elastic traffic demand function qmx,y,t is assumed to be(27)qmx,y,t=qmaxm1-γ1Φmx,y,tgmt,m=1,2In (27), qmaxm is the maximum flow demand, where qmax1=300, qmax2=250, and γ1 is a positive scalar taking 00002 $−1. The function gmt is the variation factor for demand during the modeling period and is defined by(28)g1t=t,t∈0h,2h,2-t2,t∈2h,4h,0,t∈4h,5h,g2t=t,t∈0h,1h,1,t∈1h,3h,4-t,t∈3h,4h,0,t∈4h,5h,The free-flow speed function in (4) is given by (29)Ufmx,y=Umaxm1+γ2dmx,ywhere Umaxm represents the maximum speed of each copy of flow, dmx,y is the distance between a location (x,y)∈Ω and the center of the mth CBD, and Umax1=Umax2=65km/h, γ2=00004km-1. In (4), parameter β=2×10-6. For (6), k=75 $/h, θ=5, and (30)π1ρ~x,y,t=π0ρ2x,y,tρ1x,y,t+ρ2x,y,t2(31)π2ρ~x,y,t=π0ρ1x,y,tρ1x,y,t+ρ2x,y,t2where π0=0.0025h/km. Equations (30) and (31) represent some of the other extra costs other than time costs, such as travelers choosing other paths to avoid congestion.We assume that the initial densityρ0mx,y=0 in the model ∀(x,y)∈Ω, and the cost to enter each CBD Φ0mx,y,t=0,∀x,y∈Γcm,t∈T. The total travelers attracted to two CBDs at time t and the corresponding cumulative demand at time t are calculated as follows:(32)qΩmt=∬ΩqmΦmx,y,t,tdxdyQΩmt=∫0tfCBDmξdξwhile the total cumulative inflow to the CBD through Γcm at time t is expressed by(33)FCBDmt=∮Γcm∫0tFm∙nx,y,ξdξdsTo verify the convergence of the numerical method for the proposed model, we test three different meshes of triangular elements (Mesh1: 3,501 nodes and 6,759 elements; Mesh 2: 5,281 nodes and 10,267 elements; Mesh 3: 12,168 nodes and 23,916 elements). Figure 3 shows the curves of cumulative inflow for each copy of flow in 3 cases (Mesh 1–3) and Figure 4 shows the cumulative demand. In Figure 3, the relatively stable increase of the curve indicates that more and more travelers are entering the CBD. At t=4, the curve tends to be flat, which means that travelers have basically entered the CBD, and there are almost no traffic flows in the area. In Figure 4, similar to Figure 3, the curve tends to be flat at about t=3.5. It can be seen from Figures 3 and 4 that the curves for Meshes 1–3 in each subfigure nearly coincide, which indicates good convergence among the different finite element schemes of the numerical method for the proposed model.Figure 3
The plot of the cumulative inflowFCBDm(t).
(a) (b)Figure 4
The plot of the cumulative demandQCBDm(t).
(a) (b)The following numerical examples are based on an unstructured triangular grid with 3,993 nodes and 7,743 elements. Figure5 shows the process of the 1th copy of flow density distribution and dissipation at four time phases. In Figure 5, it can be seen that in the first phase (t=1h), the density of Copy 1 is extremely low and does not appear near CBD 2, where it is almost a free flow state. Then when t=2h, as more travelers arrive at their destination, more serious traffic congestion has occurred around CBD 1 and more travelers wait in line for entry to the CBD 1. Copy 1’s density is higher in the eastern part of CBD 1’s boundary, because more travelers from the east side of the CBD to the CBD and some are waiting in line to enter the CBD. In order to reduce the total travel cost, travelers try to choose the route close to the CBD 1 despite the risk of traffic congestion. In the third phase (t=3h), we can find that traffic changes around CBD 1 are more congested, and the traffic density in the entire urban area has begun to decrease relative to phase two. Most travelers to CBD 1 have reached the CBD 1 boundary and try to enter CBD 1, and the rest of the travelers are on the distributed area of the road network. At last phase, the density of CBD 1 in the whole area is very low, because travelers have entered CBD 1. At this time, the traffic demand of CBD 1 is extreme low and Copy 1 recovers to the free-flow state. After t=4.5h or so, there is no traffic on the road network.Figure 5
The different time density of CBD1.
(a)
t=1 h (b)
t=2 h (c)
t=3 h (d)
t=4 hFigure6 shows the density changes of Copy 2 at different times in two cases. Shown on the left is the density changes of Copy 2 after taking into the cost of CO emissions at four phases, where on the right is density changes of Copy 2 in the four phases only travel time costs calculated. At the beginning, in two cases, there is a slight traffic congestion on the west side of CBD 2. Compared with case 2, there are more travelers from the border of CBD 1 to CBD 2 in case 1. In the first phase, there is no high traffic density in the whole area. In the second phase (t=2h), in two pictures, serious traffic congestion has occurred around CBD 2, and case 1 are more serious on the western side of CBD 2. The travelers are waiting to enter CBD 2 at this time. However, compared with case 1, there are more travelers in case 2 from the upper boundary of the region to CBD 2, and traffic congestion occurs more evenly around CBD 2. In third phase (t=3h), traffic congestion around CBD 2 is further aggravated in this two cases. Similar to second phase, compared to case 1, the density of travelers in the whole area in case 2 is higher, and traffic congestion around CBD 2 is almost uniform and the scope of congestion is higher in case 2. Owing to the CO emission cost, the travel cost of each traveler has increased. Therefore, some travelers cancel the travel decision. Further, with the increase of unit cost, travelers tend to travel the shortest route to reach the destination, so the density of Copy 2 relatively more gather on the region. At last phase, in case 1, travelers almost all have entered CBD 2, and traffic density in the area is decreasing until all travelers reaching CBD 2. On the contrary, in case 2, though most of the regions have lower demand, there is still a slight traffic congestion around CBD 2, because some travelers are not entering the CBD 2 entirely. In case 2, the density range of Copy 2 is larger and the congestion in the previous phase is more severe, so travelers will take longer time to enter CBD 2.Figure 6
Density contrast in two cases.
(a)
t=1 h (b)
t=1 h (c)
t=2 h (d)
t=2 h (e)
t=3 h (f)
t=3 h  (g)
t=4 h (h)
t=4 hThe flow vectorF plotted in Figure 7 clearly illustrates the movement of travelers based on reactive conditions of the continuum road network. Without congestion effect, the travelers will travel to the CBDs along a straight-line trajectory, which can be seen at first and last phase. Because there is a low flow density and flow congestion in the entire area, many travelers’ trajectories are close to a straight line. However, in the second and third phases, the travel trajectory of some travelers around the CBD has completely deviated from the straight line because travelers want to avoid congestion by using less congested areas of the city. In the third phase, clear direction of traveler’s travel path may not be observed near one of the CBD nodes for both cases because of serious congestion.Figure 7
Flow vector plot.
(a)
t=1 h (b)
t=2 h (c)
t=2 h (d)
t=2 h (e)
t=3 h (f)
t=3 h (g)
t=4 h (h)
t=4 hFigure8 shows the contours of total travel cost for Copy 1 heading to CBD 1 at different times (t=1h, t=3h), and Figure 9 shows the contours of total travel cost for Copy 2 heading to CBD 2 at different times (t=2h, t=4h). In Figure 8, at the beginning of the period, with the low density of Copy 1 in the area, the each copy flow is in the free-flow state. Therefore, the total cost contours around CBD 1 is almost equal, and it becomes a series of circular type. When t=3h, with the increase of traffic demand, the whole region becomes more congested and then travel delays will increase substantially. The serious congestion show in the eastern part of CBD 1, so the cost near the east of CBD 1 has increased. At the eastern edge area, the maximum cost has increased.Figure 8
The contour plot of the travel cost for Copy1 traveling to CBD 1.
(a)
t=1 h (b)
t=3 hFigure 9
The contour plot of the travel cost for Copy2 traveling to CBD 2.
(a)
t=2 h (b)
t=4 hIn Figure9, the distribution of cost contour around CBD 2 is not symmetrical at t=2h, which indicates that traffic congestion has occurred around CBD 2. Further, the cost contour around CBD 1 is highly distorted and much denser than the other parts of the study area because the area around CBD 1 is highly congested. At last phase (t=4h), the demand for traffic in the region is decreasing, and the each copy flow becomes the free-flow state, so there is a symmetrical distribution of the total cost contours around CBD 2 and the maximum cost in the west side of the region is also reduced.Figure10 shows the total travel cost changes of Copy 1 in different times in two cases. Shown on the left is the total travel cost changes of Copy 1 in case 1 at four phases. Based on four phases, it is obvious that the total cost in case 1 is higher than that in case 2 and, the farther away it is from CBD 1, the more obvious this phenomenon is. In case 1, the total cost from the easternmost side of the area to CBD 1 is even more than twice that of case 2, such as the last phase (t=4h), because cars will emit more CO emissions with the increase of the travel distance, which will further increase the total travel cost. It can be seen that in crowded cities, the cost of exhaust emissions greatly increases the travel cost of users, which will make users more cautious in route selection, and even some users who travel long distance will choose other ways to travel, such as public transportation, which to some extent reduces traffic congestion and CO emissions.Figure 10
The contour plot of the travel cost in two cases.
(a)
t=1 h (b)
t=1 h (c)
t=2 h (d)
t=2 h (e)
t=3 h (f)
t=3 h (g)
t=4 h (h)
t=4 h


## 5. Conclusions

We have considered a city with several CBDs that are competing for homogeneous autonomous vehicles. These vehicles are distributed continuously over the whole area. Within the city region, the road network is relatively dense and even continuum. The DTA with elastic demand for a polycentric urban area analyzes the characteristics of urban traffic flow in continuum time-varying transportation systems. During the entire trip, the emission costs ofCO is considered into the total travel costs to let travelers take environmental factors into account when making path choices. The model satisfies the reactive dynamic user equilibrium principle and it describes the route selection behavior of each copy of flow. In this paper, the total cost includes the travel time and emission cost. The model is solved by a cell-centered finite volume method for conservation law equations and a fast sweeping method for Eikonal-type equations on unstructured meshes. Numerical results verify the applicability of the proposed model to investigate the characteristics of multiple copies of urban traffic flow and the effectiveness of the solution algorithm. Finally, we use the characteristics of density, fluxes, and so on to analyze the characteristics of the traffic flow of vehicles in the city and compare the difference of the traffic flow density after the increase of the CO emission cost.Compared to the discrete model approach, the continuum model approach requires less data throughout the model set-up process, and the discrete model requires detailed data for all road links, including intersections. Therefore the continuum model applies to urban areas with high population density and road density, like some metropolitans in Asian countries. In these regions, vehicles with very high density have a large amount of exhaust emissions which influence the environment heavily. Thus,CO emissions are included in the cost. In this paper, travelers are divided into two types in a highly intensive transportation system, and each type of traveler makes its own path selection in a reactive manner on the basis of desire to minimize their total instantaneous travel cost from origin to destination. As instantaneous decisions concerning route choice based on continuous information can lead to equilibrium paradoxes, in future work, we will consider the paradoxes existing in continuum models to capture the elements.



---

*Source: 8345979-2018-11-13.xml*