Joint deployment of charging stations and photovoltaic power plants for electric vehicles

Transportation Research Part D 79 (2020) 102247

Contents lists available at ScienceDirect

Transportation Research Part D journal homepage: www.elsevier.com/locate/trd

Joint deployment of charging stations and photovoltaic power plants for electric vehicles Zhixiong Luoa, Fang Hea,b, Xi Linc,b, Jianjun Wud, Meng Lic,b,

T

⁎

a

Department of Industrial Engineering, Tsinghua University, Beijing 100084, PR China Tsinghua-Daimler Joint Research Center for Sustainable Transportation, Tsinghua University, Beijing 100084, PR China Department of Civil Engineering, Tsinghua University, Beijing 100084, PR China d State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong University, Beijing 100044, PR China b c

ARTICLE INFO

ABSTRACT

Keywords: Electric vehicles PV power plants Charging stations Time-dependent pricing Surrogate-based optimization

The focus of this study is to jointly design charging stations and photovoltaic (PV) power plants with time-dependent charging fee, to improve the management of the coupled transportation and power systems. We first propose an efficient and extended label-setting algorithm to solve the EV joint routing and charging problem that considers recharging amount choices at different stations and loop movement cases. Then, a variational inequality problem is formulated to model the equilibrium of EV traffic on transportation networks, and an optimal power flow model is proposed to model the power network flow with PV power plants and optimally serve the EV charging requirements. Based on the above models for describing system states, we then formulate a model to simultaneously design charging stations, PV plants, and time-dependent charging fee. A surrogate-based optimization (SBO) algorithm is adopted to solve the model. Numerical examples demonstrate that the proposed SBO algorithm performs well. Additionally, important insights concerning the infrastructure design and price management of the coupled transportation and power networks are derived accordingly.

1. Introduction Electric vehicles (EVs)1 are gaining increasing attention worldwide as a potential approach to reducing carbon emissions from transportation systems (Yu et al., 2018). Nevertheless, the energy sources of EVs could still significantly contribute to environmental deterioration. Based on some well-to-wheel emission analyses, EVs may discharge more greenhouse gases compared to hybrid EVs and even conventional oil-fueled vehicles in areas mainly powered by coal-fired power plants (Elgowainy et al., 2009). Thus, it is recognized that the popularization of EVs should be synchronized with the deployment of renewable energy generation. For example, the government in China is simultaneously making significant investments in EVs and renewable energy. According to the government’s plan, 4.8 million charging piles and over 12,000 charging stations are to be constructed by 2020 (State Council of the People’s Republic of China, 2015). Meanwhile, the renewable energy capacity reached 706 million kW in September 2018, in which photovoltaic (PV) energy accounted for 165 million kW, which is expected to continue rising (Xinhua Net, 2018). Among the existing types of renewable energy, solar has a higher potential for supporting the daily movements of EVs. There are several reasons that solar energy is well suited to provide energy support for EVs. First, as pointed out by Denholm et al. (2013), the

Corresponding author at: Department of Civil Engineering, Tsinghua University, Beijing 100084, PR China. E-mail address: [email protected] (M. Li). 1 For the remainder of this study, EV only refers to battery-powered EVs; other types of EVs will be particularly identified. ⁎

https://doi.org/10.1016/j.trd.2020.102247

Available online 02 February 2020 1361-9209/ © 2020 Elsevier Ltd. All rights reserved.

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

charging demands of EVs match well with the generation of solar energy. Second, the power plants for solar energy are spatially distributed, and they can be installed near the charging stations to serve the charging requirements of EVs. Third, distributed PV generators can be equipped with advanced grid-connected inverters to support reactive power control technology, thus keeping the power quality from being deteriorated by large and unevenly distributed charging demands (Zhang et al., 2018). From such a viewpoint, constructing PV plants to provide energy support for charging stations and enable travel by EVs is a promising solution to truly alleviating carbon emissions from transportation systems. However, even with the advantages presented above, jointly considering the planning of charging stations and PV power plants remains a significant challenge. The main difficulties can be summarized as follows: (a) The interactions among EV drivers’ choices, network traffic, and power flow are relatively complicated. Generally, the locations and capacities of charging stations will exert influences on the EV drivers’ equilibrated routing and charging choices, thus reshaping the traffic on transportation networks. Moreover, the optimal power flow on power networks is determined not only by the locations of charging stations and PV generators, but also by the charging requirements at each charging station. Sophisticated quantitative tools are required to capture such interactions. (b) The traffic demands of EVs and the power generation of PV plants both exhibit high variances in the time horizon, e.g., the daytime/nighttime differences or weekday/weekend differences (Xydas et al., 2016; Moon et al., 2018; Mellit and Pavan, 2010). On one hand, such variations bring difficulties in the planning process by introducing heterogeneities; on the other hand, to manage the demand and supply of energy provision in different periods more effectively, some time-dependent managing schemes (e.g., time-dependent charging fee) could be incorporated into the system (He et al., 2016), which further complicates the problem. (c) Regardless of how well the infrastructure systems and time-dependent managing schemes are designed, the power generated by PV plants may not be exactly consumed in situ by the EV charging activities. There are always cases where redundant power should be sold or additional power should be purchased. Such power transmissions should also be modeled adequately. In recent years, great efforts have been exerted to tackle the issue of deploying public charging stations on transportation networks. The methodologies can be roughly classified into the following five approaches: clustering, flow-capturing, simulation-based, activity-based, and equilibrium-based. The clustering approach utilizes traffic flow data to conduct a cluster analysis, based on which a location model is developed (Ip et al., 2010; Momtazpour et al., 2014). The flow-capturing approach assumes full knowledge of the origin–destination (O–D) information and traffic flow distribution, and then determines the locations of charging stations by minimizing the travel cost under the feasible path constraints (Hodgson and Berman, 1997; Wu and Sioshansi, 2017). Although they are simple, the two approaches above cannot capture the influence of infrastructure locations on drivers’ behaviors. The simulation-based approach mainly evaluates the impact of charging stations on EV movements with some behavioral rules, and such rules can be drawn from real data (Dong et al., 2014; Andrews et al., 2013; Li et al., 2017). Similarly, activity-based approach studies the routing and scheduling problem of individual EVs (Recker and Kang, 2010; Kang and Recker, 2014). Finally, the equilibrium-based approach incorporates the route choice behavior of drivers into a traffic assignment paradigm, and bi-level optimization models are built to determine the optimal planning of charging stations (He et al., 2013, 2014, 2015; Wang et al., 2019). Among these planning models for charging facilities, only a few have taken the power network into consideration. As an example of this workflow, He et al. (2013) incorporated the transportation network equilibrium with a coupled power network to deploy public charging stations for plug-in hybrid EVs. In their work, a direct current power flow model was utilized to describe the optimal power flow (OPF) of the power network, and a locational marginal price was encapsulated to determine the associated electricity price. However, the above researches don't include distributed energy system. From the perspective of power system, the basic technique is optimal power flow (OPF) equation. Dommel and Tinney (1968) first gave conventional computational solutions for the OPF with an alternating current (AC) network. Gan et al. (2015) turned the AC OPF equations into a second-order cone program (SOCP) by relaxing one constraint in the power flow equation, and SOCP can then be solved by commercial optimizers. With the development of smart grid technology, i.e. microgrids, distributed generator, distribution storage and smart grid control technology, it is possible to better plan, operate and control the power system (Zhang et al., 2015; Thomas et al., 2018; Pena-Bello et al., 2019); however, most studies in these fields simply focus on the power system itself without considering its interactions with transportation system. Some other researchers noticed the great benefits of jointly considering EV charging and photovoltaics generators, but they mainly focus on the operation at single places, i.e. individual houses (Munkhammar et al., 2015; Salpakari et al., 2017), commercial buildings (Tulpule et al., 2013), and charging stations (Figueiredo et al., 2017). As a result, those studies oversimplify the description of transportation system, and they ignore the system optimization at network scale. Only a few paper studied the deployment of power generators and charging facilities on a road network in recent years. For instance, Zhang et al. (2017) constructed a mixed-integer SOCP (MISOCP) to decide the optimal planning of fast charging stations, including their sites and sizes. A further extension of this work was also proposed by Zhang et al. (2018) with an accelerated generalized benders decomposition to decide the location and size of charging stations and PV power plants. However, the behavior of driver choices was simplified for the solvability of the optimization model. Based on the above reviews, we observe that the current literature has rarely focused on investigating the joint deployment of charging stations and PV plants with a comprehensive capture of the interplay between the coupled transportation and power networks under time-variant traffic demands and energy supplies, and it has not proposed time-dependent economical schemes for better managing the coupled networks. This study is aimed at filling this void. Specifically, to achieve the targeted research goals, we conducted the following tasks as the major components of this study: 2

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

(a) A novel extended label-setting algorithm to solve the joint routing and charging problem (JRCP) for EV drivers considering driving time, recharging time, recharging price, and queuing time at charging stations is developed. Specifically, the algorithm can simultaneously optimize the route and recharging amount choices. Moreover, it can capture the possible loops in drivers’ route choice, which are basically ignored in previous research. (b) A variational inequality (VI) formulation is proposed to model the equilibration of EV drivers’ routing and charging choices on transportation networks considering both traffic congestions on roads and queuing delays at charging stations due to high charging demands. A column-generation approach, coupled with the proposed routing model, is developed accordingly to solve the equilibrium model. (c) An MISOCP model is constructed to model the OPF with an AC power network. Power selling and purchasing are incorporated into the model. (d) A joint optimization model for the design of charging stations, PV power plants, and time-dependent charging fee is formulated with the above equilibrium model and OPF model. Such optimization model is difficult to solve by conventional optimization techniques; thus, a surrogate-based optimization (SBO) algorithm is incorporated to find high-quality solutions within an acceptable computational time. (e) Extensive numerical examples with various test networks are conducted. The algorithm performance, optimal planning, and pricing along with some managerial insights are illustrated and discussed. The remainder of this paper is organized as follows. Section 2 presents our basic modeling considerations and the coupled transportation and power systems, which consist of JRCP of EVs, traffic equilibrium model, and OPF model. Section 3 describes the joint optimal deployment (JOD) model and introduces the SBO framework to solve it. Section 4 illustrates the results of numerical examples, and finally, Section 5 concludes the study. 2. Model for coupled transportation and power system 2.1. Basic considerations In this research, EV drivers, transportation agency, power agency, and city manager are considered as stakeholders. Their relationship is shown as follows (see Fig. 1). In this system, the transportation agency decides the deployment of charging stations including the location and capacity (number of charging piles) of each charging station, as well as the time-dependent charging fee at each charging station for better managing the charging activities. Its object is to minimize the transportation system cost, including travel time, queuing time and charging time. EV drivers are assumed to make their routing and charging choices to minimize their total costs, including travel time on road, queuing time at charging station, and charging cost at charging station. And for the drivers charging at a charging station, the charging cost consists of charging time, charging fee, and electricity fee. Specially, the electricity price is assumed to be pre-determined in this research due to the government regulation, so it is not manipulatable. Then, the power agency decides the deployment of PV power plants, including the locations and sizes, and the power transaction with outside power network, i.e., selling power during periods with high solar output or purchasing power during other periods. The objective of power agency is to minimize the cost of power system, including power transaction cost with outside power network and the penalty of undesired voltage variance, because more stable voltage means higher power quality (Zhang et al., 2018). The city manager works as the coordinator of transportation and power agencies to minimize the cost of the coupled transportation-power system. For the modeling in detail, we consider a coupled transportation and power distribution network system as shown in Fig. 2 below. The transportation network is represented by a directed graph G (N , A ) , where N is the set of nodes and A is the set of links; the charging stations are distributed on some nodes in the transportation network. EV travel demands are generated from one set of nodes to another, and the set of O–D pairs is denoted by W . The EV drivers make routing and charging decisions based on their experience

