Bounding the efficiency gain of differentiable road pricing for EVs and GVs to manage congestion and emissions

Increasing concerns about air pollution and the promise of enhancing energy security have stimulated the growth of electric vehicles (EVs) worldwide. Compared with gasoline vehicles (GVs), EVs have no emissions and are more environmentally friendly to the sustainable transportation system. Since these two types of vehicles with different emission externalities and observable differences, in this paper, we propose a differentiable road pricing for EVs and GVs to simultaneously manage congestion and emissions by establishing a two-class bi-objective optimization (TCBO) model. First, we investigate whether the differentiable road pricing can induce user equilibrium pattern into a unique pareto-efficient pattern. Then performance of the bi-criteria system optimal is measured by bounding the deviation gap of the Pareto frontier. Specifically, we bound how far the total system travel time and total system emissions at a given Pareto optimum can deviate from their respective single-criterion based system optimum. Finally, we investigate the maximum efficiency gain of the bi-criteria system achieved through implementing differentiable road pricing by comparing total system travel time and total system emissions under two states. After defining two types of price of anarchy (POA), the theoretical bound for the worst possible ratio of total system travel time\ total system emissions in user equilibrium state to the total system travel time\ total system emissions in Pareto-efficient state is derived out. In order to validate the feasibility of theoretical bound, we conduct case studies to calculate the numerical bound of POA based on two Chinese cities: Shenzhen and Lasa. Overall, quantifying the maximum efficiency of differentiable road pricing is beneficial for improving the network designing, policy implementation and social efficiency with regard to congestion and emissions caused by EV and GV users.


