Optimal Location of Fast Charging Stations for Mixed Traffic of Electric Vehicles and Gasoline Vehicles Subject to Elastic Demands

Abstract

With the rapid development of electric vehicles (EVs), one of the urgent issues is how to deploy limited charging facilities to provide services for as many EVs as possible. This paper proposes a bilevel model to depict the interaction between traffic flow distribution and the location of charging stations (CSs) in the EVs and gasoline vehicles (GVs) hybrid network. The upper level model is a maximum flow-covering model where the CSs are deployed on links with higher demands. The lower level model is a stochastic user equilibrium model under elastic demands (SUE-ED) that considers both demands uncertainty and perceived path constraints, which have a significant influence on the distribution of link flow. Besides the path travel cost, the utility of charging facilities, charging speed, and waiting time at CSs due to space capacity restraint are also considered for the EVs when making a path assignment in the lower level model. A mixed-integer nonlinear program is constructed, and the equivalence of SUE-ED is proven, where a heuristic algorithm is used to solve the model. Finally, the network trial and sensitivity analysis are carried out to illustrate the feasibility and effectiveness of the proposed model.

1. Introduction

The market share of new energy vehicles without fossil fuels has been increasing rapidly in recent years, especially electric vehicles (EVs), which provide better performance, higher efficiency, and zero emissions [1,2]. To promote the usage of EVs, governments around the world have successively issued a series of active policies to provide subsidies for EV purchases and to deploy public charging infrastructures at convenient locations [2,3]. However, the market acceptance of EVs seems to fluctuate with the decreasing subsidy. On the one hand, limited driving range, insufficient public charging infrastructures, and longer charging time are the main reasons for at-home charging. On the other hand, the deployment of EV charging stations may induce more trips [4,5], that is, demands between each origin to destination (OD) pair that respond to the available fast charging are elastic, so does the route choice. Therefore, it is especially important to study how to effectively and reasonably deploy EV public charging stations within the hybrid network of rapidly increasing EVs and dominated gasoline vehicles (GVs), which should help reduce the range anxiety of EV users and maximize the coverage rate of EVs.

According to previous studies, a variety of factors affect the location of charging stations, including, among others, preference and travel choice behavior of users [6,7], travel demand of users [8], information of the en-route energy consumption of EVs [9,10,11,12], information of the remained range of EVs [13,14], charging speed and demand of battery [7,15], road environment [5,16], land supply [17], elastic charging demands [18]. However, it is not practical to build a complicated location model and to solve the location problem by incorporating so many macro and micro factors into the planning and decision making of EV charging stations. This paper mainly studies the uncertainties due to the choice behavior of EV and GV users under elastic demands and uncertain path constraints. A bilevel model is proposed in this paper, where the upper model deals with a maximum flow-covering (MFC) problem, while the lower model is a stochastic user equilibrium model with elastic demands (SUE-ED). It is worth noting that factors, such as road congestion, charging speed, range limit of EVs, and capacity of charging stations, are already considered in users’ travel choices. The lower model considers the elastic demands and the distance constraints of EVs, which have a significant influence on route choice and the distribution of link flow. Finally, a generalized Lagrangian function is constructed to prove the equivalence between the stochastic user equilibrium model and the elastic demand model.

The remainder of this paper is structured as follows: Section 2 reviews the relevant literature. Section 3 affirms the problem hypothesis, analyzes three charging paths of the EVs, and gives the symbolic description used in the article. Section 4 establishes a bilevel model. Section 5 and Section 6 introduce the algorithms used to solve the model and numerical analysis applying to the Nguyen–Dupuis network, respectively. Finally, the research results and future research directions are discussed in Section 7.

2. Literature Review

The conventional and dominant three types of location optimization models include the point demand model [19,20], the flow demand model [20,21,22,23], and the multi-objective optimization model [24,25,26,27]. The P-Median model assumes that “charging demand is generated in the road network node” and is widely used as one of the point demand models. The flow demand model is formulated on the basis of the Flow Capturing Location Model (FCLM). For the first time, some researches proposed the Flow Refueling Location Model (FRLM), where mileage limitations were explicitly considered in facility location issues [28,29,30,31]. One branch of FRLM aims to maximize the demand coverage by locating a fixed number of charging facilities, which is called the maximum coverage location problem. Although the multi-objective optimization models have the advantage of addressing more complicated experimental requirements, they are not good at dealing with the uncertainty planning problem [32].