Fig. 1. Relationship between stakeholders. 3

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Fig. 2. Illustration of a coupled transportation and power network.

including travel time, charging time, charging fee, electricity fee, and queuing time at charging stations. The travel time on a link is a monotonic and convex function of the associated traffic volume. The charging time at a charging station is a linear function of the charging amount (Hõimoja et al., 2012; Pourazarm et al, 2016). The electricity fee is a linear function of the charging amount, and the rate is electricity price. The charging fee is a linear function of the charging amount where the station-specific rate is determined by our succeeding optimization model. The queuing time of a charging station is a function of the associated total charging demand (the functional form will be introduced later in subsection 2.3). Besides, drivers are assumed to be fully aware of the information, thus a classic traffic equilibrium state is formulated in transportation network. The power distribution network is represented by an undirected graph G (M , L ) , where M is the set of distribution buses and L is the set of distribution lines. All the distribution buses except for the substation are candidate locations for building PV power plants. The substation is the root bus, indexed by zero, which transmits the energy between the distribution network and outside network. For any distribution bus m M , Nm denotes the set of nodes in the transportation network, whose electricity load is supplied by bus m M . For example, the blue dotted line indicates that the electricity loads of node 5 and node 9 are supplied by bus 3 (see Fig. 2). For any node n N , the electricity load consists of four parts: EVs charging load, commercial load, industrial load, and residual load. Note that both the EV travel demands and rates of power generation of PV plants are time-dependent. To capture their time variances, we set a pool of scenarios q Q , and each scenario denotes a particular time interval. The objective of the joint design and pricing model in this study are to maximize the long-term benefits over all the divided time intervals. In the following, we successively introduce an extended label-setting algorithm that solves the JRCP under given traffic and price settings, the traffic equilibrium model that captures the overall traffic state, and the location model of PV generators with OPF given the electricity loads of charging stations. For notational brevity, all subscripts on scenario s are ignored in the JRCP and traffic equilibrium model. 2.2. Joint routing and charging problem This subsection introduces the algorithm for solving the JRCP of EVs. Given travel time, charging cost (including charging time, electricity fee and charging fee) per charging amount, and queuing time at charging stations, the algorithm optimizes an EV’s routing and recharging amount choices at different charging stations to fulfill its trip without running out of battery and minimize its generalized travel cost. Similar routing problems are comprehensively investigated in the literature (e.g., He et al., 2014, 2015; Xu et al., 2017). However, as illustrated by the following example, there exists a case of loop movements that previous methods may not be able to handle appropriately, as most of them try to avoid a loop and take it as an error, i.e., label-correcting algorithm (Ahuja et al., 1993; He et al., 2014, 2015). To the best of the authors’ knowledge, there is still no research on solving the JRCP of EVs that considers loop movements. As shown in Fig. 3, considering an EV with a maximum range of 160 km and an initial range of 60 km starting at the origin (node 1); the vehicle cannot reach the destination (node 4) without charging. The vehicle can charge at node 3 for a range of 40 km and then drive to node 4, in which case node 1 and node 2 are both bypassed two times. We observe that “loops” occur in the path, and the occurrence of such loops nullifies the use of a link-based shortest-path formulation as described by He et al. (2014). Origin

Destination

4

1

80 km

Charging Station 5 km 5 km

Fig. 3. Illustration of an EV path consisting of loops.

4

2

5 km 5 km

3

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

To tackle the problem of possible loops shown above, we propose an extended label-setting algorithm (Ahuja et al., 1993) to solve the JRCP as follows. Step0: Step1: Step2: Step3: Step4: Step5:

Initialize three empty sets: R 1, R 2, R3 Start from the origin nodes, generate the its label(s), i.e., lo . Let R 1 = {lo} For labels in R 1, generate new labels for all their directly subsequent nodes. Then let R2 being the set of those new labels. R3 R1, and perform the dominance test for all labels in R2 using R3 R 2 . Delete all the dominated labels in R2 and let R 1 R 2 , R2 Let R3 If R 1 is an empty set, go to step 5, otherwise, go to Step 2. Choose the minimum-cost label corresponding to the destination node in R3 and figure out the EV routing using backward induction method.

.

To illustrate the JRCP algorithm explicitly, we show the definition of a label, the label generation, and dominance rules. For brevity, in this algorithm, we take the unit of state of charge (SOC) as “km” and the unit of charging price as time in “min.” We define a label as a vector lc = {Np, Nc , Tc, Ec , Wc, Sc, tc, ec } , in which Nc is the current node and Np is the previous node. If Nc is the origin node, we let Np = 0. Tc , Ec , Wc , and Sc are the travel cost, SOC, queuing time, and unit charging cost at the previous charging station with positive charging amount, respectively; tc and ec are the travel cost and electricity consumption from the last charging station to the current node, respectively. If a vehicle does not charge before Nc , we let Tc , Ec , Wc , and Sc to be zero, and tc and ec to be respectively the travel cost and electricity consumption from the origin to Nc . The generating and dominance rules are respectively introduced in parts 2.2.1 and 2.2.2 below. 2.2.1. Label generating rule In Step 1, at the origin node, the label is lo = {0, o, 0, 0, 0, 0, 0, 0} . If the origin node is a charging station, we generate another label lo' = {0, o, 0, Einitial , Wo, So, 0, 0} , where Einitial denotes the initial SOC of a vehicle traveling from origin "o" to destination "d . "Wo and So denote the queuing time and unit charging cost at the origin, respectively. In Step 2, without loss of generality, we show the label generating rule for one of the labels in R1 and one of its directly subsequent node, i.e., lc = {Np, Nc , Tc, Ec , Wc, Sc, tc, ec } and Ns , where (Nc , Ns ) A , and for brevity we define a = (Nc , Ns ) . Let Esafe and Emax denote the minimum and maximum allowable SOC, respectively. Then, if the vehicle does not charge before Ns , it can reach Ns under the condition Einitial (ec + ea) Esafe . On the other hand, if the vehicle charges before Ns , it can reach Ns under the condition Emax (ec + ea) Esafe . When the vehicle can reach Ns , we generate label(s) based on the following four cases. Case I:. When Ns is a normal node (not a charging station), then generate a label ls,1. Case II:. When Ns is a charging station and the first recharging station for the vehicle, then generate labels ls,1 and ls,2 . Case III:. When Ns is a charging station, but not the first recharging station for the vehicle with Ss ls,3. which means that the vehicle charges just enough at the last charging station to get to Ns .

Sc , then generate labels ls,1 and

Case IV:. When Ns is a charging station, but not the first recharging station for the vehicle with Ss > Sc (contrary to case III), then generate labels ls,1 and ls,4 , which means that the vehicle charges to Emax at the previous charging station and then charges at Ns . The details of labels ls,1, ls,2 , ls,3, and ls,4 are listed as follows. ls,1 ls,2 ls,3

{Nc , Ns , Tc , Ec , Wc, Sc, t + ta, ec + ea} {Nc , Ns , t + ta, Einitial ec ea , Ws , Ss , 0, 0} {Nc , Ns , Tc + Wc + tc + ta + Sc max{Esafe + ec + ea

ls,4

{Nc , Ns , Tc + Wc + tc + ta + Sc (Emax

Ec ), Emax

Ec , 0}, Esafe , Ws , Ss , 0, 0} ec

ea, Ws, Ss, 0, 0}

From the label generating rule, we can easily know that this rule can find the path with loops. The difficult part of the above generating rule is that if the unit charging price of the former charging station is higher or equal to that of the latter charging station, then the vehicle charges just enough amount of electricity at the former charging station to reach the latter charging station; otherwise, the vehicle charges to full SOC (maximum SOC) at the former charging station. We name this rule HMLM (Higher cost in former charging station - Minimum charging amount - Lower cost in former charging station - Maximum charging amount). The following proposition demonstrates that at least one optimal solution to the JRCP satisfies HMLM. Proposition 1.. The set of optimal solutions to the JRCP contains one satisfying HMLM. The proof of proposition 1 is presented in Appendix B. For readers’ comprehension convenience, two toy examples are provided in Fig. 4. For an O-D pair 1–4, the maximum SOC, initial SOC, and safe SOC are equivalent to 120 km, 120 km, and 20 km, respectively. In instance a, the charging cost at charging station 2 and 3 are 1 unit/km and 2 units/km, respectively. Instance b is the opposite. In instance a, the optimal policy is to charge 60 km (reaching the maximum SOC) at charging station 2. In instance b, the optimal policy is to charge 20 km (minimum amount to get to charging station 3) at charging station 2. The SOC along the routing is shown in the blue line, which is consistent with “HMLM”. 2.2.2. Label dominance rule Recall that a label at node Nc is defined as lc = {Np, Nc , Tc, Ec , Wc, Sc, tc, ec } . The SOC at node Nc ranges from [max{Esafe, Ec ec}, Emax ec ] when the vehicle has charged before; otherwise, the SOC is Einitial ec . The following statement defines 5

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Fig. 4. Toy example for HMLM rule.

a core terminology, i.e., dominance, used in JRCP. Definition 1.. At node Nc N , label lc' is said to dominate label lc if for any label at the destination node Nd originating from lc , we can always find another label at Nd originating from lc' with smaller or equal total travel cost. Clearly, in the JRCP process, we can remove those labels dominated by others without taking the risk of losing the optimal one. The proposition below then provides a specific rule for label dominance. Proposition 2.. Sc' Sc , ec' Ec' ec

lc' = {N p' , Nc , Tc' , Ec' , Wc' , Sc' , tc' , ec' } dominates Ec , and ec' ec .