Introduction
Selfish travelers usually pursue their own optimal strategies and will not achieve the social optimum, which has led to great losses for the whole society. It is shown that the system optimum can be decentralized as a user equilibrium by charging the optimal toll on each link, which is equal to the marginal social cost and marginal private cost [1][2]3]. Since Pigou [4] first introduced marginal-cost pricing, it has been one of the most efficient ways to manage traffic demand and has been a long-standing research topic. The existing literature has investigated marginal cost pricing in multiclass traffic assignment. Netter [5] and Dafermos and Sparrow [6] discussed marginal cost pricing for different vehicle types, such as trucks, passenger cars and buses, in terms of vehicle size. Moreover, enhanced awareness of environmental protection has led to the wide adoption of ecofriendly electric vehicles all over the world. Global EV sales exceeded 5.1 million in 2018, which was almost double the number of new electric car sales in 2017 [7]. The International Energy Agency (IEA) forecasts that the number of electric vehicles on the road around the world will hit 125 million by 2030 [8]. Since EVs and GVs can be physically distinguished and have different impacts on total system emissions, we should consider differentiable marginal social pricing for EVs and GVs.
Furthermore, the multiclass problem has been extended into multiclass multicriteria problems by integrating travel time cost and monetary cost simultaneously. For example, Yang and Huang [9] and Guo and Yang [10] proved the existence of valid tolls that can decentralize multiclass multicriteria network equilibrium flow patterns when travelers value time differently, and the system cost can be measured in terms of either time or monetary units.
Environmental issues and sustainable transportation have attracted ever-increasing attention worldwide in recent years. Congestion and emissions are two types of disutility that negatively impact our daily life in various ways, from being a minor health inconvenience to affecting job satisfaction. An increasing number of studies consider both congestion and emissions as two aspects of traffic assignment problems. Yin [11] that emission-oriented road pricing for single class users can be used to induce the flow pattern with minimum emissions and established a bi-objective optimization model to identify the Pareto optimal solutions for total travel time and total emissions. Wen [12] proposed a bilevel pricing model to minimize emissions and the total travel time based on the dynamic user equilibrium. Wang [13] reviewed the development of dynamic traffic assignment models considering emission road pricing. Chen and Yang [14] established the necessary and sufficient conditions for road pricing to decentralize a Pareto-efficient flow pattern in terms of congestion and emissions.
The methods used to bound the efficiency gain of a pricing scheme have been proven to be very important and meaningful for evaluating pricing schemes. The upper bound of the ratio of the social welfare between the competitive equilibrium and the social optimum is established to quantify the inefficiency of the oligopolistic equilibria via higher tolls and less demand [15]. The efficiency loss of a general second-best pricing scheme is obtained by bounding a given road pricing scheme [16]. The upper bound of the User equilibrium-Cournot-Nash (UE-CN) mixed equilibrium traffic assignment is derived by analytic derivation when road pricing is considered (not considered) part of the total travel costs [17].
The efficiency gain of differentiable road pricing can be bounded by considering the notion of the price of anarchy (POA), which is a concept used in economics and computer science to qualify the inefficiency of the Nash equilibrium by seeking the worst possible ratio between the total cost incurred by untolled travelers and the system optimum in terms of congestion and emissions under differentiable tolls. After Koutsoupias and Papadimitriou [18] proposed the concept of the price of anarchy, Roughgarden [19] proved that the price of anarchy is independent of the network topology when the cost function is linear and nonlinear.
Research on qualifying the POA under different types of cost functions can be classified into two types. The first type is related to the theoretical bound of the POA and includes a series of studies, such as [18][19][20][21][22][23][24]. Correa et al. [25] proved the bound of POA through a geometric method, and the results are consistent with those in [24]. By introducing valid tolls, Karakostas and Kolliopoulos [26] indicated that POA decreases when link time functions are continuous and strictly increasing. Investigating a multiclass multicriteria problem, Han and Yang [27] derived the bound of the POA in different cases: the tolls were both considered and not considered as a system cost under the time criterion and cost criterion functions, respectively. The second type is related to the numerical bound of the POA, which can be calculated by considering large-scale networks. Youn et al. [28] estimated the POA by analyzing the road networks of Boston, London and New York, and the simulation results showed that this value can reach 1.3. Grange et al. [29] provided the results of case studies on the POA in three networks, considering the differences between the demand levels of each network. Liao et al. [30] carried out a simulation and calculated the value of the POA based on the Sioux Falls network. O'Hare et al. [31] proposed four types of mechanisms that explain the relationship between the POA and travel demand in urban networks when the link delay functions are increasing and separable. Zhang et al. [32] carried out a case study and calculated the value of the POA using actual road network traffic data from Eastern Massachusetts.
In summary, road pricing has been considered in the context of a multicriteria (time vs cost) equilibrium, different vehicle types in terms of size, and a bicriteria (congestion and emission) system optimum. However, to the best of our knowledge, no attempt has been made to look into differentiable road pricing for EVs and GVs to simultaneously manage congestion and emissions. Moreover, the existing literature has never qualified the efficiency gain of such differentiable road pricing for EVs and GVs in terms of congestion and emissions. To fill these gaps, this paper will propose differentiable road pricing for EVs and GVs by establishing a two-class bi-objective optimization (TCBO) model to simultaneously manage congestion and emissions. Since EVs and GVs can be physically distinguished, differentiable road pricing for EVs and GVs is feasible. After investigating the decentralization of differentiable road pricing, the theoretical upper bound of POA T will be derived to qualify the efficiency gain of differentiable road pricing in terms of system total travel time. Finally, we report the numerical bounds of POA T and POA E by conducting simulations based on the actual networks in Shenzhen.
Here we use Fig 1 to explain the organization of this paper: Section 2 will provide the notations and preliminaries of this research. Section 3 will investigate whether the differentiable road pricing for EVs and GVs can decentralize a pareto-efficient flow pattern as user equilibrium. To be more specific, after implementing differentiable road pricing, travelers still select the route based on Wardrop principle, but pareto-efficient pattern can be achieved from the perspective of system manager. In Fig 1, any user equilibrium flow pattern such as UE 1 and UE 2 can be induced into a pareto-efficient pattern (point C). Section 3 will derive the deviation gap of a point on the pareto frontier and explore the following questions: when the point C reach minimum total system travel time (Point A), how far could the corresponding total system emissions deviate from its minimum value? Conversely, when total system emissions are minimized (Point B), how far could the corresponding total system travel time deviate from its minimum value? Section 5 will qualify the efficiency gain of differentiable road pricing by comparing the total system travel time/ total system emissions between the state without tolls and the state with differentiable tolls based on the concept of price of anarchy (POA). In Fig 1, since differentiable road pricing can induce point UE 1 and point UE 2 to point C located on pareto-efficient frontier, we will bound the ratio of total system travel time/ total system emissions at point UE 1 to total system travel time/ total system emissions at point C by defining POA. After obtaining the theoretical bound of POA, the case studies based on two cities will be conducted to calculate the numerical bound of POA and validate the effectiveness of general bound. Section 6 will conclude with a summary of the main findings and suggestions for future studies. Table 1 summarizes the parameters, variables and abbreviations used in this paper to describe the mathematical models related to the differentiable road pricing scheme.