In a hybrid network with both EVs and GVs, the layout of EV charging stations and the distribution of EV flow affect mutually. Elastic demands in a hybrid network result in many more uncertainties. Most of the previous studies used the User Equilibrium Model (UE) to locate charging facilities. In these studies, Xu et al. dealt with the user equilibrium problem in a hybrid transport network with battery switching stations and road grade constraints [33]. Jing et al. gave a comprehensive discussion of the equilibrium network model [34]. Jiang et al. introduced the path distance constraints into the UE model [13]. Zheng et al. proposed a bilevel model where the upper layer minimized travel costs, and the lower layer aimed to find the path-constrained EVs equalization flow [35]. The bounded rationality of EV users led to much more complicated energy consumption [7,36,37], and, therefore, route choice behavior reflected more unobserved heterogeneity, which resulted in various elastic demands.

The network design problem with elastic demands, where the induced or transferred OD demand is the subject of responses in traveler itinerary choices to enjoy the improvement of new infrastructure, have several formulations with various motivations. Ge et al. [38] is one of the early attempts to consider both the proportion of EVs and the charge rate of EVs when determining the elastic charging demands from the total number of vehicles on road connections. The elastic demand was formulated based either on feedback of congested travel and congested station on route choices [39] or on the assumption that charge demand between OD pairs follows a nonlinear inverse cost function without considering the pre-generating paths and charging combinations [40].

It should be noted that it is hard to consider all these constraints simultaneously, say, the elastic demand of the road network, the capacity of the charging stations, and the range limit of EVs by using the SUE-ED model in location problems of public EV charging stations. In this study, a novel bilevel public charging station location model that combines a flow-capturing location model and a multi-objective optimization model is proposed.

3. Problem Assumptions and EV Paths Analysis

The proposed location model has two optimization objectives. The upper layer employs an improved maximum coverage model to maximize the coverage of the total EV flow by deploying a given number of charging stations on the links where the EV flow is the largest. The charging stations can serve more EV users by deploying charging facilities on the links where most EV drivers pass, which is an effective way to improve the utilization rate of public charging facilities and alleviate range anxieties of electric vehicle drivers [4]. The lower layer uses an SUE-ED model with incomplete information. The assignment result of SUE is a decisive factor for the placement of the charging stations. There is an inherent difference in the driving behavior of each user (EV users and GV users), in the hybrid network. In particular, range limit, charging time, and location of charging facilities have an important impact on the path choose behaviors of EV users.

3.1. Notations

i: The set of traveler types in the network. i presents an element, mainly including EVs and GVs. When specific instructions are needed, subscripts g and e are used to indicate variables or parameters related to GVs and EVs, respectively;

W: The set of OD pairs. w is an element of the set;

$R{ }_w$: The set of paths between all OD pairs;

p: The number of charging stations subject to the budget;

U: The utility value of the charging stations, resulting in a reduction in travel cost;

$K·U$: The waiting cost in charging stations, where K is a constant coefficient;

$l_k^w$: Length of path k between OD pair w;

$l_k^{ri}$: Length of path k between node r and I;

$l_k^{ij}$: Length of path k between node i and j;

$l_k^{js}$: Length of path k between node j and s;

$R_e$: The range limit of EVs;

ε: Charging power. Unit: hour/mile;

$t_{ck}^w$: The minimum charging time;

$v_a$: Traffic flow on link a, including EVs and GVs, so $v_a=v_{ae}+v_{ag};$

$t_a$($v_{ag}$,$v_{ae}$): Actual travel time through link a;

$t_a^0$: Free travel time through link a;

$H_a$: Capacity of link a;

$x_a$: Binary parameter. If the charging facility is on link a, $x_a$ = 1, otherwise $x_a$ = 0;

${\delta }_{ak}^w$: Binary parameter. If the path k crosses the link a, ${\delta }_{ak}^w$ = 1, otherwise ${\delta }_{ak}^w$ = 0;

$q{ }_i^w$: Elastic travel demand of type-i traveler between OD pair w;

${\overline{q}}_i^w$: Maximum potential demand of type-i traveler between OD pair w;

$f{ }_{ki}^w$: Traffic flow of type-i traveler on path k between OD pair w;

$p{ }_{ki}^w$: The probability that type-i traveler chooses the path k between OD pair w;

$c{ }_k^w$: The actual travel time of the path k between OD pair w, and $c{ }_k^w={\displaystyle \sum }{t}_a( v{ }_a){\delta }_{ak}^w;$

${\overline{c}}_{ki}^w$: The generalized travel cost of type-i traveler on path k between OD pair w;

$C{ }_{wi}$: The expected perceived travel cost of the type-i traveler between OD pair w;

$D{ }_{wi}$(·): The traffic demand function of the type-i traveler between OD pair w;

$c^w(\overrightarrow{x})$: Vector of the actual travel cost of all paths between OD pair w;