lc = {Np, Nc , Tc, Ec , Wc, Sc, tc, ec }

if

Tc' + Wc' + tc'

Tc + Wc + tc ,

The proof of proposition 2 is presented in Appendix B. Now, with Propositions 1 and 2, we can guarantee that (1) the extended label-setting algorithm is capable of generating the optimal solution, and (2) the optimal solution would never be dominated in the dominance test. Thus, it is sufficient to conclude that the algorithm can output one optimal solution to the JRCP upon termination. 2.3. Traffic equilibrium of EVs With the joint routing and charging protocol presented in the previous subsection, we are ready to propose the traffic equilibrium model for EVs. Let W denote the set of an O–D pair, Dw the demand of O–D pair w W , and Pw the feasible path set for O–D pair w W ; P = w W Pw . Note that a feasible path is one that a vehicle can reach to the destination with/without charging through the path. Here, we define that the path is the joint link choices and charging station choices, with the charging amount choice at each charging station being determined based on the JRCP mentioned above. Then, the feasible domain of path flow V can be formulated as

V: p Pw

f pw

f pw = Dw

0

w

w

W, p

W

(1) (2)

Pw

denotes the flow on path p Pw for O–D pair w W . Constraint (1) ensures that the sum of the path flow between O–D Here, pair w W is equal to the corresponding demand Dw . Constraint (2) guarantees that all path flows are non-negative. Within the feasible domain V , the traffic equilibrium condition (TEC) can be formulated as a complementary system below. TEC:

f pw

cpw

w cmin

f pw (cpw

0 w cmin )

w =0

W, p w

(3)

Pw

W, p

(4)

Pw w and cmin

denotes the cost of path p Pw, whilew W denote the minimum path cost of O–D pair w W . Constraint (3) Here, w ensures that cmin is the minimum path cost of O–D pair w W . Constraint (4) is a complementary constraint, meaning that for any path p Pw whose cost is larger than the minimum path cost of O–D pair w W , the path flow is zero. To solve the traffic V by solving the VI and then prove the equivalence of the VI equilibrium problem, we formulate the equilibrium as the solution f and TEC. VI:

cpw

w W

p Pw

cpw f pw

f pw

0

In the VI, cpw is the path cost for p

f

V

(5)

Pw under path flow f 6

V . For any given path flow f

V , the formulation of the

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

corresponding cpw is

cpw =

ta (va ga )

w p, a

+

a A

w p, n

+ Epw, n (sn + pn,0 + pn ))

(6)

w W

p Pw

4

va ga

ta (va ga ) = ta0 1 + 0.15

va =