Preliminaries and definitions
In this paper, a mixed urban network G = (N, A) consisting of pure electric vehicles (EVs) without emissions and gasoline vehicles (GVs) is considered, where N and A represent the set of nodes and links, respectively. We assume that the travel demand d w for each OD pair w 2 W is fixed, and the percentage of EVs is P. The set of feasible path flow patterns for the EVs and GVs are expressed in Eqs (1)- (4): O where d w;E ar =d w;G ar indicates (0 or1) whether route r (r 2 R w ) chosen by EV/GV users passes through link a (a 2 A). Moreover, the relationship between EV and GV users' path travel cost and link travel cost can be expressed as shown in Eqs (5) and (6): When there is no road pricing scheme, a traveler's link travel cost is c E a ðvÞ ¼ kt a ðvÞ. Under differentiable tolls τ ¼ ðt G a ; t E a ; a 2 AÞ, an EV user's cost is c E a ðvÞ ¼ kt a ðvÞ þ t E a , and a GV user's cost is c G a ðvÞ ¼ kt a ðvÞ þ t G a . This research investigates whether path flow vector f � ¼ ðf E� rw ; f G� rw Þ and minimum OD cost vector μ ¼ ðm E w ; m G w Þ exists, so that for each r 2 R w , w 2 W, the mixed equilibrium conditions for EVs and GVs hold (Eqs (7) and (8)).
where m E w and m G w represent the EV and GV users' minimum path travel costs for OD pair w, w 2 W. The BPR (US Bureau of Public Roads) link time function is used in this paper and assumes that the value of time for EV and GV users (κ) is the same.
Assuming that the value of time for EV and GV users is the same, this paper uses the BPR (US Bureau of Public Roads) link time function (9): Since more than 90% [32] of the carbon monoxide (CO) in urban centers is from vehicular emissions, Wallace et al. [33] used a nonlinear macroscopic model to estimate the emission load of link-based vehicular CO, and the coefficients are equivalent to those in TRANSYT-7F. The emission functions of such a form (10) are widely used; for example, they are used in [11,[34][35][36].
As long as t a (v a ) is strictly increasing and the speed L a /(t a (v a )) < 1/0.7962 = 75 km/h, In other words, the total CO emission function e a (v a ). v a is strictly convex. The relationship between e a (v a ) and t a (v a ) is given in Eq (11): where ROP a � l a = 0.2038 � exp(0.7962(l a /t a (v a ))), l a denotes the length (km) of a link and t a (v a ) denotes the link travel time. If the maximum speed is 75 (km/h), then η(v a ) � 0.5513, and thus, η min can be obtained. If the minimum speed is set to 0 (km/h), then we have η(v a ) � 0.2038. Due to the development of clean energy technology, this paper focuses only on emissions in urban traffic networks and ignores emissions from the power grid. Although both EVs and GVs contribute to the total system travel time, since only GVs contribute to the total system emissions, a bi-objective optimization program with the objective function (12) can be obtained: Let v � 2 Ω � v denote the set of Pareto-efficient solutions. The Pareto-efficient link flow pattern is defined as follows: Definition 1. Pareto-efficient flow pattern When a link flow pattern v � 2 Ω � v is Pareto-efficient optimal, neither the total system travel time T (v � ) nor the total system emissions E (v � ) can be reduced without increasing the other.
, and at least one of the inequalities holds: Basically, we calculate the emissions and travel time based on volume, and the objective is to minimize the total cost of both based on weighting factors for each. Since both T(v) and E (v) are convex, their Hessian matrices are both positive definite, and thus, the Pareto-efficient solution can be obtained by solving the weighting bi-objective mathematical problem (13): subject to constraints (1)-(4).
where m E w and m G w are Lagrange multipliers corresponding to constraints (1) and (2), respectively. It is evident that conditions (14)- (15) are consistent with the user equilibrium conditions of the EVs and GV users (7) and (8).
Definition 2. Differentiable marginal road pricing This research proposes a differentiable link-based marginal road pricing, differentiable road pricing or differentiable tolls for short, to manage congestion and emissions simultaneously in the urban networks with EVs and GVs. We investigate whether differentiable marginal pricing for EVs and GVs can induce the user equilibrium flow pattern into pareto optimum flow pattern by considering the congestion and emissions simultaneously. The optimal toll on a link is equal to the difference between the marginal social cost and the marginal private cost, which can internalize the user externalities including congestion part and emissions part, and thus achieve a Pareto-optimal efficient pattern in the network with EVs and GVs.
From the KKT conditions of bi-objective optimization program (1)-(4), (13), we can conclude that for any optimal solutions (f where m E w and m G w are Lagrange multipliers correspond to constraints (1) and (2), respectively. It is evident that conditions (14)-(15) are consistent with user equilibrium condition of EVs and GV users (7) and (8). According to the first order condition of the bi-objective program (1)-(4), (13), the differentiable road pricing for EV and GV can be expressed as (18) and (19):

Decentralization of differentiable tolls
In this section, we investigate whether the differentiable road pricing can induce user equilibrium into two-criteria pareto optimal in terms of congestion and emissions, and explore whether the differentiable can decentralize a Pareto-efficient flow pattern as a unique user equilibrium flow pattern.

can be decentralized as user equilibrium by nonnegative differentiable tolls.
Proof. Since Pareto-efficient solutions can be obtained by solving the bi-objective problem can be decentralized as a user equilibrium pattern by differentiable tolls.
First, consider the following linear programming (LP) problem: The dual problem of programming (26)-(29) is shown as follows: where l G a and l E a , a 2 A are Lagrange multipliers related to constraints (21) and (22) and m G w ðw 2 WÞ and m E w ðw 2 WÞ are vectors of multipliers associated with the constraints (23) and (24) Then, inequalities (27) and (28) can be rewritten as inequalities At the optimal points of the primal and dual problems, the following complementary slackness conditions will satisfy: Eqs (32) and (33) are alternatives of the user equilibrium conditions (7)-(8). Since we have proven that optimal solutions of the LP problem can be decentralized as a UE flow pattern by differentiable tolls, we investigate whether v � is the optimal solution of the LP problem through proof by contradiction.
Assume that v � is not the optimal solution of the LP programming (26)- (29), and v 0 is the optimal solution instead: v � 6 ¼ v 0 . According to the definition of Pareto-efficient flows, we can Similarly, e a (�) increases with any link flow pattern: Since inequalities (36) and (37) conflict with the definition of Pareto-efficient link flows (Definition 1), the assumption is false, and v 0 is not the optimal solution, implying that v � can be decentralized as user equilibrium. Thus, Proposition 1 can be concluded.
Since this bi-objective optimization problem involves with different vehicle types, and thus the link flow interactions among them cannot be ignored. Considering the interactions among link emission function e a (�) and link travel time function t a (�), the Jacobian of vector t a (v) is a matrix r v t(v) with the dimension of (2n) × (2n), and the Jacobian of vector e a (v) is a matrix r v e(v) with the dimension of n × 2n. Since only adjacent links have flow interactions, r v t(v) and r v e(v) are sparse matrices in most cases (Appendix 1).
Since the link flow interactions are asymmetric, this two-class mixed equilibrium link flow pattern cannot be obtained by the solution of an equivalent mathematical program but can be formulated as a variation inequality (VI) problem (38): Proof. First, we prove that Taking EV as an example, considering the user equilibrium condition (14), (39) holds. For GVs, inequality (40) holds. By summing inequality (39) and inequality (39) and summarizing the inequalities over all paths, r 2 R w , the OD pairs ω 2 W, and we obtain inequality (41): The right-hand side (RHS) of inequality (42) can be rewritten as: The left-hand side (LHS) of inequality(42) can be rewritten as (44): We prove that the solution of variation inequality (VI) (38) is a UE flow pattern. From inequality (38) and (42), we have: where the right-hand side (RHS) of inequality (45) is a constant and the left-hand side (LHS) is a linear program. Inequality (45) shows that � f is the optimal solution of the linear program, and the first-order conditions of the minimization program is (23)- (25) and (46)-(49).
Taking GV users as an example and c E rw are continuous and the solution sets in (38) are compact, the solution to variation inequality (38) always exists.
To analyze the stability and validity of the equilibrium, Proposition 3 is given to explore the uniqueness of the equilibrium for the two-criteria problem.
Proposition 2. Differentiable road pricing for EVs and GVs can decentralize the unique vehicle type-specific link flow pattern.

Proof.
Let v 1 and v 2 denote any two types of aggregate link flow patterns. Assuming that v 1 and v 2 are in user equilibrium, according to variation inequality (38), we have: After negating inequalities (HYPERLINK 50) and (51) summarizing them, we obtain: Since the link cost function of EVs and GVs is strictly increasing, we obtain: From inequality (52) and inequality (53), we obtain: Since v 1 and v 2 denote any two types of aggregate link flow patterns, we have c E a ðv 1 Þ 6 ¼ c E a ðv 2 Þ and c G a ðv 1 Þ 6 ¼ c G a ðv 2 Þ. According to Eq (54), we obtain v E1 a ¼ v E2 a and v G1 a ¼ v G2 a , thus Proposition 2 can be concluded.

Bound of the deviation factors
The solution set of problem (12) naturally forms a Pareto frontier (Fig 1), in which each point represents a Pareto optimum. Considering this two-criteria system under differentiable road pricing, we explore the gap deviation between the total system travel time (T) and total system emissions (E) in the Pareto optimum and their respective single-criterion-based system optimum.
By changing the percentage of EVs, T(v) remains constant, while E(v) varies within a range; this range can be bounded via Lemma 1: Lemma 1. For any feasible link flow pattern v under differentiable tolls, total system emissions E(v) can be bounded by total system travel time T (v): Proof. See in Appendix 2. Proposition 3. It holds that Proof. According to the definition of deviation factors, we have: By substituting v ME and v So into Lemma 1, we have: Thus, Proposition 3 can be concluded.

Numerical example
A numerical example is conducted. Travel demand d w is set as 1000, the link-based travel time function is t a (v a ) = 0.5 +0.15 (v a /c a ) 4 , and the link-based emission function is e a (v a ) = 0.2038 �t a (v a ) � exp(0.7962 �(l a /t a (v a ))). The percentage of EVs (P) are changed from 0.1 to 0.9, and the results for a � E and a � T are calculated. These results are summarized in Table 2. After obtaining the numerical results for the deviation factors, we selected three of these factors to plot the results in Fig 2, where the values of P are 0.1, 0.5 and 0.9. It is shown that the theoretical results of the deviation factors given in Proposition 3 can always bound the calculated results.

Bounding the efficiency gain of differentiable road pricing
To qualify how much this two-class two-criteria system will be improved after imposing differentiable road pricing, the efficiency gain of differentiable road pricing will be bounded by the concept of the price of anarchy. Here, the price of anarchy for total system travel time (POA T ) and the price of anarchy for total system emissions (POA E ) are used to evaluate the performance of differentiable road pricing while managing congestion and emissions, respectively.

Theoretical bound of price of anarchy
Definition 3. Price of anarchy for the total system travel time The price of anarchy for the total system travel time (POA T ) is the ratio of the total system travel time in user equilibrium state without tolls to that in pareto-efficient optimal state with differentiable tolls and is defined in Eq (60): where � v E and � v G represent the link flow patterns of EVs and GVs without tolls, andv E andv G represent the link flow patterns of EVs and GVs with differentiable tolls. The upper bound of POA T is given in Proposition 4. Proposition 4. Considering the upper bound of POA T , it holds that: Proof. See in Appendix 3.
The theoretical bound of POA T is only related to parameter β in BPR function (9), it is proved that the upper bound of POA T is independent on the urban network structure and percentage of electric vehicles. Definition 4. Price of anarchy for the total system emissions The price of anarchy for total system emissions (POA E ) is the ratio of total system emissions in user equilibrium state without tolls to total system emissions in pareto-efficient optimal state with differentiable tolls. It is defined in (62): where � v E and � v G represent the link flow patterns for EVs and GVs without tolls,v E andv G represent the link flow patterns for EVs and GVs with differentiable tolls.
Proposition 5. Given a non-overlapping network, it holds that: Prof. See in Appendix 4. The theoretical bound of POA E is related to the parameter β in BPR function (9), η max and η min in Eq (11). It is proved that the upper bound of POA T is independent on the urban network structure and percentage of electric vehicles.
We have algebraically derived out the general bound of POA T and POA E , in order to validate the validity of Proposition 4 and Proposition 5, we will calculate the numerical bound of POA T and POA E by conducting case studies based on the actual network and travel demand of two cities in China: Shenzhen and Lasa.

Input data.
The travel demand data used in this paper are based on the resident trip survey (The latest 2018) conducted twice per decade during the past 20 years and a traffic big data platform integrating mobile phone, public transport, taxi, navigation system, GPS system, and traffic volume data provided by Shenzhen Urban Planning Traffic Center (SUPTC).
The survey data include both 'roadside interviews' and 'home interviews, such as trip behavior, the nature and intensity of traffic, the freight structure, income, employment, etc. The surveys were well defined and divided into 'zones' so that the origins and destinations of trips could be geographically observed. The variables used [37] for the interviews are given in Table 3.
Basic travel demand data were collected by conducting a one-year home-interview survey during the afternoon peak hour (18:00-19:00). In transportation planning, fixed-demand models are often used to predict traffic flows during morning and evening peak hours because travel demand during rush hours mainly consists of commuters and is thus regarded as fixed.
In this research, total travel demand in Shenzhen are divided according to different administrative districts to estimate travel demand and generate the OD matrix ( Table 4). The original input data used for the OD matrix are shown in S2 Text. There are 491 OD pairs with fixed demand in Shenzhen (S2 File). The number of OD pairs in different districts is provided in Table 5.
The network topology and size are different in each district (Fig 3), and the number of links and nodes are summarized in Table 6.
In general, there are 6 road types in Shenzhen, including expressways, arterial roads, sub arterial roads, and approach and side roads. These are represented by different colors in Fig 3, which displays the basic network structure and road types in Shenzhen.
The BPR link time function given in (9) is used in the case study, and the values of parameters α, β and t 0 a of 6 different road types are shown in Table 7. The network structure in Shenzhen is shown in Fig 3. The mileage of each road type in 10 districts of Shenzhen is given in Table 8:

PLOS ONE
After preparing the data for conducting case study well, we will introduce the major procedure of the case study.

Major procedure of case study.
To calculate the numerical bounds of POA T and POA E , this paper conducted a case study using the software TransCAD 6.0. The Geographic Information System Developer's Kit (GISDK) [38] is a collection of software tools that come with TransCAD that makes it possible to feasibly revise the function in TransCAD software through integrating other programs. GISDK was used for the second development. The procedure includes the following steps: Step 1: Given the road networks of Shenzhen, categorize all the links over the networks into 10 groups according to the districts; Step 2: Use the OD table (Q) among 10 districts in Shenzhen to assign the total travel demand; Step 3: Using the objective function file, where vs2010 compiling and customized DDL will be used, differentiable tolls are added to the link impedance function.
Then, we conduct the first case study based on 10 districts in Shenzhen and evaluate the efficiency gain of differentiable road pricing in managing congestion and vehicular emissions. The code of the following procedure is provided in S1 Text.
The major procedure used to simulate the Shenzhen networks Step 1: Let � = 0.25; multiply Q by � to obtain the new ODtable �Q For � = 0.25:2 Step 2: Let P = 0; divide ΦQ into two groups: �Q 1 (EV) and �Q 2 (GV) For P = 0:1 Step 3: Initialize the network Step 4: Create the net file for route selection and the net set Step 5: Traffic is assigned based on a certain objective function Step 6: Update the network file Step 7: Calculate the total system travel time and total system emissions across all links in each district. Save the output figures and the calculation results Step8: P = P+0.1 End Step9: Let � = � +0.25 End 7×9 simulations are conducted in each district, and 630 total simulations are conducted across the 10 networks of Shenzhen.

Results analysis.
Results of POA T . Google Map is used to identify different road types in the road network of Shenzhen. The geographic information is input into a geographic information system (GIS), which is a software to compile the geographic information into a database and has the ability to illustrate the geographic information through graphs and charts. The input map data of Shenzhen are provided in S1 and S3 Files, S2 Text. Here, we choose the situation under original demand level (ϕ = 1) to display the distribution of total system travel time with and without road pricing. VOC represents the link travel cost varying from 0 (green) to 1.2 (red), which is represented by different colors. The number of link flows is represented by the line weight of the road link.
The results of case study show that both Tð� vÞ and TðvÞ will not change with P, which is consistent with Lemma 1: by changing the percentage of EVs, T(v) remains constant.
After implementing differentiable tolls, the link-based trip distribution will change considerably, as shown in Fig 4. When ϕ = 1, the total travel time changes the least in Nanshan and the most in the Longhua district. The total travel time across 10 districts without tolls Tð� vÞ is 69083 h, and that with differentiable tolls TðvÞ is 60767 h.
Since the travel demand data were collected by using a one-year home-interview survey during the afternoon peak hour, there are great differences in the congestion conditions during different time periods of the day. To better describe the variations in travel demand, the weighting factor of traffic demand level ϕ (multiplying the weighting factor by the original traffic demand) is introduced. Moreover, the percentage of EVs (P) among each OD pair with fixed demand is also introduced to simulate the mixed urban networks with both EVs and GVs.
In actual urban networks, the weighting factor of traffic demand level ϕ varies within a certain range. By changing ϕ within the range of 0.5~2, the results of link-based traffic flow solutions without and with differentiable tolls can be obtained (Table 9). It can be observed that both the total system travel time without tolls Tð� vÞ and the total system travel time with differentiable tolls TðvÞ increase with the travel demand. It can be observed that when the demand level is 1, Tð� vÞ is the largest in Baoan district (15452 h), and TðvÞ is the largest in Longgang district (13343 h); when the demand is doubled, both Tð� vÞ and TðvÞ are the largest in Longgang district, at 42728 h and 34791 h, respectively. The results of case study (S1 Data) show that both Tð� vÞ and TðvÞ do not change with the percentage of electric vehicles (EVs), which is consistent with Lemma 1. The results are shown in Table 9 After obtaining the calculation results of Tð� vÞ and TðvÞ, the results of POA T can be calculated according to Eq (60). By changing the weighting factor of demand level ϕ from 0.5 to 2, the results of POA T under different ϕ can be obtained. The results are summarized in Table 10.
As shown in Table 10, the upper numerical bound of POA T is 1.3393. When F is 0.5, 0.75, 1 and 1.5, the upper bound of POA T occurs in the Longhua district; when F is 1.25, the upper bound of POA T occurs in Yantian. When F is 1.75 and 2, the upper bound of POA T occurs in Baoan. The upper bounds of POA T under different F in each district is shown in Table 11.
The variation regulation of POA T in accordance with F is similar across 10 districts and first increases and then decreases, reaching the maximum value when the value of F is between 1.2 and 1.5. All 10 districts show a similar trend, which implies that there must be a mechanism that drives this regulation. It can also be concluded that the variation in POA T in terms of F is independent on the network structure. The results of case study show that the upper bound of POA T across 10 districts is 1.3393 ( Table 11).
The value of POA T in Nanshan is always the smallest among the 10 districts, implying that the travelers in this district are the least sensitive to the implementation of the differentiable road pricing scheme (Fig 5). Since Nanshan is not only the most economically powerful district in Shenzhen but also the center of scientific research, business, enterprise, and education, travelers commuting from Nanshan are less sensitive to the charging road pricing. If we focus on the maximum value of each curve, the maximum POA T is in Luo Hu, while Futian and Nanshan have the lowest values among 10 districts. These three districts are within the boundaries of the Shenzhen Special Economic Zone (SEZ), in which both the GDP and population rank among the top 3 in Shenzhen. It is showed that the efficiency gain of differentiable road pricing in the districts within the SEZ is less than that outside of the SEZ.
Results of POA E . When the percentage of EVs (P) continuously changes from 0.1 to 0.9, the distribution of vehicular emissions under scenarios without and with differentiable tolls are shown in Fig 6.  Here, we display the results of Eð� vÞ and EðvÞ under original demand level (F = 1). It is shown that the differentiable road pricing scheme could efficiently manage vehicular emissions.
Unlike the calculation results of the total system travel time, where both Tð� vÞ and TðvÞ do not change with P (Fig 4), the results in Fig 6 show that both Eð� vÞ and EðvÞ decrease with P. It is shown in Fig 6 (A) and S1 Data that Eð� vÞ will decrease from 1778098 g/h to 175342 g/h when P increases from 0 to 0.9. After implementing the differentiable tolls, EðvÞ will decrease from 1647737 g/h to 162181 g/h with an increase in P (Fig 6 (A)). If we fix P, the surface of Eð� vÞ is always above that of EðvÞ. Take a look at the case when P is 0.5, the result of Eð� vÞ and EðvÞ is 887407 g/h and 822227 g/h, respectively. When F = 1, both Eð� vÞ and EðvÞ change with P (Table 12). The other cases are shown in Fig 7 and S1 Data.
According to the calculation results displayed in Fig 7, there are more vehicular emissions in Luohu, Futian, Nanshan, Baoan and Longhua than in the other 5 districts. After implementing the differentiable road pricing scheme, total system emissions are reduced in each district to an extent. The upper surface shown in Fig 7 represents total emissions without tolls, and the lower curve surface represents total system emissions with differentiable tolls. It is shown that differentiable tolls can efficiently decrease total system emissions. There are more vehicular emissions in Luohu, Futian, Nanshan, Baoan and Longhua than in the other 5 districts. All 10 networks in Shenzhen display a similar trend: If we fix P, E will increase with ϕ; if we fix ϕ, E will decrease with P.
After obtaining the calculation results of Eð� vÞ and EðvÞ, the numerical results of POA E can be obtained according to formula (62), which is provided in S3 Data. The upper bound of  Table 13. The calculation results of POA E can be obtained by multiplying the original matrices by different ϕ. Here, ϕ changes from 0.5 to 2 and P changes from 0.1 to 0.9; 7×9 simulations are carried out in each district, and 630 simulations are carried out across the 10 networks in Shenzhen.  As is shown in S3 Data, the percentage of electric vehicles (P) has direct impacts on total system emissions (E) but has little impact on the value of POA E . To draw out the smooth surface graph shown in Fig 8, each matrix with 7×9 dimensions can be expanded into 350×450 dimensions through the cubic spline interpolation method. The P-F-POA E surface graph ( Fig  8) is plotted to show the relationship among P, F and POA E .
The similarities across the 10 networks in Shenzhen imply that there must be general mechanisms to explain the variations. However, the individual data of each district show that the values of POA E are very different in each district, and they are too varied to be statistically significant. For example, Fig 8(d) has only 2 peaks, while the figures for other districts have 3 peaks considering the relationship between POA E and F. Above all, the variation trends of POA E regarding F and P tend to be similar across 10 networks in Shenzhen.
We also have some interesting findings from the simulation results: the variation trends of POA T and POA E in terms of ϕ are similar across the different networks, implying that these  values are independent of the network topology. The efficiency gain of differentiable road pricing in managing system total travel time is less in the districts within the Shenzhen Economic Special Zone (SEZ) than in the other districts, and the efficiency gain of differentiable road pricing in managing system total emissions is larger in the districts within the SEZ than in the other districts, which implies that differentiable road pricing is more effective for managing congestion outside the SEZ and managing emissions within the SEZ.

Feasibility validation. 1) Confirmation experiment. Proposition 4 and
Proposition 5 indicate that the bound of price of anarchy for total system travel time (POA T ) and (POA E ) are independent on the network structure, and these regulations observed on the 10 networks of Shenzhen also exist on networks of other cities. In order to validate the theoretical bound and generality of these insights and findings, we conducted another case study based on a city in China, Lasa. This case study is conducted based on the same procedure introduced before.
• Results of Total system travel time

PLOS ONE
Likewise, we only choose the scenario with original demand level (ϕ = 1) to display the link-based traffic flow pattern without differentiable tolls � v, and link-based traffic flow pattern with differentiable tollsv, respectively, and other results can be found in S2 Data.
Similarly, the total system travel time decreases a lot after implementing the differentiable road pricing scheme. Then we display the results of only Eð� vÞ and EðvÞ for the original demand level (ϕ = 1). After obtaining the calculation results of Tð� vÞ and Tð� vÞ, the results of POA T can be calculated according to Eq (60). By changing the weighting factor of demand level F from 0.5 to 2, the calculation results of POA T under different F can be obtained. The results are summarized in Table 14.
As is shown in Fig 10 and Table 14, the upper numerical bound of POA T in Lasa is 1.2653 when the value of F is 1.5. Similar to the calculation results of Shenzhen (Table 10), the numerical bound of POA T in Lasa is bounded by the theoretical bound in Proposition 4. The variation regulation of POA T in Lasa is similar to that in Shenzhen, which validates that the bound of POA T is independent on the urban traffic network structure. and other results can be found in S2 Data. • Results of total system emissions In Fig 11, if we fix the percentage of EVs (P) and compare Fig 11(A) with Fig 11(B) respectively, it is showed that the total system emissions in B are always smaller than that in A, indicating that the differentiable road pricing scheme could efficiently manage vehicular emissions. If we observe Fig 11(A) and 11(B), total system emissions decrease with P in both cases.
After obtaining the calculation results of Eð� vÞ and EðvÞ, the numerical results of POA E can be obtained according to formula (62), which is provided in S4 Data. The upper numerical bound of POA E in Lasa is 1.2678 when F is 1.25 and P is 0.9. the variation trends of POA E regarding F and P tend to be similar with Shenzhen (Fig 12).
2) Sensitivity analysis on theoretical bound of POA T . As is shown in Proposition 4, theoretical bound of POA T is only related to the parameter β in BPR function, the value of parameter β is determined by the road type ( Table 7). The calculation results show that the numerical upper bound of POA T in Shenzhen is 1.3393, and the numerical upper bound of POA E in Lasa is 1.2653. According to Table 7, the value of β is decided by the road type and can be set within the range of 2.1-2.6. Then we substitute the value of β into Proposition 4 and conduct parameter analysis on the values of POA � T (Table 15). Sensitivity analysis on parameter β (Table 15) (Table 16). The calculation results show that the numerical upper bound of POA E in Shenzhen is 1.3393, and the numerical upper bound of POA E in Lasa is 1.2678. According to Table 7, the value of β is decided by the road type and can be set within the range of 2.1-2.6. Sensitivity analysis on parameter β (Table 16) show that the numerical values of POA T can always be bounded by POA � E , this analysis validate the effectiveness of Proposition 5. From the experimental results of both cities, the variation trends of POA T and POA E are similar across the different networks, implying that the bound of price of anarchy is independent on the network structure and percentage of electric vehicles.