${\theta }_i$: a non-negative parameter that characterizes the uncertainty of type-i traveler’s understanding of the path travel time.

3.2. Propose Assumptions

This paper focus on a hybrid network, denoted as G = (N,A), where N is the set of nodes, A is the set of links. W is used to represent a set of OD pairs, and w is an element, w ∈ W. $R{ }_w$ represents the set of paths between w. The main differences between GVs and EVs are the limitation of travel distance and the composition of travel costs. Without loss of generality, the following assumptions are made:

(a)

There are two types of users in the mixed network: GVs and EVs, where EVs have an identical range limit;

(b)

The travel demand of GVs and EVs between each OD pair is elastic. And two types of vehicles have incomplete information about the travel cost;

(c)

Each EV is fully charged at its origin;

(d)

The level of anxiety and risk-taking behaviors of EV drivers are not considered in this model;

(e)

The charging time of electric vehicles is linearly related to mileage;

(f)

The charging facilities will be placed in the middle of links whose EV flow ranks p;

(g)

Deploying the charging facility on the road/path will increase the attractiveness of the route, which is called the utility of the charging facilities U, reducing the path travel cost;

(h)

Charging facilities have a fixed charging capacity. It is assumed that the waiting cost is proportional to the attraction value, expressed by $K·U$.

3.3. Analysis of EV Paths

When performing SUE-ED traffic assignment, the GVs can choose any route because of no range limit. However, the selection of feasible path sets is required before the distribution of EV flow. To describe the three types of EV paths more clearly, the concept of sub-path is described as follows.

Suppose that the origin is r and the destination is s in a pair OD w, and there are two charging stations i, j, which are also regarded as nodes, located in the arc of the path k (Figure 1). When $k^{ri}$, $k^{ij}$, $k^{js}$ no longer includes other charging facilities except i, j, they are called sub-paths [41].

Based on the relationship between the travel mileage limit $R_e$ and the sub-path distance, the three travel paths of EVs are discussed, as shown in

Table 1. Three travel paths of electric vehicles (EVs).

Scenarios	Conditions	Travel Cost
1. A path that can be completed without charging	$l_k^w$ ≤ $R_e$	${\overline{c}}_{ke}^w=c_k^w$
2. A path that cannot be completed after charging	$l_k^{ri}$ ≥ $R_e$ or $l_k^{ij}{R}_e$ or $l_k^{js} R_e$	${\overline{c}}_{ke}^w=\infty$
3. A path that can be completed after charging	$l_k^w\ge R_e$ and $l_k^{ri}$ ≤ $R_e$ and $l_k^{ij}$ ≤ $R_e$ and $l_k^{js}$ ≤ $R_e$	${\overline{c}}_{ke}^w=c_k^w$ + $t_{ck}^w$ + (K − 1)U

:

Scenario 1. A travel path that can be completed without charging. When $l_k^w$ ≤ $R{ }_e$, the EV drivers can reach the destination without charging. The generalized travel cost is ${\overline{c}}_{ke}^w=c_k^w$ (no charging stations) or ${\overline{c}}_{ke}^w=c_k^w-U$ (with charging stations).

Scenario 2. A travel path that cannot be completed after charging. When $l_k^{ri}$ ≥ $R_e$ or $l_k^{ij}\ge R_e$ or $l_k^{js}\ge R_e$, the path k is an infeasible path. And the generalized path travel cost becomes infinite (${\overline{c}}_{ke}^w=\infty$) and the probability that the drivers select the path is zero.

Scenario 3. A travel path that can be completed after charging. The path k is defined as a feasible path, if $l_k^w\ge R_e$ and $l_k^{ri}$ ≤ $R_e$ and $l_k^{ij}$ ≤ $R_e$ and $l_k^{js}$ ≤ $R_e$. The drivers can complete the trip by charging at least once and the minimum charging time is $t_{ck}^w$ = ε…($l_k^w-R_e$). The generalized travel cost is expressed as ${\overline{c}}_{ke}^w$ = $c_k^w+t_{ck}^w-U+K·U$ = $c_k^w$ + $t_{ck}^w$ + (K − 1)U, consisting of four parts: travel time, charging time, charging facilities utility, and charging stations’ waiting cost.

For GVs, the general path travel cost is ${\overline{c}}_{kg}^w$ = $c_k^w$.

4. Establish a Double-Layer Model

In this section, a bilevel optimization model for the location of charging facilities is proposed. The upper level model determines the location of the charging facilities by selecting the top p links. The lower level model calculates users’ generalized travel cost, and randomly allocates the flow demand between the OD pairs to the filtered paths through known charging facilities location.

4.1. Upper Level Problem