(wn (vn, En sn, kn ) n N

a

w w p, a f p

a

(7)

A

sn vn diag (En EnT ) 2vn En

wn (vn, En sn, kn) =

A

(8)

6 i=0

i, kn (sn vn

En )i , vn En > 0

0,

(9)

vn En = 0

In the above formulation, ta (va ga ) is the travel time of link a A , which is a strictly increasing function of the flow va on the link with capacity ga as a parameter. pw, a denotes that link a A is utilized by path p Pw, w W if it is equal to 1 and 0 otherwise; wn (vn, En sn) is the queuing time of node n N , which is a function of two column vectors; (vn, En ) , sn , and kn respectively denote the traffic flows and charging demands, the charging time for unit charge, and the number of charging piles at node n ; pw, n is equal to 1 if node n N is utilized by path p Pw, w W and 0 otherwise; Epw, n is the charging amount of path p Pw, w W at node n N , and it is a predetermined parameter because the charging amount can be determined given the path; pn,0 and pn are the electricity price and unit charging fee, respectively, for a unit charging amount at node n N . Without loss of generality, we adapt the form of the Bureau of Public Roads function in this research in constraint (7). Constraint (8) is used to calculate the flow of link a A . As noticed from Section 2.2, the charging amount at each node is a general distribution. We model the queue at charging stations as the M G k systems and incorporate the average waiting time approximation proposed by Kimura (1996):

tk, G =

Var (G ) + E (G )2 1+ 2E (G )(k E (G ))

k 1

(k

i=1

where G is the distribution of service time, charging station, we have v

1

1) ! (k E (G )) i! ( E (G )) k i

(10)

is the arriving rate, and k is the number of servers. Taking this into the queuing of

diag (E EnT )

sn2 n v n Var (G ) + E (G ) 2 E (G 2 ) n 1 = = v E 2E ( G ) 2E ( G ) 2sn nv n n 1

=

sn vn diag (En EnT ) 2vn En

(11)

where diag () is a function for obtaining the diagonal elements of a determinant and turning them into a column vector. For the rest of this approximation (k

1 E (G ))

1+

k 1 (k 1) ! (k E (G )) i=1 i ! ( E (G ))k i

1

, we take a polynomial function (

6 i=0

i, k (

E (G ))i ) to fit it in order to

improve the solution efficiency. Thus, the formulation of queuing time at node n is as shown in constraint (9). Next, we analyze the mathematical properties of the VI and prove that it is equivalent to the traffic equilibrium. There exists at least one solution for the VI, as the feasible path flow domain V is a compact set and the mapping of VI is continuous (Harker and Pang, 1990). However, there may not exist an equivalent convex programming formulation, as the mapping of VI is asymmetric (Wang et al., 2019). The equivalence of VI to the equilibrium condition is verified as follows. Proposition 3.. The VI proposed in (5) is equivalent to the equilibrium conditions (3)–(4).

We provide the proof of proposition 3 in Appendix B. We then adopt the asymmetric VI solution approach proposed by Aghassi et al. (2006). This approach makes use of the linear dual property and turns the VI into a nonlinear programming (NLP) problem. The formulation is as follows: NLP:

minz = f , ,c

s. t .

w W

p Pw

p Pw

cpw f pw

f pw = Dw

w W

w

w Dw

W

(12)

f pw

0

w

W, p

Pw

(13)

w

cpw

w

W, p

Pw

(14)

In the NLP, cpw is the path cost for drivers from O–D pair w W and choosing path p Pw and cpw is determined by constraints (6)–(9). w denotes the dual variable corresponding to constraint (12), whose physical meaning is the minimum path cost for O–D pair w W . If the solution results in z = 0 , the solution satisfies the VI and then a traffic equilibrium is obtained. Because the number of paths is intractable in a realistic network, we employ the column generation framework (Lawphongpanich and Hearn, 1984) to gradually expand the usable set and solve the NLP within the restricted path set during each iteration. Let Pw¯ denote the restricted Pw ; then, the restricted path set for all O–D pairs is P¯ = w W Pw¯ . We define the NLP under a path set for O–D pair w W and Pw¯ restricted path set as follows: 7

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

NLP(P¯ ):

minz =

s. t .

cpw f pw

¯ p Pw

w W

f , ,c

f pw = Dw

¯ p Pw

w Dw

w W

w

W

(15)

f pw

0

w

W, p

Pw¯

(16)

w

cpw

w

W, p

Pw¯

(17)

The overall framework for solving the NLP problem with a column generation framework is as follows. CG-NLP: Step0:

Initialize the network, set path flow f 0 be 0. Set the initial path set P¯ = . Set the iteration count n = 0 . Determine all traffic states including link travel time and average queuing time with f n . Solve JRCP with current traffic states for all O-D pairs and obtain the shortest path set P¯n . If (P¯ P¯n) = P¯n , terminate the CG-NLP. Then the f n and its corresponding P¯ are the solution of NLP. Otherwise, P¯ (P¯ P¯n) and go to step 2.

Step1:

Step2:

Solve the NLP(Pq¯ ) and obtain the solution f . Let n

n + 1, f n

f . And go to Step 1.

We can guarantee the convergence of the CG-NLP if the NLP(P¯ ) can be exactly solved during each iteration. After solving the traffic equilibrium, we acquire path flow f and the corresponding path set P , which consists of the route choice and the charging choice. 2.4. Optimal power flow In this section, a power agency is assumed to manage the power system, focusing on the power flow, the purchase/sell of electricity and the deployment of PV generators. Like the drivers and the transportation agency, the power agency is aware of the spatial and temporal charging information. This part is devoted to formulating the model for optimizing the locations and sizes of the PV generators given the electricity load and construction budget for the power agency. Let Q denote the scenario set comprising different scenarios of different electric loads and different power output levels of PV generator. Let M and L denote the bus set and link set of the power grid, respectively. Inspired by the power system’s mathematical formulation techniques proposed by Gan et al. (2015) and Zhang et al. (2018), we formulate this model as an MISOCP. MISOCP:

min

(ce+, qp0,+q

q Q

x , y , s, v, l

s. t .

(c3 x mpv + c4 ympv )

m M {0}

Smn, q = sm, q +

(Shm, q h M

0 = s0, q +

|Smn, q

lmn, q vm, q

lmn, q

|Imn |2

| Vm|2

vm, q

_

z h0 lh0, q)

q

q

q

q

Q , (m , n )

Q , (m , n)

|Vm |2

q

Q, m

M {0}&m

Q

Q, m

s0, q = p0, q

s0,load q

q

Q

smpv, q = pmpv, q + jqmpv, q pv q ym

ympv

|pmpv, q |2 + |qmpv, q |2

q

q q

ympv

(19)

q

Q , (m , n )

L

(21) (22)

L

(23)

M

(24)

M {0}

(25) (26)

Q, m

Q, m

Q, m

n

(20)

L

Q, m

q

|qmpv, q|

(18)

|z mn |2 lmn, q

smload ,q

pmpv, q

v0|

Bpv

z hm lhm, q)

sm, q = smpv, q

0

|vm, q

0

vn, q = 2Re (zmn Smn, q)

|2

m M

m

(Sh0, q h M

vm, q

ce, qp0, q ) +

(27)

M

(28)

M

(29)

M q

Q, m

(30)

M 8

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

smload ,q =

n Nm

p0, q = p0,+q

p0,+q

(LnC, q + LnI, q + LnR, q + LnEV ,q )

M

(31)

q

Q

(32)

0

q

Q

(33)

M {0}

(34)

x mpv = {0, 1} 0

Q, m

p0, q

0, p0, q

ympv

q

m

pv xmpv ymax

m

(35)

M {0}

The object function consists of two parts. The first part is the cost/profit of purchasing/selling power. and p0, q are the purchasing power and selling power, respectively, at root bus 0 in period q Q ; and ce+, q and ce, q are the corresponding unit power price2. The second part is used to control the power quality which is usually measured by voltage deviations3 (Zhang et al., 2018). vm, q and v0 are respectively the square of the nodal voltage at bus m and root bus 0. The coefficient is used to balance the two parts of the object.4 Constraint (18) is the budget constraints. x mpv is a binary variable that is equal to 1 if the PV generator is built at bus m and 0 otherwise. ympv is the corresponding PV load capacity at bus m ; c3 and c4 are respectively the fixed cost to build a PV generator and the flexible cost to add a per-unit load capacity; Bpv is the budget for the PV generator building. Constraint (19) and constraint (20) are respectively used to calculate the total power demand at bus m and root bus 0. Smn, q is the apparent power flow in distribution line (m , n) in scenario q Q , where n is closer to the root bus 0 in the distribution power network; sm, q is the total apparent power injection at bus m , and s0, q is the power consumption of the distribution power system; z hm is the impendence of line (h, m ), with z hm denoting its conjugate; lhm, q is the square of the line (h, m ) current in scenario q Q . Constraint (21) shows the voltage drop at line (m , n) . Re () is a function to denote the real part of a complex number. Constraint (22) adopts the convex relaxation proposed by Gan et al. (2015), which is the relationship between voltage and current. Constraint (23) and constraint (24) ensure that the current and voltage never exceed the line limits. Imn is the current capacity of line (m , n) ; Vm and Vm are respectively the lower and

p0,+q

_

upper limit of bus voltage. Constraint (25) and constraint (26) respectively denote the power consumption or production at bus m and Q ; p0, q is the amount of bus 0. smpv, q and smload , q are respectively the PV power injection and the electricity load at bus m in scenario q power selling/purchasing at root bus 0. Constraint (27) guarantees the use of the reactive control technology in the PV generator. pmpv, q and qmpv, q are respectively the active part and reactive part of the PV power injection. Constraint (28) ensures that the active power production at a PV generator does not exceed the maximum active power that it can produce in scenario q . q is the PV power transfer ration in scenario q Q . Constraint (29) guarantees the maximum reactive compensation power. Constraint (30) guarantees that the apparent power production is less than the PV generator capacity. Constraint (31) is used to calculate the load at bus m . LnC, q , LnI, q , and LnR, q are respectively the commercial load, industrial load, and residual load at transportation node n , and we assume that they are given in this research; LnEV , q is the EVs’ charging load at transportation node n . Constraint (32) separates the active power transaction at node 0 into two parts: purchasing power and selling power. Constraint (33) guarantees the non-negativity of the two pv parts. Constraint (34) ensures that x mpv are binary variables. Constraint (35) ensures that the maximum size of a PV generator is ymax , the maximum load capacity of a single PV generator. Note that the current, voltage, power, resistance, and reactance are all complex numbers. As a result, all the constraints should be built on the computation rule of complex numbers. Readers may refer to Gan et al. (2015) and Zhang et al. (2018) for a more theoretical explanation of the modeling of power systems. The model is an MISOCP and can be solved easily using common solvers, e.g., CPLEX. 3. Joint optimal design model 3.1. Model formulation This part focuses on formulating the joint optimization on the deployment of charging stations and PV power plants considering a time-dependent charging fee strategy to minimize the coupled transportation and power system cost. The system cost consists of drivers’ travel time on the road, queuing and charging time at charging stations, power purchasing (or selling) cost (or profit) at the substations with outside power system, and the penalty of undesired voltage deviation. Note that the money flow between different stakeholders, i.e., electricity fee and charging fee, is not considered in the total system cost. By combining the VI and MISOCP, the JOD model can be formulated as follows. JOD:

min

x , y, p, s, v, l, f , c

s. t .

q Q

p Pw

Constraints(18)

w W

f pw, q cpw, q

Epw, n, q (pn, q,0 + pn, q ) + (ce+, qp0,+q n N

ce, qp0, q ) +

m M

|vm, q

v0|

(35)

2 Since purchasing price (ce+, q ) is higher than the selling price (ce, q ), simultaneously purchasing and selling power is not allowed, as it results in a larger cost. 3 The absolute value can be easily replaced by vm, q by adding two linear inequality vm, q > vm, q v0 and vm, q > vm,q + v0 4 We assume , designed by the power supply quality requirement in real life, is given in this study.

9

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

(c1 xncs + c2 yncs )

n N

yncs

cs x ncs ymax

x ncs

n

= {0, 1} *

n

pmin

pn, q

pmax

LnEV ,q

=

q

p Pw

cpw, q =

(38)

N

p Pw

yncs

w W

(37)

(39)

N

w W

LnEV ,q

n

N,

f pw, q

w n, p

Q, n

cpw, q f pw, q

ta (va, q ga )

w p, a

q

(40)

Q q

Q, n

N

(41) (42)

N f pw, q

0

q

Q , fq

(wn (vn, q, En, q sn, yncs )

+

a A

Vq w p, n

(43)

+ Epw, n, q (sn + pn, q ))

q

Q

(44)

n N

w W

p Pw

4

va, q

ta (va, q ga ) = ta0 1 + 0.15 va, q =

(36)

N

n

yncs

Bcs

q

ga

w w p, a f p, q

q

Q, a Q, a

sn vn, q diag (En, q EnT, q)

wn (vn, q, En, q sn, yncs ) =

2vn, q En, q

A

(45)

A 6 i=0

i, yncs (sn vn, q

(46)

En, q )i , vn, q En, q > 0

0,

vn, q En, q = 0

q

Q, n

N (47)

where Vq in constraint (43) is denoted as follows: p Pq

f pw, q

f pw, q = Dw, q

0

w

w

W, p

W

(48) (49)

Pw

In the above formulation, f pw, q is the flow of path p Pw of O–D pair w W in scenario q Q with cpw, q being the corresponding path cost and Epw, n, q being the corresponding charging amount at transportation node n N ; pn, q,0 and pn, q are the electricity price and unit charging fee, respectively, at transportation node n N in scenario q Q , with pmin and pmax being the charging-congestion price limitation set by the city manager; p0,+q and p0, q are the purchasing power and selling power, respectively, at root bus 0 in scenario q Q ; ce+, q and ce, q are the corresponding unit power price, respectively; vm, q and v0 are respectively the square of nodal voltage at bus m and root bus 0; x ncs a binary variable that is equal to 1 if a charging station is built at node n and is 0 otherwise; yncs is the corresponding charging piles at node n ; Bcs is the budget for the charging station with c1 and c2 being the unit cost of charging station and charging pile, respectively. The objective function calculates the social cost. The first part is travel time, queuing time, and charging time. The second part is power transaction cost/profit. The third part is the penalty for undesired voltage deviation. Constraint (36) is the budget constraint for the charging station. Constraint (37) shows the maximum charging piles in a charging station. Constraints (38) and (39) are commonly used binary and non-negative integer constraints. Constraint (40) guarantees that the charging fee is within the limit published by the government. Constraint (41) calculates the charging demand at each node in each period. Constraint (42) ensure the charging demand at each node is less than its available charging ability which is linear to the size of charging station yncs and is the maximum charging demand one charging pile can serve. Constraint (43) is the VI that denotes the traffic equilibrium state under a given design of charging station and scenario-dependent price for each scenario. Constraints (44)–(47) are used calculate each path’s cost specifically as explained in subsection 2.3. Constraints (48)–(49) determine the feasible domain of the VI in constraint (43). 3.2. Solution procedure As formulated above, in the JOD, the functional forms of the objective function and the decision variables are relatively complicated; thus, it is difficult to apply conventional optimization techniques, e.g., descent direction algorithm and branch-and-bound algorithm, here. Furthermore, given a charging station deployment and a time-dependent pricing scheme, it could be time-consuming to evaluate the associated system cost by solving the NLP and MISOCP. Thus, our goal is to obtain a high-quality solution within an acceptable computational time. An SBO algorithm is incorporated to solve the JOD model. The SBO is widely used in design of transportation facilities, i.e., charging schemes (Rodriguez-Roman and Ritchie, 2019; Chen et al., 2014), and uses surrogate modeling techniques to look for local or global optima (Queipo et al., 2005; Han and Zhang, 2012). The surrogate modeling process builds a 10

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

surrogate model by utilizing sample data, which are collected by running time-expensive codes, i.e., the JOD model. In the above framework, the most important part is the heuristic method to adjust the solution. There are many types of surrogate models, such as the response surface method (RSM) and Kriging and radial basis function (RBF) (Han and Zhang, 2012; Chen et al., 2014). However, the RSM requires a large number of initial points (quadratic RSM needs O (n2) initial points), which means a lot of computation time. x2

Kriging takes the formulation of a Gaussian kernel function (x ) = e 2 , where it is difficult to fit the near point mutation, and usually causes a singular matrix when building the surrogate model. In this research, we employ the RBF model and take the format (x ) = x 3 . The approximate object function of the surrogate model is defined as follows: n

F (v ) =

( vi

i

i=1

v ) + P (v )

(50)

where i is the weight coefficient corresponding to the i th variable vi in the variable set; ( vi v ) is the basis function, whose variable is the Euclidean distance between vi and v , and P (v ) is the global trend function. Here, we simply make P (x ) = 0 . The approximate object function should have the following equations:

F (vi ) = yi n i=0

i

i = 1,

(51)

,n

=0

(52)

where yi is the i th result in the result set. Solving the linear equation by Eqs. (50)–(52), the RBF is as follows:

F (v ) = where ,

0

+

, and

0

T

1 (F

are defined as follows:

= [ ( vi = (1T

(54)

n

= [ ( vi v )]i

0

(53)

0 1)

vj )]ij 11) 11T

(55)

n×n

(56)

1F

Here, we take the location and capacity of the charging stations and the time-dependent charging fee rate at each charging station as our variables in the variable set, and take the total system cost as the result in the result set. Note that the location and size of the PV generator are not part of the variables, as they can be easily determined by solving the MISOCP given the charging load of the transportation systems. The major part of the SBO is the surrogate modeling, which can be easily completed by Eqs. (54)–(56). The next part is to find a new and better variable. This can be determined by solving the mathematical programming problem below:

minF (v ) = v

0

T

+

s. t . Constraints (36)

1 (F

0 1)

(40)

This mathematical programming problem is nonlinear but simple, and can be solved by common nonlinear programming solvers or some heuristic algorithms, e.g., CONOPT 3.14 or genetic algorithms. This result is taken as the plan of charging stations and price strategy for the next iteration. Based on the surrogate modeling techniques mentioned above, we proposed the SBO algorithm to solve the proposed model. The steps of our proposed SBO for solving the JOD are as follows. SBO: Step1: Step2:

Randomly generate K1 plans (charging station location, size and time-dependent charging fee rate) as the initial set H1 and solve the NLP and MISCOP M (Maximum Iteration Times) do For i = 1 1. Build the surrogate model based on Set Hi and optimize the surrogate model to get N different new plans Pi1,1, Pi1,2,

2. Slightly and randomly perturb the best K2 plans in Set H0 to get N different new plans 3. Solve the NLP and MISOCP based on Pi1, j and Pi2, j , Let Hi + 1 Step3: Step4:

Hi

{Pi1,1, Pi1,2,

, Pi1,N }

{Pi2,1, Pi2,2,

j = 1,

Pi2,1,

Pi2,2,

, Pi1,N

, Pi2, N

,N

, Pi2,N }

End For Set the charging station location and size being the best planning in HM + 1. Randomly generate K3 plans on time-dependent charging fee rate as initial set

H1'and solve the NLP and MISCOP model For i = 1

M ' (Maximum Iteration Times) do

1. Build the surrogate model based on Set Hi' and optimize the surrogate model to get a new plan N different new plans Pi3,1, Pi3,2, 2. Slightly perturb the best K2 plans in Set H0 to get N different new plans Pi4,1, Pi4,2, 3. Solve the NLP and MISOCP based on Pi3, j and Pi4, j , 4. Let Hi'+ 1

Step5:

End For

Hi'

{Pi3,1, Pi3,2,

Choose the best plan in H '

M '+ 1

, Pi3, N }

{Pi4,1, Pi4,2,

j = 1,

, Pi3,N

, Pi4, N

,N.

, Pi4, N }

with optimized PV location and capacity as the final plan.

Note we can reduce the computation time by easily employing the parallel computation framework in step 1–4. We also can take the 11

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Fig. 5. Average electricity load.

route-choice set of the first 10 plans in H0' as the initial route choice set in step 3–4, which avoids too much path enumeration and reduces the iteration times of CG NLP . 4. Numerical examples This section conducts numerical examples on some test networks to evaluate the performances of proposed methodologies and derive some managerial insights. We implement our methods using MATLAB R2018a, where the NLP is solved by CONOPT 3.14 through the TOMLAB toolbox; the MISCOP is solved by CPLEX 12.8 via YALMIP toolbox and the surrogate model is solved by genetic algorithm of MATLAB R2018a. All the examples are run on a desktop PC with Intel(R) Core(TM) i7-8700 CPU and 16 GB RAM. In all tests below, we set that the charging time for one kwh equals two mins, and one kwh supports an EV travelling for five kms. The price for purchasing/selling energy is 0.55/0.45 yuan/kwh. According to the empirical time-dependent daily electricity load and solar power generating level (PGE, 2019), we divide a day into six periods (each period denotes a scenario) as follows (see Figs. 5, 6, and Table 1). The Electric load denotes the average background electricity load (including commercial, industrial, and residential load) at each bus. Solar generating level denotes the ratio of active power generation over PV generation capacity. The weight of all periods are the same and equals 1/6, as they are equally divided. We first run our algorithm on a small Nguyen-Dupius (ND) network and further test it on a larger Sioux-Falls (SF) network. 4.1. Nguyen-Dupius network The ND and the power transmission network are showed in Fig. 7. The ND network contains 13 nodes and 19 links, whose characteristics is showed in Table 2. The coupled power transmission network is a 110KV distribution network, with bus 0 being a substation and connected to a 220KV/110KV transformer with 150MVA capacity. The characteristics of this network is shown in Table 3. The coupled relation of the two networks is showed in Table 4. To give an example, the charging demands of nodes 1, 4, and 5 are supplied by bus 1. The O-D matrix is showed in Table 5. The solution we obtain consists of the charging stations’ locations and sizes, PV generators’ locations and capacity, and the timedependent price, which yields an average social cost of 2.33 × 106 yuan/hour. Fig. 8shows the locations of deployed charging stations

Fig. 6. Average solar generating ability. 12

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Table 1 Electricity load and solar generating ability in 6 periods. Period

1

2

3

4

5

6

Time Electric load (MW) Solar generating level

22:00–2:00 1.4789 0

2:00–6:00 1.2033 0

6:00–10:00 1.6040 0.2828

10:00–14:00 1.9943 0.7984

14:00–18:00 1.9987 0.5203

18:00–22:00 1.9944 0.0038

Fig. 7. The ND transportation network (left) and the power transmission network (right). Table 2 Property of ND network. Link Start

Link End

Link Capacity

Free Travel Time (min)

Link Length (km)

1 1 4 4 5 5 6 6 7 7 8 9 9 10 11 11 12 12 13

5 12 5 9 6 9 7 10 8 11 2 10 13 11 2 3 6 8 3

2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000

40 40 40 76 40 40 40 40 40 40 40 40 76 40 40 40 40 112 40

40 40 40 80 40 40 40 40 40 40 40 40 80 40 40 40 40 120 40

Table 3 Property of power transmission network. Line Sequence

Capacity (MVA)

Length (km)

Resistance (Ohm)

Reactance (Ohm)

1 2 3 4 5 6 7

50 50 50 50 50 50 50

20 20 20 20 20 20 20

10.5 10.5 10.5 10.5 10.5 10.5 10.5

25.3 25.3 25.3 25.3 25.3 25.3 25.3

13

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Table 4 Coupled relation of ND and power transmission networks. Power.

Trans.

Power.

Trans.

Power.

Trans.

Power.

Trans.

0 1

6 1,4,5

2 3

12 7,8

4 5

10 9

6 7

2,3,11 3

Table 5 O-D matrix of ND network. Origin

1 1 4 4

Destination

2 3 2 3

Initial SOC

200 200 200 200

Maximum SOC

200 200 200 200

Safety SOC

Travel Demand

30 30 30 30

Period 1

Period 2

Period 3

Period 4

Period 5

Period 6

840 480 360 640

105 60 45 80

2100 1200 900 1600

1260 720 540 960

1470 840 630 1120

1890 1080 810 1440

Fig. 8. Charging stations and PV power plants planning result.

and PV generators. The number of charging piles in charging stations 5, 9,10, 11, and 12 are 145, 140, 145, 150, and 150, respectively. And the capacity of PV generation at bus 1, 2, and 4 are 33.06 MW, 28.37 MW, and 50.00 MW, respectively. 12 8 2, 1 5 9 13 3, 4 9 10 11 2 , and 4 9 13 3, The shortest paths for the above O-Ds are 1 respectively. We observe that all the charging stations are located in the above shortest paths, which can be explained because doing this helps reduce the social cost. With regard to the deployment of PV plants, the location plan supports our previous inference that PV plants should be deployed closely to charging stations to promote the local consumption of photovoltaic energy. Figs. 9 and 10 show the traffic flow and the power flow of the six periods, respectively. On the traffic network, the colors of the links represent the values of the corresponding ratio of road traffic volume over road capacity. On the power network, the colors represent the values of

14

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Fig. 9. Traffic equilibrium flow in six periods.

power flows, with positive values denoting the flow from bus to substation and the negative ones denoting the opposite direction. It can be observed that PV generators deployed at buses 1, 2 and 4 reach peak solar power generation levels at periods 4–5, consistent with the temporal distribution of solar power production level. The time-dependent price is showed in Table 6 (For reader’s convenience, the result is the sum of electricity price and unit charging fee). To analyze the effect of time-dependent price, we compare the system’s queuing times at charging stations under the optimal pricing strategy (referred to as Case I) with the ones without pricing strategy (referred to as Case II) in Table 7. The results reveal the significant effect of time-dependent pricing strategy in reducing the systems’ queuing times. Figs. 11–12 further show the queuing time at each charging station in Cases I and II, respectively. It can be observed that the optimal pricing strategy attempts to achieve a more efficient utilization of charging resource and reduce charging queuing time through balancing charging demands across different charging stations. For instance, without pricing, Station 12 corresponds to the highest queuing time; the optimal electricity pricing strategy sets a high charging price at it to induce drivers to charge elsewhere. Fig. 13 shows the energy transaction (a negative value represents selling power and a positive value denotes purchasing power at substations), the energy produced by PV plants and the total electricity load in six periods. We observe that the distribution power network can sell energy in periods 4 and 5, mainly because of the high solar power production level in periods 4 and 5. 4.2. Sioux-Falls network We solve the JOD model on the Sioux-Falls network. The CPU time is approximately 19 h in the computing environment mentioned above. Fig. 14 shows the optimal plan of charging stations and PV plants. The number of charging piles in charging stations 8, 10 and 14 are 120, 135, and 120, respectively. The PV generation capacity of buses 1, 3, 5, and 7 are respectively 50.00 MW, 50.00 MW, 41.43 MW, and 50.00 MW, respectively. Table 8 shows the optimal time-dependent price of different charging stations. Fig. 15 depicts the convergence curve of our surrogate-based optimization (SBO) algorithm. To evaluate the performance of SBO, we propose a valid lower bound (see Appendix E) and adopt the genetic algorithm (GA) (Yin, 2000) as a comparison. The GA has the same initial solution and the same population size (In this research, the population size is 24) in each iteration as the SBO. The result is shown in Fig. 15. The result shows that the SBO outperforms the GA and can give a proper solution after a certain number of iterations. We further conduct the sensitivity analysis on the budgets for charging stations and budget for PV power plants (see Fig. 16). It can be observed from the figure that the system cost decreases approximate linearly with the increase of PV power plants budget. However, the system cost decreases negatively exponentially with the increase of charging stations budget. The possible reason is that the energy generated by PV power plants is linear to the its capacity. And the queuing time in the charging station is nearly reciprocal to the number of servers (Kimura, 1996). Besides, considering the budget for charging stations and PV power plants

15

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Fig. 10. Power flow state in six periods. Table 6 Time-dependent price (yuan/kwh) of charging stations.

Charging Charging Charging Charging Charging

Station Station Station Station Station

5 9 10 11 12

Period 1

Period 2

Period 3

Period 4

Period 5

Period 6

0.50 0.76 0.75 0.50 1.00

0.50 0.50 0.55 0.50 0.59

0.50 0.50 0.55 0.50 1.00

0.50 0.60 0.50 0.50 1.00

0.52 0.50 0.50 0.50 1.00

0.60 0.50 0.60 0.50 0.98

Period 1

Period 2

Period 3

Period 4

Period 5

Period 6

2320 1269.81 1851.82

290 0.06 0.06

5800 85078.32 86120.66

3480 3130.77 7799.53

4060 7384.34 11580.91

5220 40850.94 42748.22

Table 7 Total queuing delay (min) of Cases I and II.

Total traffic demand Queue Delay (Case I) Queue Delay (Case II)

16

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Fig. 11. Queuing time of Case I.

Fig. 12. Queuing time of Case II.

Fig. 13. Energy transaction, solar energy and total electric load in six periods.

are CNY4.5 × 107 and CNY1.5 × 109 , respectively, we can expect that transferring the budget of PV power plants to the budget of charging stations reduces the system cost. 5. Conclusions This paper is primarily concerned about the joint optimization on the planning of charging stations and PV plants with timedependent charging fee to better manage coupled transportation and power networks. Specially, we propose an extend label-setting

17

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Fig. 14. Planning results of charging stations and PV plants. Table 8 Time-dependent price (CNY: yuan/kwh) of charging stations.

Charging Station 8 Charging Station 10 Charging Station 14

Period 1

Period 2

Period 3

Period 4

Period 5

Period 6

1.00 2.00 1.00

1.01 2.00 1.00

1.17 1.00 1.00

2.00 1.89 1.00

2.00 1.00 1.10

2.00 1.00 1.00

algorithm to solve the EVs’ JRCP that considers loop movement cases as well as EVs’ recharging amount choices at different stations. The queuing at charging stations is modeled as a M G k system. Then, a VI problem is formulated to model the equilibrium of EV traffic on transportation networks, and an optimal power flow model is proposed to model the power network flow with PV plants and optimally serve the EV charging requirements. Based on the above models, we formulate the joint optimization problem as a bilevel programming model and adopt a SBO framework to solve it. Numerical studies demonstrate the effectiveness of the optimal design of charging stations, PV plants and time-dependent charging fee strategy. Though it captures the time-variance of traffic demand and solar energy generating levels by considering a pool of periods, this research assumes that traffic flow reaches an equilibrium in each period, i.e., the routing choice and charging choice must be made in only one period. As a result, this research ignores that drivers may make choices across more than one period, i.e., changing the scheduled departure time due to congestion on roads or in charging stations, and that would result in a multi-period equilibrium model. Our current modeling cannot capture such behavioral features and we will incorporate this consideration in future study, which will also lead to more advanced strategies to better manage travelers’ choices of starting times for charging. Another direction in future study is to improve the solution efficiency and enhance the proposed framework’s potential for solving more realistic problems.

18

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Fig. 15. Convergence curve and lower bound.

Fig. 16. Sensitivity analysis on construction budgets.

CRediT authorship contribution statement Zhixiong Luo: Methodology, Software, Validation, Writing - original draft, Visualization. Fang He: Conceptualization, Writing review & editing, Funding acquisition. Xi Lin: Methodology, Writing - review & editing. Jianjun Wu: Writing - review & editing. Meng Li: Conceptualization, Supervision, Writing - review & editing, Funding acquisition. Acknowledgements This research is partially supported by grants from National Natural Science Foundation of China (71871126, 51622807, U1766205, 71931003). The authors would like to thank Prof. Zechun Hu for advice on power networks.

19

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Appendix A. Notation and definition See Table A1. Table A1 Notations and definitions. Notations

Definitions

N A W M L Nm Q

Set of nodes in transportation network Set of links in transportation network Set of O-D pairs Set of distribution buses in power network Set of distribution lines in power network Set of nodes in the transportation network whose energy is supplied by busm Sets of all scenarios Weight of scenarioq Q

q

lc Np Nc Tc Ec Wc Sc tc ec Einitial Emax Esafe

A label at node C Previous node in a label

Current node in a label (namely node c) Total cost when reaching last charging station (with positive charging amount) before getting node C SOC when reaching last charging station (with positive charging amount) before getting node C Queuing cost at last charging station (with positive charging amount) Unit charging cost of last charging station before getting node C Total travel cost from last charging station to node C Total electric consumption from last charging station to node C Initial SOC of an electric vehicle (EV) Maximum SOC of an EV Safety SOC of an EV due to range anxiety

Pw f pw

Set of path for O-D pairw Flow of a path p Pw

Dw cpw

Travel Demand of O-D pairw Total cost of a path p Pw

w cmin va ga

Minimum total cost of all path p

Pw of O-D pairw

Indicators that equals 1 if node n

N is utilized by path p

w p, a w p, n

W W W

Flow of transportation linka A Capacity of transportation linka A Indicators that equals 1 if link a A is utilized by path p

Pw, w

wn Epw, n

Queuing time of noden N Charging amount of path p Pw, w

vn En sn pn

Column Vector denoting the flows at noden N Column Vector denoting the charging demand at node n N Charging time for unit charging amount at charging station n Charging price for unit charging amount at charging station n

w

Pw Pw¯ P¯

W at the node n

Price of purchasing power at the substation in scenario q

ce, q

Price of selling power at substation in scenario q

p0,+q

Power purchasing in scenario q

p0, q

Power selling in scenario q

x mpv ympv

Bpv c3 c4 Smn,q sm,q z hm lhm, q

W , otherwise, 0. W , otherwise, 0.

N

N N

Q

Q

Q

Q

The square of nodal voltage at bus m The square of nodal voltage at root bus 0, set to be (110KV)2 Binary variable denoting whether building PV plants at bus m PV load capacity at bus m

M

M

Budgets of PV plants

Fixed cost to build a PV generation Flexible cost of per-unit load capacity Apparent power flow in distribution line(m , n) Total apparent power injection at busm The impendence of distribution line(h , m) Square of the line (h , m) ’s current

Imn Vm /Vm

Current capacity of line(m , n) Lower / Upper limit of nodal voltage at busm

smpv, q

PV’s power injection at busm

_

Pw, w

w Dual variable of constraints (15), denoting the cmin Set of path for O-D pair w W Restricted path set for O-D pair w W Restricted path set for all O-D pairs

ce+, q

vm, q v0

M

(continued on next page) 20

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Table A1 (continued) Notations

Definitions

load sm ,q

Electricity load at busm

pmpv, q /qmpv, q

Active / Reactive part of the PV power injection at busm Power transfer ratio in scenario q

q

Q

LnC, q , LnI,q , LnR, q

Ratio of maximum reactive power to the load capacity Commercial, industrial, and Residual load at transportation node n

LnEV ,q

EVs’ charging load at transportation node n Maximum load capacity of the single PV generation

pv ymax

Binary variable indicating whether building charging station at node n Integer variable denoting the number of charging piles at node n

x ncs yncs

Budget of charging stations Unit charging fee limit at charging station

Bcs pmin / pmax

Appendix B. Proof of propositions Proposition 1.. The set of optimal solutions to the JRCP contain one satisfying HMLM. Proof.. The JRCP consists of route choices, charging station choices and charging amount choices. Given any fixed route, assume the optimal charging station choice is in a sequence of c1, c2, , cn , and the charging amount and unit charging price at ci are respectively yi and Si . The electricity consumption from ci to ci + 1 is ei, i + 1. Let 0 and n + 1 respectively denote the origin node and the destination node (see Fig. B1). Then, the optimal charging amount choice at each charging station can be decided by the following linear program (LP). LP n

min

i=1

s. t .

Si yi

(B-1) k

Einitial +

y i=1 i

Einitial +

k y i=1 i

yi

i = 1,

0

k 1 i=0

k i=0

ei, i + 1

ei, i + 1 Esafe

Emax

k = 1,

k = 1,

,n

(B-2)

,n

(B-3) (B-4)

,n

The objective function (B-1) is the total charging cost, as the path travel cost and queue cost at charging station are both determined. Constraint (B-2) ensures that the SOC after leaving each charging station is less than or equal the maximum SOC; and constraint (B-3) ensures that the SOC when getting each charging station is larger or equal than the safety SOC. Constraint (B-4) states k that the charging amount at each charging station is nonnegative. Then, let z k = i = 1 yi , k {1, 2, , n}andz 0 = 0 , thus yi = z i z i 1, i {1, 2, , n} . The above LP can be reformulated as follows: RLP: n 1

min

i=1

s. t .

z i (Si

k 1

Einitial + z k

zi

1

i=0

k

Einitial + z k

zi

Si + 1) + zn Sn

ei, i + 1

i=0

0

i = 1,

ei, i + 1 Esafe

min

n 1 i=k

z k (Sk

Emax

k = 1,

k = 1,

n

(B-6)

,n

(B-7) (B-8)

,n

For notation brevity, let z = Esafe + blem can be rewritten as: FRLP:

(B-5)

_ k

k e i = 0 i, i + 1

Einitial and z¯k = Emax +

k 1 e i = 0 i, i + 1

Sk + 1) + zn Sn

Einitial for any k = 1,

, n . Thus, the pro-

(B-9)

Fig. B1. The illustrated route. 21

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

s. t .

zk zk To

z

zk

z¯k

z

k = 1,

_ k

zk

0

1

(B-10)

n

n

(B-11)

k = 1,

(B-12)

,n

the problem feasible, we must have Emax ek, k + 1 Esafe, k {1, 2, , n} . We next show that {1, 2, , n 1} must hold. As Sk > 0 , we claim that at the optimal solution to FRLP without constraint (B-12) is given

make z¯k , k

_ k+1

k = 1,

by:

z k = z¯k k = 1, z k = z k = 1,

,n ,n

_ k

1; Sk < Sk + 1 1; Sk Sk + 1

zn = z

_ n

If z

z¯k ,

_ k+1

k

{1, 2,

1} , then above solution is also the optimal solution to the original FRLP; otherwise, the above

,n

solution is outside the feasible domain of FRLP, and that implies the optimal solution to FRLP must be bounded by constraint (B-12), i.e., for at least one k {1, 2, , n 1} we have z k z k 1 = 0 = yk , suggesting that the vehicle does not charge at the k th charging station; this contradicts our setting that c1, c2, , cn is the optimal charging station choice for this route. Then, we prove that z z¯k , k {1, 2, , n 1} must hold, and the optimal solution is given above. Thus, for any route, at least one of the optimal

_ k+1

charging policies satisfies HMLM, which further implies that one of the optimal solutions to JRCP satisfies HMLM. Proof completes■ Proposition

Sc'

Sc , ec'

lc' = {N p' , Nc , Tc' , Ec' , Wc' , Sc' , tc' , ec' } dominates

