**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:

**Figure 2.** Framework of the bilevel model to solve the optimal location of charging stations (CSs).

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 *ta*(0) = *t* <sup>0</sup> *<sup>a</sup>* , calculate the effective path impedance *c <sup>w</sup> <sup>k</sup>* , and obtain the initial link flows *Vag*(1) and *Vae*(1). Let the upper generation counter *<sup>z</sup>* <sup>=</sup> 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 *Re*, 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 *c w ke*, elastic flow demand between OD pairs, the probability that the path is selected to obtain updated road link flows *Vag*(2) and *Vae*(2). Then set the iteration counter *n* = 1.

Step 4.3: Repeat the random loading of step 4.2 to obtain additional link flow {*yag*}, {*yae*}.

Step 4.4: Use the predetermined step size sequence {α*n*}: *<sup>a</sup>* = 1/*n, n* = 1,2, ... ,∞.

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

$$V\_{a\underline{\mathcal{S}}}^{n+1} = V\_{a\underline{\mathcal{S}}}^{n} + (\frac{1}{n})(y\_{a\underline{\mathcal{S}}}^{n} - V\_{a\underline{\mathcal{S}}}^{n});\ V\_{a\underline{\mathcal{e}}}^{n+1} = V\_{a\underline{\mathcal{e}}}^{n} + (\frac{1}{n})(y\_{a\underline{\mathcal{e}}}^{n} - V\_{a\underline{\mathcal{e}}}^{n}).$$

Step 4.6: If the convergence condition is met *<sup>a</sup>* (*vn*<sup>+</sup><sup>1</sup> *<sup>a</sup>* <sup>−</sup>*vn a* ) 2 *<sup>a</sup> vn a* ≤ ω, where ω indicates the convergence accuracy, then proceed to step 5. {*Vn*+<sup>1</sup> *ag* }, {*Vn*+<sup>1</sup> *ae* } 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), 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 θ*<sup>e</sup>* and θ*g*, respectively. The parameters were set as follows during the trial calculation. *qw <sup>e</sup>* ( *C we*) = 400 − 7 *C we*, *qw <sup>g</sup>* ( *C wg*) = 400 − 7 *C wg*, θ*<sup>e</sup>* = 0.1, θ*<sup>g</sup>* = 0.1, *U* = 5, *Re* = 20, ε = 1, *K* = 0.5, *p* = 3. The trial results are shown in Table 3.

**Figure 3.** Nguyen–Dupuis network.




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

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. 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.

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


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.

**Figure 4.** The flow demands of origin-destination (OD) pair (1,3).

The results of sensitivity analysis about range limits, charging speed, charging facilities utility, and charging stations waiting cost are shown in Figures 5–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 *Re* = 45.

**Figure 5.** The sensitivity analysis of the range limit.

Then the effect of charging speed on the location of the charging facilities is examined in Figure 6. Formula *t w ck* = <sup>ε</sup>·(*<sup>l</sup> w <sup>k</sup>* <sup>−</sup> *Re*) 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.

**Figure 6.** The sensitivity analysis of 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.

**Figure 7.** The sensitivity analysis of charging facilities' utility.

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.

**Figure 8.** The sensitivity analysis of waiting cost coefficient *K*.

### **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.

**Author Contributions:** Conceptualization, H.G. and K.L.; methodology, H.G. and K.L.; software, H.G. and X.P.; validation, H.G., K.L. and X.P.; formal analysis, H.G., K.L. and X.P.; investigation, H.G., K.L. and C.L.; resources, C.L.; writing—original draft preparation, H.G. and K.L.; writing—review and editing, H.G. and K.L.; visualization, H.G.; supervision, K.L. and C.L.; project administration, K.L.; funding acquisition, K.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant Nos. 51378091 and 71871043) and the open project of the Key Laboratory of Advanced Urban Public Transportation Science, Ministry of Transport, PRC.

**Acknowledgments:** This work was carried out by the joint research program of the Institute of Materials and Systems for Sustainability, Nagoya University.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