The upper model is designed to maximize the total covered EV link flow by deploying a given number of charging facilities. The EV flow is covered when the charging facility is on the link. That is

$$\mathrm{Max}{{\displaystyle \sum }}_a{v}_{ae}{x}_a$$

(1)

$$\mathrm{Subject}\mathrm{to}{{\displaystyle \sum }}_{a\in A}{x}_a=p$$

(2)

Equation (1) is the objective function of the upper model, and Equation (2) is the budget constraint indicating the number of charging stations p in a given network.

4.2. Lower Level Problem

The Bureau of Public Road (BPR) function in the lower model is employed [42].

$$t_a(v_{ag}, v_{ae})=t_a^0\left\{1+0.15\times {(\frac{v_{ag}+v_{ae}}{H_a})}^4\right\}$$

(3)

Since $U$, $t_{ck}^w$, k are parameters that are independent of the equilibrium flow of the link, the Jacobian matrix [43] in the lower layer problem is as follows:

$$\frac{\partial {\overline{c}}_{kg}^w}{\partial v_{ae}}=\frac{\partial {\overline{c}}_{ke}^w}{\partial v_{ag}}=0.6{\displaystyle \sum }{t}_a^0{\delta }_{ak}^w\frac{{(v_{ag}+v_{ae})}^3}{H_a{}^4}$$

(4)

It proves that the Jacobian matrix is symmetrical, and the lower layer model can be established as a convex function problem. It is assumed that all types of travelers make path selection in a random manner. According to the random utility theory, the probability that type-i traveler chooses the path k between OD pair w is:

$$p{ }_{ki}^w=\frac{exp(-{\theta }_i{\overline{c}}_{ki}^w)}{{{\displaystyle \sum }}_{k\in R{ }_w}exp(-{\theta }_i{\overline{c}}_{ki}^w)}\forall k,w,i$$

(5)

Assume that elastic travel demand $q{ }_i^w$ of type-i traveler is a strictly monotonically decreasing function of the expected minimum travel time between OD pair w, and with an upper bound.

$$q{ }_i^w=D{ }_{wi}( C{ }_{wi})\le {\overline{q}}_i^w \forall w,i$$

(6)

According to the discrete selection theory [44,45], define $C{ }_{wi}$ as:

$$C{ }_{wi}(c^w(\overrightarrow{x})=\mathrm{E}\left[\underset{\mathrm{k}\in {\mathrm{R}}_{\mathrm{w}}}{\min }\left\{{\overline{c}}_{ki}^w\right\}| c^w(\overrightarrow{x})\right]=-\frac{1}{{\theta }_i}\ln {{\displaystyle \sum }}_{k\in R{ }_w}exp(-{\theta }_i{\overline{c}}_{ki}^w) \forall w,i$$

(7)

For the type-i traveler, its SUE condition can be expressed as:

$$f{ }_{ki}^w=q{ }_i^{w}p{ }_{ki}^w\forall k,w,i$$

(8)

In the mixed network, SUE-ED problem can be described by the equivalent mathematical programming model:

$$\mathrm{Min}\mathrm{Z}(x,f,q)={{\displaystyle \sum }}_{a\in A}{{\displaystyle \int }}_0^{v_a}{t}_a(w)\mathrm{d}w+{{\displaystyle \sum }}_{i\in I}\frac{1}{{\theta }_i}{{\displaystyle \sum }}_{w\in W}{{\displaystyle \sum }}_{k\in R{ }_w}f{}_{ki}^w(\ln f{}_{ki}^w-1)-{{\displaystyle \sum }}_{i\in I}\frac{1}{{\theta }_i}{{\displaystyle \sum }}_{w\in W}{{\displaystyle \sum }}_{k\in R{ }_w}q{}_i^w(\ln q{}_i^w-1)-{{\displaystyle \sum }}_{i\in I}{{\displaystyle \sum }}_{w\in W}{{\displaystyle \int }}_0^{q{ }_i^w}{D}_{wi}^{-1}(w)\mathrm{d}w$$

(9)

$$\mathrm{Subject}\mathrm{to}:q{}_i^w={{\displaystyle \sum }}_{k}f{}_{ki}^w,\forall w,i$$

(10)

$$f{}_{ki}^w\ge 0,\forall k,w,i$$

(11)

$$q{ }_i^w\ge 0,\forall w,i$$

(12)

$$v_{ai}={{\displaystyle \sum }}_{w\in W}{{\displaystyle \sum }}_{k\in R_w}f{ }_{ki}^w{\delta }_{ak}^w,\forall a,i$$

(13)

$$R_e-l_k^{ij}\ge 0, \forall k$$

(14)

Equation (10) is the flow conservation constraint. Equation (11) is the non-negative constraint of path flow for type-i travelers. Equation (12) is the non-negative constraint of type-i traveler’s OD flow demand. Equation (13) is the correlation between link flow and path flow. The novelty of this problem is that the introduction of sub-paths in the Equation (14) can generate a feasible set of paths in advance from the finite paths between each OD pair. The superscripts i and j include the origin r and the destination s of all OD pairs in Equation (14).

It is necessary to prove the equivalence between the solution of the proposed program Equation (9) and the solution of the SUE-ED model. The generalized Lagrangian function [46] of the mathematical model is constructed as follows:

$$\mathrm{L}=\mathrm{Z}+{{\displaystyle \sum }}_i{{\displaystyle \sum }}_w{\lambda }_{wi}(q{ }_i^w-{{\displaystyle \sum }}_{k\in R_w}f{ }_{ki}^w) -{{\displaystyle \sum }}_i{{\displaystyle \sum }}_w{{\displaystyle \sum }}_k{\mu }_{ki}^{w}f{ }_{ki}^w-{{\displaystyle \sum }}_{i\in I}{{\displaystyle \sum }}_{w\in W}{\nu }_{wi}q{ }_i^w-{{\displaystyle \sum }}_w{{\displaystyle \sum }}_k{\alpha }_{ke}^{ij}\ast (R_e-l_k^{ij})$$

(15)

${\lambda }_{wi}$, ${\mu }_{ki}^w$, ${\nu }_{wi}$, ${\alpha }_{ke}^{ij}$ are the Lagrange multipliers.

$$\begin{array}{cc}\mathrm{If}{l}_k^{ij}\le R_e& {\alpha }_{ke }^{ij}=0\\ \mathrm{If}{l}_k^{ij}\ge R_e& {\alpha }_{ke}^{ij}\ge 0\end{array}$$

(16)

According to the Kuhn Tucker conditions [47], Equation (15) must satisfy the following conditions at the extreme point:

$$\frac{\partial L}{\partial f{ }_{ki}^w} =0,{\mu }_{ki}^{w}f{ }_{ki}^w=0,{\mu }_{ki}^w\ge 0 \forall k,w,i$$

(17)

$$\frac{\partial L}{\partial q{ }_i^w} =0,{\nu }_{wi}q{ }_i^w=0,{\nu }_{wi}\ge 0 \forall w,i$$

(18)

The partial derivative of $f{ }_{ki}^w$ for Equation (15) is derived as follows:

$$\frac{\partial L}{\partial f{ }_{ki}^w}={{\displaystyle \sum }}_a{t}_a(v_a){\delta }_{ak}^w +\frac{1}{{\theta }_i}\ln f{ }_{ki}^w-{\lambda }_{wi}-{\mu }_{ki}^w=0 \forall k,w,i$$

(19)

Note that if $f{ }_{ki}^w$ = 0, then $\frac{\partial L}{\partial f{ }_{ki}^w}$ does not exist. The above formula is only valid when $f{ }_{ki}^w>$ 0. So ${\mu }_{ki}^w$ = 0.

$$f{ }_{ki}^w=\exp \left[-{\theta }_i({\overline{c}}_{ki}^w-{\lambda }_{wi})\right]=\exp \left[-{\theta }_i{\overline{c}}_{ki}^w\right]·\exp \left[-{\theta }_i{\lambda }_{wi}\right]$$

(20)

$$q{ }_i^w=\exp (-{\theta }_i{\lambda }_{wi}){\displaystyle \sum }_{k\in R_w}\exp (-{\theta }_i{\overline{c}}_{ki}^w)$$

(21)

$$p{ }_{ki}^w=\frac{exp(-{\theta }_i{\overline{c}}_{ki}^w)}{{{\displaystyle \sum }}_{k\in R{ }_w}exp(-{\theta }_i{\overline{c}}_{ki}^w)}\forall k,w,i$$

(22)

Equation (22) shows that type-i traveler follows the Logit model to select the travel path, which satisfies the SUE-ED condition described in Equation (8).

Calculate the partial derivative of $q{ }_i^w$ for Equation (15) as follows:

$$\frac{\partial L}{\partial q{ }_i^w}=-\frac{1}{{\theta }_i}\ln q{ }_i^w-D_{wi}^{-1}(q{ }_i^w)+{\lambda }_{wi}-{\nu }_{wi}=0\forall w,i$$

(23)

Because $\ln q{ }_i^w$ should exist, so ${\nu }_{wi}$ = 0. Then derive from Equation (21):

$${\theta }_i{\lambda }_{wi}=\ln q{ }_i^w-\ln {\displaystyle \sum }_{r\in R_w}\exp \left[-{\theta }_i{\overline{c}}_{ri}^w\right]$$

(24)

Comparing Equation (23) and Equation (24), the inverse function of the flow demand function is as follows:

$$D_{wi}^{-1}(q{ }_i^w)=-\frac{1}{{\theta }_i}\ln {\displaystyle \sum }_{r\in R_w}\exp \left[-{\theta }_i{\overline{c}}_{ri}^w\right]\forall w,i$$

(25)

Then $q{ }_i^w=D_{wi}(C_{wi})$. This shows that Equations (9)–(14) can be used to represent a multi-user SUE problem under elastic demand. The BPR function $t_a$($v_a)$ is a strictly monotonically increasing function of the link flow $v_a$. The objective function is a strict convex function about the link flow vector v and the path flow vector f. At the same time, the constraints of Equations (9)–(14) are linear equality constraints and non-negative constraints, so its solution space is a convex set. According to the optimization theory, the strict convex function defined on the convex set only has one optimal solution.

5. Solution Method

The proposed bilevel model is an NP-hard (non-deterministic polynomial hard) problem because the lower level model is a mixed-integer nonlinear program with nonconvex path choosing, subject to uncertainties in elastic demands, where it is difficult to find a polynomial time complexity algorithm [48]. Thus, a heuristic method is employed to solve the problem of SUE-ED and MFC in an iterative manner. The detailed processes in Figure 2 are as follows:

Step 1: Input the road network parameters, such as the length of the link and the capacity of the links, and find all the paths between the OD pair w, and record them as the initial path set.

Step 2: Initialize. There is no charging facility in the road network and relax the range limit of EVs. Based on the travel time of the zero-flow link $t_a(0)=t_a^{ 0}$, calculate the effective path impedance $c{ }_k^w$, and obtain the initial link flows $V_{ag}(1)$ and $V_{ae}(1)$. Let the upper generation counter z = 1.

Step 3: Sort all the links flow of EVs in ascending order to find the top p ranked links, and arrange the charging facilities in the middle of the p links. Then increase the upper iteration counter by 1.

Step 4: Perform SUE-ED traffic distribution on the network with the location of the charging facilities obtained by step 3. The detailed steps are as follows:

Step 4.1: Check feasible paths and update the path set. If the length of any sub-path is greater than $R_e$, remove the path from the initial set of paths. The program stops if there is no feasible path between any OD pair. If there is at least one feasible path between each OD pair, proceed to the next step.

Step 4.2: Carry out random loading of traffic flow in the network with charging facilities. Calculate the generalized path travel cost ${\overline{c}}_{ke}^w$, elastic flow demand between OD pairs, the probability that the path is selected to obtain updated road link flows $V_{ag}(2)$ and $V_{ae}(2).$ Then set the iteration counter n = 1.

Step 4.3: Repeat the random loading of step 4.2 to obtain additional link flow {$y_{ag}$}, {$y_{ae}$}.

Step 4.4: Use the predetermined step size sequence {${\alpha }_n$}: ${}_a=1/n$, n = 1,2,…,∞.

Step 4.5: Calculate the current flow of links by Method of Successive Average [49].

$V_{ag}^{n+1}$ = $V_{ag}^n$ + ($\frac{1}{n}$)($y_{ag}^n$ − $V_{ag}^n$); $V_{ae}^{n+1}$ = $V_{ae}^n$ + ($\frac{1}{n}$)($y_{ae}^n$ − $V_{ae}^n$).

Step 4.6: If the convergence condition is met $\frac{\sqrt{{{\displaystyle \sum }}_a{({\overline{v}}_a^{n+1}-{\overline{v}}_a^n)}^2}}{{{\displaystyle \sum }}_a{\overline{v}}_a^n}\le \omega$, where ω indicates the convergence accuracy, then proceed to step 5. {$V_{ag}^{n+1}$}, {$V_{ae}^{n+1}$} are the sets of balanced link flow for the GVs and EVs, respectively. Otherwise, set n = n + 1 and go to step 4.3.

Step 5: Repeat step 3 and update charging facilities’ location. The program stops until the location of the charging facilities is no longer changed; that is, the maximum coverage flow remains stable; otherwise, go to step 4.

6. Numerical Analysis

The model is applied to the Nguyen–Dupuis network [32,45] for a case study in this section, which has been widely used in transportation network researches in the past decades. Due to the small scale of the Nguyen–Dupuis network, the paths can be enumerated to better analyze the relationship between the location of the charging station and the feasible path and the traffic volume of the section. At the same time, the elastic demand between OD pairs, the effects of charging speed, range limitation, charging facilities utility, and waiting cost on the location of charging facilities were evaluated. The test network (shown in Figure 3) consisted of 13 nodes, 19 segments, 25 paths (shown in

Table 2. Path compositions and lengths in the Nguyen–Dupuis network.

OD Pair	Path Number	Path Composition	Length
(1,2)	1	1-5-6-7-8-2	29
	2	1-12-8-2	32
	3	1-5-6-7-11-12	33
	4	1-12-6-7-8-2	35
	5	1-5-6-10-11-2	38
	6	1-12-6-7-11-2	39
	7	1-5-9-10-11-2	41
	8	1-12-6-10-11-2	44
(1.3)	9	1-5-6-7-11-3	32
	10	1-5-9-13-3	36
	11	1-5-6-10-11-3	37
	12	1-12-6-7-11-3	38
	13	1-5-9-10-11-3	40
	14	1-12-6-10-11-3	43
(4,2)	15	4-5-6-7-8-2	31
	16	4-5-6-7-11-2	35
	17	4-9-10-11-2	37
	18	4-5-6-10-11-2	40
	19	4-5-9-10-11-2	43
(4,3)	20	4-9-13-3	32
	21	4-5-6-7-11-3	34
	22	4-9-10-11-3	36
	23	4-5-9-13-3	38
	24	4-5-6-10-11-3	39
	25	4-5-9-10-11-3	42

), and 4 OD pairs (1, 3), (1, 2), (4, 3), (4, 2). The traveler’s OD demand function uses a linear function. The familiarity of the two types of travelers to the network was measured by the parameters ${\theta }_e$ and ${\theta }_g,$ respectively. The parameters were set as follows during the trial calculation. $q_e^w( C{ }_{we})=400-7 C{ }_{we}$, $q_g^w( C{ }_{wg})=400-7 C{ }_{wg}$, ${\theta }_e$ = 0.1, ${\theta }_g$ = 0.1, U = 5, $R_e$ = 20, ε = 1, K = 0.5, p = 3. The trial results are shown in

Table 3. The charging facility locations and Electric Vehicle (EV) flows

Link Node	1		2		3		4	20
	Location	EV Flow	Location	EV Flow	Location	EV Flow	Location	EV Flow	Location
(1,5)	×	367.5	×	219.3	√	197.4	×	172.5	×
(4,5)	×	364.9	×	195.4	×	172	×	143.7	×
(4,9)	×	164.5	×	54.8	×	27.4	×	0	×
(5,6)	×	538.5	√	350.1	√	337.1	√	316.2	√
(5,9)	×	193.9	×	64.6	×	32.3	×	0	×
(6,7)	×	499.7	√	337.3	√	363.7	√	380	√
(7,8)	×	196.7	×	180.5	×	210.3	×	233.2	×
(8,2)	×	248.7	×	197.8	×	219	√	233.2	√
(1,12)	×	196.7	×	100.5	×	83.2	×	63.8	×
(11,3)	×	380.1	×	200.4	×	175.2	×	146.8	×
(12,6)	×	144.8	×	83.2	×	74.6	×	63.8	×
(12,8)	×	51.9	×	17.3	×	8.7	×	0	×
(13,3)	×	157.1	×	52.4	×	26.2	×	0	×
(6,10)	×	183.7	×	95.9	×	48	×	0	×
(7,11)	×	302.9	×	156.8	×	153.4	×	146.8	×
(9,10)	×	201.4	×	67.1	×	33.6	×	0	×
(9,13)	×	157.1	×	52.4	×	26.2	×	0	×
(10,11)	×	385.1	√	163	×	81.5	×	0	×
(11,2)	×	308	×	119.5	×	59.7	×	0	×

.

Table 3 shows the results of 20 iterations. The total coverage of the charging stations was 0, 906.7, 919.8, 924.9…929.5. In the first iteration, only six paths were available among the four OD pairs, thereby causing a change in the probability of path selection. It is not difficult to find that the EV link flow tended to be zero in the infeasible paths, and the EV flows gradually reached equilibrium in the links with charging facilities at last. The final locations of the charging facilities were (5,6); (6,7); (8,2). The changes in the location of charging stations (CSs) and feasible paths are listed in

Table 4. The changes in the location of charging stations (CSs) and feasible paths.

OD Pair	Location 1	Location 2	Location 3
	(5,6); (6,7); (10,11)	(5,6); (6,7); (1,5)	(5,6); (6,7); (8,2)
(1,2)	Three paths, including 1,4,5	Two paths, including 1,4	Two paths, including 1,4
(1,3)	Three paths, including 9,11,12	Two paths, including 9,12	Two paths, including 9,12
(4,2)	Two paths, including 15,18	One path, 15	One path, 15
(4,3)	Two paths, including 21, 24	One path, 21	One path, 21

. After the first iteration, the feasible paths for the EVs were reduced because of range limitations. Among the four OD pairs, only six paths were left to choose, which led to a change in the probability of path selection, so the charging facility location also changed. The reason for the change in CSs from (5,6); (6,7); (1,5) to (5,6); (6,7); (8,2) is mainly due to the fact that there were three paths, including link (8,2), but just one path including link(1,5). In particular, only one feasible path, including the link (8,2), was left in the OD pair (4,2), so all EV users could only select this path.

The elastic demand analysis between OD pairs (1,3) is shown in Figure 4. Traffic demands are mainly affected by travel costs. The maximum demand for both electric cars and gasoline were set as 400. At the beginning, there were no CSs on the network, and the range limit of EVs was relaxed, so the flow demands of EVs and GVs were the same, 265.81. In the second iteration, the reduction in travel demand for EVs and GVs was mainly due to the increase in travel costs caused by the increased link flow. For EVs, the feasible paths were less than GVs because of the range limit, and the travel cost was not only affected by the link flow as GVs, but also affected by the utility of charging facilities, charging speed, and waiting cost at charging stations. Therefore, the travel demand of EVs was much lower than GVs, which inspired us to increase the range of EVs and the accessibility of charging facilities to promote the travel demand and development of EVs. The elastic fluctuation of flow demands affect the distribution of link flow and further affect the travel cost and the location of charging facilities, which will react to the travel demand. In the continuous iteration and mutual influence, the final equilibrium was reached. The change trend of elastic demand in other OD pairs was consistent with the OD pair (1,3), which was matched with the actual travel demand.

The results of sensitivity analysis about range limits, charging speed, charging facilities utility, and charging stations waiting cost are shown in Figure 5, Figure 6, Figure 7 and Figure 8.

In Figure 5, the impact of EV range limit on the locations of charging facilities is revealed. When the range limits change, the charging facilities were located differently. As the range limit increased, there were more feasible paths within a certain range, and the EV covered flow first increases and then decreased, and finally increases when $R_e=45$.

Then the effect of charging speed on the location of the charging facilities is examined in Figure 6. Formula $t_{ck}^w=\epsilon ·(l_k^w-R_e)$ indicates that the charging speed directly affects the charging time and general travel cost, which results in different charging facility layouts. Different values of ε can present different charging methods. ε = 0.1 means fast charging and ε = 10 means slow charging. EV drivers tend to choose a path with a long travel distance without charging facilities, rather than a path with a very slow charging speed.

In Figure 7, the utility of the charging facilities reflects the attractiveness of the link for EV users. The greater the utility value, the higher the level of anxiety for EV users, so they are more likely to choose a path with charging facilities. If there are many types of EVs in the network, the utility value may not be very large for EVs with large battery capacity.

The relationship between the charging stations waiting cost coefficient K and the charging station position is shown in Figure 8. When 0 < K < 1, the charging cost is reduced. When K > 1, the charging cost increases, and the total coverage of the electric vehicle begins to decrease. When deploying a fast charging station, a slightly larger capacity should be considered to reduce the occurrence of K > 1.

7. Conclusions

This paper studies the deployment of public charging stations in a hybrid network to maximize the service efficiency of the charging facilities. A bilevel model was proposed to depict the interaction between the mixed link flow and the location of the charging stations. The link flow obtained by the lower SUE-ED model is the key to determine the location of the CSs. In the upper level model, the CSs were arranged on the link flow ranking p to achieve the maximum coverage. Four important factors were taken into accounts in the lower level SUE model, including the range limit of EVs that affected the path choice, the elastic travel demand closely related to the distribution of link flow, the road congestion effect, and the capacity of charging facilities in travel costs. These four elements interact with each other and continuously iterate to reach equilibrium, which was more consistent with the actual travel situation. A hybrid integer nonlinear programming method based on the method of successive average (MSA) was constructed to prove the equivalence and the uniqueness of the SUE-ED model with range constraints. Finally, a network trial was conducted to examine the impact of elastic demand between OD pairs, the range limit, charging speed, the charging facilities’ utility, and waiting cost on the location problem.

It should be noted that the actual road network systems are very diverse, especially in urban areas, and significantly differ from the example of the Nguyen–Dupuis Network. More realistic factors are required to be considered in the future when designing the location of charging stations, not only by technical and operational factors but also by social factors. And assumptions can be relaxed appropriately in future work. The nonlinearity of the charging time, the uncertainties in EVs energy consumption, as well as the bounded rationality of EV travelers, should be considered.