2..

Ec'

Ec and

ec

ec .

ec'

lc = {Np, Nc , Tc, Ec , Wc, Sc, tc, ec }

if

Tc' + Wc' + tc'

Tc + Wc + tc ,

Proof:. We observe that we can prove the proposition if: a routing and charging solution represented by label lc reaches Nc with travel cost TNc and remaining SOC E Nc , then there always exists a routing and charging solution represented by label lc' with reaches Nc with travel cost TN' c TNc and remaining SOC E N' c E Nc ; this claim is straightforward since no matter how to make routing and charging decisions after lc , the vehicle with label lc' can make the same decision following that of lc to ensure that the travel cost at the destination is less or equal. Note that TNc = Tc + Wc + tc + Sc yc , E Nc = Ec + yc ec , TN' c = Tc' + Wc' + tc' + Sc' yc' , E N' c = Ec' + yc' ec' , where yc and yc' are charging amount at the previous charging stations of lc and lc' respectively; for any yc , we let yc' = min(yc , Emax Ec' ) yc , then we have:

TN' c = Tc' + Wc' + tc' + Sc' yc'

Tc + Wc + tc + Sc yc = TNc

E N' c = Ec' + yc'

ec' = Ec' + yc

E N' c = Ec' + yc'

ec' = Emax