Conclusions and future study
In this paper, we theoretically and empirically bound the efficiency of differentiable road pricing for EVs and GVs. Considering the asymmetric link flow interactions between EVs and GVs, Proposition 1 and Proposition 2 prove that the differentiable road pricing can decentralize a Pareto-efficient flow as the unique user equilibrium. The deviation gap between the single-criterion-based system optimal on Pareto frontier is bounded by two deviation factors and given in Lemma 1 and Proposition 3. Finally, the theoretical bound of POA T and POA E is given in Proposition 4 and Proposition 5. In order to validate the effectiveness of theoretical bound, we conduct case studies based on two cities of China: Shenzhen and Lasa. Empirical results show that the numerical bound of POA T is 1.3393 in Shenzhen and 1.2653 Lasa, respectively. The numerical bound of POA E is 1.4414 in Shenzhen and 1.2678 in Lasa, respectively. After sensitivity analysis on different parameters related to the bound of POA, we  validate the effectiveness of Proposition 4 and Proposition 5. Moreover, the variation regulation of POA experiences a similar trend in different networks, which implies that there must be a mechanism driving this regulation. From both theoretical derivation and empirical analysis, it can be concluded that the bound POA T in is independent on the network structure and percentage of electric vehicles. To develop this research further, the model will be extended to the multi-class model with elastic demand, where users' various value of time and travel demand will be considered. This model can also be extended to include other air pollutants except for carbon monoxide (CO).
deviation gap in this two-criteria system, we define two deviation factors: where α T (v) measures the deviation gap between the total system travel time of any Pareto-efficient flow pattern T(v) and the minimum system travel time T min , and α E (v) measures the deviation gap between the total system emissions between any feasible flow pattern E(v) and the minimum system emissions E min . For any Pareto optimal link flow v, a � T and a � E represent the upper bounds of α T (v) and α E (v), respectively; a T ðvÞ � a � T ; a E ðvÞ � a � E :a � T and a � E are defined in Eqs (76) and (77): After investigating the bound of a � T ðvÞ and a � E ðvÞ, we obtain Lemma 1. Then, Proposition 3 can be easily obtained through Lemma 1.
Here, 0/0 by convention. According to the KKT condition of program (21)-(25), we can obtain that:v