ec' ec'

Ec + yc Emax

ec

ec = E Nc (ifyc

Emax

E Nc (ifyc > Emax

Ec' )

Ec' )

Thus, we above conditions, we can ensure that lc' dominates lc . Proof completes. ■ Proposition 3.. The VI proposed in (5) is equivalent to the equilibrium condition (3)-(4). Proof.. The Eqs. (1)–(5) are listed as follows.

V: f pw

f pw = Dw

p Pw

0

w

W, p

TEC:

cpw

w cmin

f pw (cpw

w cmin )

=0

VI:

w

w W

(1) (2)

Pw

0

w w

p Pw

W

W, p

W, p

cpw f pw

(3)

Pw

(4)

Pw

f pw

0

f

V

(5)

Suppose that there exists a solution that satisfy the VI constraints. And the path flow for OD pair w W are {f pq , f pq , , f pq , , f pq } , where n w is the number of path for OD pair w W . Consider the constraints (1) and constraints (2), there

f pw

1

2

j

j

i

nw

must exists a positive variable, i.e., f pw > 0 . We can generate another feasible path flow f pw for OD pair w W such that f pw = 0 , i i f pw = f pw + f pw and keep the rest path flow unchanged. For other O-D pairs, we make no change. Taking this new feasible path flow i

in to VI constraint (5), the result is

(c pwj

c pwi ) f pq

i

0,

i, j = 1, ..,n w , w

(B-13)

W

> 0 , thus Note that must holds, which indicates that, for an OD pair w W , if the path flow is positive, then the w = minc pwi , the above inequality yields that path cost must be the minimum cost among all the path cost. Taking cmin f pw i

c pwj

c pwi

i

22

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

c pwi

w cmin

(c pwi

0

i = 1, ..,n w , w

w cmin ) f pw = 0 i

(B-14)

W

i = 1, ..,n w , w

(B-15)

W

which is exactly the traffic equilibrium condition. On the other hand, the traffic equilibrium condition yields that when c pwi > min (c pwi ), f pw = 0 and when f pw = 0 , c pwi min (c pwi ) . Considering the feasible region V , then p Pq

i = 1,

, nw

cpw

f pw

i

f pw

i

i = 1,

, nw

0 can be easily verified. Thus the VI constraints (5) hold. Proof completes. ■

Appendix C. Evaluation of extended label-setting algorithm C.1 Solution quality We first compare the solution result of those two methods, i.e., the mix-integer linear programming (MILP) model proposed by He et al. (2014) and the extended label-setting algorithm on the Sioux Falls network (see Fig. C1). The travel time (before comma, unit: minute) and equivalent energy consumption (after comma, unit: km) are marked near the link. We test three O-D pairs 1–13, 5–24, and 7–20, with equivalent initial SOC being 55 km, 45 km, and 35 km, respectively. The equivalent maximum and minimum SOC are

Fig. C1. Sioux Falls network (left) and result comparison (right).

23

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

150 km and 25 km, respectively. The queuing time at station 11 and 16 are 2 min and 14 min respectively. The charging speed is 2 km amount of energy per minutes. Integer programming model is solved by CPLEX 12.8 through Yalmip tool box. All the computation is conducted in MATLAB R2018a. The results are shown in Fig. C1, where the two methods give the same results for the first two O-D pairs. For O-D pair 7–20, MILP model gives the path choice being 7–16-17–19-20 with charging amount being 13.4 km at a total cost of 34.96 min. On the other hand, extended label-setting algorithm gives the path choice being 7-18-16-18-20 with charging amount being 11.6 km at a total cost of 33.04 min. The comparison result indicates that extended label-setting algorithm can capture the charging loop behavior and incur a lower cost for EV drivers. C.2 Solution time We further test the 76 O-D pairs selected from Stabler (2019) based on the Sioux Falls network (see Fig. C1) with same parameter setting (link travel time, energy consumption, maximum and minimum SOC), except that all EVs’ equivalent initial SOC is set to be 35 km. The result shown in Fig. C2 indicates that the extended label setting algorithm can efficiently solve the joint routing and charging problem for EV drivers.

Fig. C2. Boxplot of the solutions time.

Appendix D. Parameters of Sioux-Falls network tests See Tables D1–D4.

24

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Table D1 Property of SF network. Link Start

Link End

Link Capacity

Free Travel Time (min)

Link Length (km)

1 1 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 8 8 8 8 9 9 9 10 10 10 10 10 11 11 11 11 12 12 12 13 13 14 14 14 15 15 15 15 16 16 16 16 17 17 17 18 18 18 19 19 19 20 20 20 20 21 21 21 22 22

2 3 1 6 1 4 12 3 5 11 4 6 9 2 5 8 8 18 6 7 9 16 5 8 10 9 11 15 16 17 4 10 12 14 3 11 13 12 24 11 15 23 10 14 19 22 8 10 17 18 10 16 19 7 16 20 15 17 20 18 19 21 22 20 22 24 15 20

3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000

54 36 54 45 36 36 36 36 18 54 18 36 45 45 36 18 27 18 18 27 90 45 45 90 27 27 45 54 36 72 54 45 54 36 36 54 27 27 36 36 45 36 54 45 27 27 45 36 18 27 72 18 18 18 27 36 27 18 36 36 36 54 45 54 18 27 27 45

54 36 54 45 36 36 36 36 18 54 18 36 45 45 36 18 27 18 18 27 90 45 45 90 27 27 45 54 36 72 54 45 54 36 36 54 27 27 36 36 45 36 54 45 27 27 45 36 18 27 72 18 18 18 27 36 27 18 36 36 36 54 45 54 18 27 27 45

(continued on next page) 25

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Table D1 (continued) Link Start

Link End

Link Capacity

Free Travel Time (min)

Link Length (km)

22 22 23 23 23 24 24 24

21 23 14 22 24 13 21 23

3000 3000 3000 3000 3000 3000 3000 3000

18 36 36 36 18 36 27 18

18 36 36 36 18 36 27 18

Table D2 Property of power transmission network. Line Sequence

Capacity(MVA)

Length(Km)

Resistance(Ohm)

Reactance(Ohm)

1 2 3 4 5 6 7

50 50 50 50 50 50 50

40 40 40 40 40 40 40

21 21 21 21 21 21 21

50.6 50.6 50.6 50.6 50.6 50.6 50.6

Table D3 The coupling of this two network. Power.

Trans.

Power.

Trans.

Power.

Trans.

Power.

Trans.

0 1 2 3

10 9 16,17 11,14

4 5 6 7

4,5 2,6,8 1,3 15,19

8 9 10 11

20 21,22 7,18 12

12 13

23,24 13

Table D4 O-D Matrix. Origin

1 2 3 4 4 6 7 12 13 18 19 19 20 23

Destination

20 13 19 18 19 23 12 7 2 4 3 4 1 6

Initial SOC

180 180 180 180 180 180 180 180 180 180 180 180 180 180

Maximum SOC

180 180 180 180 180 180 180 180 180 180 180 180 180 180

Safety SOC

Traffic Demand

20 20 20 20 20 20 20 20 20 20 20 20 20 20

Period 1

Period 2

Period 3

Period 4

Period 5

Period 6

320 160 320 400 240 400 480 440 80 320 160 160 200 400

40 20 40 50 30 50 60 55 10 40 20 20 25 50

800 400 800 1000 600 1000 1200 1100 200 800 400 400 500 1000

480 240 480 600 360 600 720 660 120 480 240 240 300 600

560 280 560 700 420 700 840 770 140 560 280 280 350 700

720 360 720 900 540 900 1080 990 180 720 360 360 450 900

Appendix E. Lower bound of JOD The methodology for determining a valid lower bound of social costs is presented, defined by the objective function in JOD, given the budget of charging stations and PV generations. The lower bound is obtained by ignoring the queuing time at charging stations and getting the lower bound from two parts: travel time and charging time plus the power transaction cost and penalty of voltage deviations. The first part calculates the system optimal state of the transportation system in LB1 and the second part decides the optimal power flow state by optimally assigning the charging demand of EVs in LB2. The validate lower bound of JOD then can be determined by B = LB1 + LB2 . LB1: 26

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

LB1 = min

q Q

v, E

vaw, q = Dqw

s. t . va, q =

a A

0

q

vaw, q Ea q

va, q ta (va, q) + s

q

vaw, q

w W

Eqw Eqw

a A

Q, w Q, a

dqw (E0,wq

Q, w

w W

Eqw (E-1)

W

(E-2)

A Esafe )

q

Q, w

(E-3)

W

(E-4)

W

In LB1, the object function consists of travel time and charging time where va, q is the traffic flow on link a A in period q Q and ta (va, q) is the corresponding travel time. Eqw is the total charging amount of O-D pair w W in period q Q and s is charging time for unit charging amount. Constraint (E-1) guarantees the flow-balance where is the node-link adjacency matrix; vaw, q is the vector denoting the flow of O-D pair w W on link a A in period q Q ; Dqw is the vector with the origin node being dqw (the demand of O-D pair w W in period q Q ), the destination node being dqw , and the other node being 0. Constraint (E-2) calculates the flow on link a A in period q Q by aggregating the traffic flow of each O-D pair. Constraint (E-3) and Constraint (E4) make sure each O-D pair’s charging amount is larger than zero and the difference between total electricity consumption and total initial usable electricity where Ea is the electricity consumption on link a A ; E0,wq is the initial electricity of O-D pair w W ; and Esafe is the safe electricity due range anxiety. LB2:

LB2 =

min

q Q

x , y, p, s, v, l, f

(ce+, qp0,+q

ce, qp0, q ) +

m M

|vm, q

v0|

s. t .Constraints (18)-(39) p Pw

f pw, q LnEV ,q = LnEV ,q

f pw, q = d w, q

0

q

q

Q, w

w W

p Pw

yncs

q

Q, w

W, p f pw, q Enw, p, q

Q, n

W

(E-5) (E-6)

Pw q

Q, n

N

(E-7) (E-8)

N

In LB2, the object function consists of power transaction cost and und undesired voltage deviation. The first part includes the cost of purchasing and the profit of selling power. p0,+q and p0, q are the purchasing power and selling power, respectively, at root bus 0 in period q Q ; and ce+, q and ce, q are the corresponding unit power price. The second part is used to control the power quality measured by voltage deviations (Zhang et al., 2018). Constraints (18)-(39) are the constraints of power system (see Section 2.4) and the constraints on charging station deployment (see Section 3.1). Constraint (E-5) ensures the flow balance of each O-D pairs where f pw, q is the flow of routing-charging choice p Pw for O-D pair p Pw in period q Q and d w, q is the corresponding traffic demand. Constraint (E-6) ensures the non-negativity of flow f pw, q . Constraint (E-7) calculates the charging demand at each charging station where Enw, p, q is the charging demand at charging station n N of routing-charging choice p Pw . Constraint (E-8) ensures that the charging demand at each charging station is less than its available charging ability, which is linear to the size of charging station yncs and is the available charging demand one charging pile can serve. Since the sum of object function in LB1 and LB2 is lower than JOD by neglecting the waiting time and the feasible domain is enlarged by neglecting the traffic equilibrium constraints and separate optimization on LB1 and LB2, the above result LB = LB1 + LB2 determines a valid lower bound of JOD. Appendix F. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.trd.2020.102247.

References Aghassi, M., Bertsimas, D., Perakis, G., 2006. Solving asymmetric variational inequalities via convex optimization. Oper. Res. Lett. 34 (5), 481–490. Ahuja, R.K., Magnanti, T.L., Orlin, J.B., 1993. Network Flows. Prentice hall, Upper Saddle River, NJ. Andrews, M., Dogru, M.K., Hobby, J.D., Jin, Y., Tucci, G.H., 2013. Modeling and optimization for electric vehicle charging infrastructure. In: In IEEE Innovative Smart Grid Technologies Conference, pp. 1–10. Chen, X., Zhang, L., He, X., Xiong, C., Li, Z., 2014. Surrogate-based optimization of expensive-to-evaluate objective for optimal highway toll charges in transportation network. Comput.-Aided Civ. Infrastruct. Eng. 29 (5), 359–381. Denholm, P., Kuss, M., Margolis, R.M., 2013. Co-benefits of large scale plug-in hybrid electric vehicle and solar PV deployment. J. Power Sources 236, 350–356. Dommel, H.W., Tinney, W.F., 1968. Optimal power flow solutions. IEEE Trans. Power Appar. Syst. 10, 1866–1876. Dong, J., Liu, C., Lin, Z., 2014. Charging infrastructure planning for promoting battery electric vehicles: an activity-based approach using multiday travel data. Transport. Res. Part C: Emerg. Technol. 38, 44–55.

27

Transportation Research Part D 79 (2020) 102247

Z. Luo, et al.

Elgowainy, A., Burnham, A., Wang, M., Molburg, J., Rousseau, A., 2009. Well-to-wheels energy use and greenhouse gas emissions of plug-in hybrid electric vehicles. SAE Int. J. Fuels Lubr. 2 (1), 627–644. Figueiredo, R., Nunes, P., Brito, M.C., 2017. The feasibility of solar parking lots for electric vehicles. Energy 140, 1182–1197. Gan, L., Li, N., Topcu, U., Low, S.H., 2015. Exact convex relaxation of optimal power flow in radial networks. IEEE Trans. Autom. Control 60 (1), 72–87. Han, Z.H., Zhang, K.S., 2012. Surrogate-based optimization. In Real-world applications of genetic algorithms. InTech. Harker, P.T., Pang, J.S., 1990. Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications. Math. Program. 48 (1–3), 161–220. He, F., Wu, D., Yin, Y., Guan, Y., 2013. Optimal deployment of public charging stations for plug-in hybrid electric vehicles. Transport. Res. Part B: Methodol. 47, 87–101. He, F., Yin, Y., Lawphongpanich, S., 2014. Network equilibrium models with battery electric vehicles. Transport. Res. Part B: Methodol. 67, 306–319. He, F., Yin, Y., Zhou, J., 2015. Deploying public charging stations for electric vehicles on urban road networks. Transport. Res. Part C: Emerg. Technol. 60, 227–240. He, F., Yin, Y., Wang, J., Yang, Y., 2016. Sustainability SI: optimal prices of electricity at public charging stations for plug-in electric vehicles. Networks Spatial Econ. 16 (1), 131–154. Hodgson, M., Berman, O., 1997. A billboard location model. Geogr. Environ. Model. 1, 25–45. Hõimoja, H., Rufer, A., Dziechciaruk, G., Vezzini, A., 2012. June). An ultrafast EV charging station demonstrator. In: International Symposium on Power Electronics Power Electronics, Electrical Drives, Automation and Motion. IEEE, pp. 1390–1395. Ip, A., Fong, S., Liu, E., 2010. Optimization for allocating BEV recharging stations in urban areas by using hierarchical clustering. In: Advanced Information Management and Service (IMS), 2010 6th International Conference on IEEE, pp. 460–465. Kang, J.E., Recker, W., 2014. Strategic hydrogen refueling station locations with scheduling and routing considerations of individual vehicles. Transport. Sci. 49 (4), 767–783. Kimura, T., 1996. A transform-free approximation for the finite capacity M/G/s queue. Oper. Res. 44 (6), 984–988. Lawphongpanich, S., Hearn, D.W., 1984. Simplical decomposition of the asymmetric traffic assignment problem. Transport. Res. Part B: Methodol. 18 (2), 123–133. Li, M., Jia, Y., Shen, Z., He, F., 2017. Improving the electrification rate of the vehicle miles traveled in Beijing: a data-driven approach. Transport. Res. Part A: Policy Pract. 97, 106–120. Mellit, A., Pavan, A.M., 2010. A 24-h forecast of solar irradiance using artificial neural network: application for performance prediction of a grid-connected PV plant at Trieste, Italy. Sol. Energy 84 (5), 807–821. Momtazpour, M., Butler, P., Ramakrishnan, N., Hossain, M.S., Bozchalui, M.C., Sharma, R., 2014. Charging and storage infrastructure design for electric vehicles. ACM Trans. Intell. Syst. Technol. (TIST) 5 (3), 42. Moon, H., Park, S.Y., Jeong, C., Lee, J., 2018. Forecasting electricity demand of electric vehicles by analyzing consumers’ charging patterns. Transport. Res. Part D: Transp. Environ. 62, 64–79. Munkhammar, J., Bishop, J.D., Sarralde, J.J., Tian, W., Choudhary, R., 2015. Household electricity use, electric vehicle home-charging and distributed photovoltaic power production in the city of Westminster. Energy Build. 86, 439–448. Pena-Bello, A., Barbour, E., Gonzalez, M.C., Patel, M.K., Parra, D., 2019. Optimized PV-coupled battery systems for combining applications: Impact of battery technology and geography. Renew. Sustain. Energy Rev. 112, 978–990. PG&E, 2019. < https://www.pge.com/nots/rates/2000_static.shtml > (accessed 09.07.2019). Pourazarm, S., Cassandras, C.G., Wang, T., 2016. Optimal routing and charging of energy-limited vehicles in traffic networks. Int. J. Robust Nonlinear Control 26 (6), 1325–1350. Queipo, N.V., Haftka, R.T., Shyy, W., Goel, T., Vaidyanathan, R., Tucker, P.K., 2005. Surrogate-based analysis and optimization. Prog. Aerosp. Sci. 41 (1), 1–28. Recker, W.W., Kang, J.E., 2010. An Activity-Based Assessment of the Potential Impacts of Plug-In Hybrid Electric Vehicles on Energy and Emissions Using One-Day Travel Data. Rodriguez-Roman, D., Ritchie, S.G., 2019. Surrogate-based optimization for the design of area charging schemes under environmental constraints. Transport. Res. Part D: Transp. Environ. 72, 162–186. Salpakari, J., Rasku, T., Lindgren, J., Lund, P.D., 2017. Flexibility of electric vehicles and space heating in net zero energy houses: an optimal control model with thermal dynamics and battery degradation. Appl. Energy 190, 800–812. Stabler B. Transportation Networks. < https://github.com/bstabler/TransportationNetworks > (accessed in 11.11.2019). State Council of the People’s Republic of China, 2015. < http://english.gov.cn/news/top_news/2015/10/13/content_281475210562536.htm > (accessed 02.06. 2019). Thomas, D., Deblecker, O., Ioakimidis, C.S., 2018. Optimal operation of an energy management system for a grid-connected smart building considering photovoltaics’ uncertainty and stochastic electric vehicles’ driving schedule. Appl. Energy 210, 1188–1206. Tulpule, P.J., Marano, V., Yurkovich, S., Rizzoni, G., 2013. Economic and environmental impacts of a PV powered workplace parking garage charging station. Appl. Energy 108, 323–332. Wang, C., He, F., Lin, X., Shen, Z.J.M., Li, M., 2019. Designing locations and capacities for charging stations to support intercity travel of electric vehicles: an expanded network approach. Transport. Res. Part C: Emerg. Technol. 102, 210–232. Wu, F., Sioshansi, R., 2017. A stochastic flow-capturing model to optimize the location of fast-charging stations with uncertain electric vehicle flows. Transport. Res. Part D: Transp. Environ. 53, 354–376. Xinhua Net, 2018. < http://www.xinhuanet.com/power/2018-12/25/c_1123901532.htm > (accessed 02.06.2019). Xu, M., Meng, Q., Liu, K., 2017. Network user equilibrium problems for the mixed battery electric vehicles and gasoline vehicles subject to battery swapping stations and road grade constraints. Transport. Res. Part B: Methodol. 99, 138–166. Xydas, E., Marmaras, C., Cipcigan, L.M., Jenkins, N., Carroll, S., Barker, M., 2016. A data-driven approach for characterising the charging demand of electric vehicles: a UK case study. Appl. Energy 162, 763–771. Yin, Y., 2000. Genetic-algorithms-based approach for bilevel programming models. J. Transp. Eng. 126 (2), 115–120. Yu, A., Wei, Y., Chen, W., Peng, N., Peng, L., 2018. Life cycle environmental impacts and carbon emissions: a case study of electric and gasoline vehicles in China. Transport. Res. Part D: Transp. Environ. 65, 409–420. Zhang, H., Moura, S.J., Hu, Z., Qi, W., Song, Y., 2017. A second-order cone programming model for planning PEV fast-charging stations. IEEE Trans. Power Syst. 33 (3), 2763–2777. Zhang, H., Moura, S.J., Hu, Z., Qi, W., Song, Y., 2018. Joint PEV charging network and distributed PV generation planning based on accelerated generalized benders decomposition. IEEE Trans. Transp. Electrif. 4 (3), 789–803. Zhang, W., Ge, W., Huang, M., Jiang, J., 2015. Optimal day-time charging strategies for electric vehicles considering photovoltaic power system and distribution grid constraints. Math. Probl. Eng., 2015.

28