# Electric Vehicle Routing Problems with Stochastic Demands and Dynamic Remedial Measures

Authors: Xianlong Ge; Ziqiang Zhu; Yuanzhi Jin
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

With the rapid development of e-commerce, logistic enterprises must better predict customer demand to improve distribution efficiency, so as to deliver goods in advance, which makes logistics stochastic and dynamic. In order to deal with this challenge and respond to the concept of “green logistics,” an electric vehicle routing problem with stochastic demands (EVRPSD) and proactive remedial measures is investigated, and an EVRPSD model with probability constraints is established. At the same time, a hybrid heuristic algorithm, combining a saving method and an improved Tabu search algorithm, is proposed to solve the model. Moreover, two insertion strategies with the greedy algorithm for charging stations and dynamic nodes are introduced. Finally, a large number of experimental data show that the heuristic algorithm proposed in this paper is feasible and effective.


---

## Body

## 1. Introduction

With the rapid development of e-commerce, the order quantity of logistic enterprises has increased rapidly, and the customers have higher requirements on the timeliness of delivery than before. The traditional collaborative pattern for processing orders (consumers, e-commerce platforms, sellers, and logistic enterprises) has been difficult to meet the requirements under a large number of orders, and it is urgent to transform the current pattern of logistics.A huge number of orders means that a lot of demand-related data can be obtained; therefore, e-commerce platforms have the ability of predicting consumer demand before the big promotion, and then logistic enterprises can transport the goods to the area near the target customers in advance. Once customers submit their order, the corresponding goods will be delivered from the nearest warehouse immediately, and it can greatly shorten the delivery time, realizing the immediate delivery. At present, some logistic enterprises have started to implement this pattern by establishing “front warehouses” that are near to the end customers, whose service radius is 3–5 kilometers, so it can shorten the delivery time from days to hours or even minutes.Front warehouses are characterized by the stochastic demand, which can be predicted by historical data, while the real demand is confirmed until a delivery vehicle arrives at a front warehouse. Therefore, the probability of insufficient goods may occur, and the vehicle may return to the depot for replenishment during the distribution process. To investigate this scenario, scholars proposed thecapacitated vehicle routing problem with stochastic demands (CVRPSD) [1]. In addition, when the vehicle starts from a regional warehouse to a front warehouse, there may be a replenishment request from new front warehouses, whose position and demand cannot be predicted in advance, and such a demand is defined as dynamic. The distribution center constantly updates the routing plan according to the emerging dynamic demand, which is defined as the dynamic demand vehicle routing problem (DVRP) [2]. Ecological logistics is another issue we need to consider toward the current serious environmental problems. At present, the general solution is to use electric vehicles instead of traditional fuel vehicles in the logistic industry. However, due to the limited battery capacity, electric vehicles often need to be charged or change battery on the way one or more times to successfully complete the whole task, which will increase the total cost to some extent. Therefore, it is necessary to develop a reasonable charging/changing strategy to reduce the distribution cost, and a new problem, electric vehicle routing problem (EVRP), is emerged [3].In order to meet the needs of real delivery scenarios, we will consider both the randomness and dynamics of customers and the mileage constraints of vehicles, thus proposing a new variant, namely, electric vehicle routing problems with stochastic demands and dynamic remedial measures. However, it should be noted that we have not considered the customer’s time window and the scenario of multiple distribution centers.This paper is organized as follows. In Section2, we will extend from the above three directions. In Section 3, we briefly describe the electric vehicle problem and give some general assumptions. In Section 4, we introduce an EVRPSD mathematical model in detail and propose an evaluation method for battery swap station insertion as well as the strategy of dynamic customer insertion. And, the detailed implementation of the CW-TS heuristic algorithm for solving the model is described in Section 5. The computational results are reported in Section 6. Finally, Section 7 contains a summary and future works.


## 2. Literature Review

There are a large body of literatures focusing onCVRPSD. Stewart and Golden [1] established a mathematical model based on the opportunity constraint and binary possibility theory for CVRPSD. The opportunity constraint model stipulated that a stochastic customer should be served only with an acceptable probability of failure, regardless of the additional cost caused by service failure, while the binary possibility theory considered the opportunity cost caused by service failure. Dror et al. [4] transformed the stochastic problem into a series of determined problems through a phased approach and obtained the solution of the problem under the consideration of the number of service failures.The methods to solve theCVRPSD model are divided into exact algorithms and heuristic algorithms. For the former, Gauvin et al. [5] designed a new branch-cut-and-price algorithm to solve CVRPSD, and a Tabu search heuristic and a bidirectional labeling algorithm were utilized to accelerate the pricing stage. Gendreau et al. [6] proposed an L-shaped algorithm for CVRPSD with uncertain customer location, which takes the location of the depot and the tour plan as decision variables. Numerical results shown that the algorithm was feasible and effective by performing on an instance containing 100 nodes. For the latter, Lei et al. [7] used the adaptive large neighborhood search algorithm to solve the routing problems of uncertain demand with time windows. Bertazzi and Secomandi [8] introduced a faster rollout search for the vehicle routing problem with stochastic demands and restocking and designed a quick rollout algorithm to solve the problem. In order to reduce the running time, they also proposed a new method to estimate the expected cost of the route when executing any VRP rollout algorithm with replenishment.The existing literature onDDVRP mainly focuses on the dynamic problem with uncertain customer demands. Qureshi et al. [9] solved DDVRP with soft time windows by means of the column generation algorithm, updating route status by using an event update strategy dynamically. Abdallah et al. [10] constructed a DDVRP model based on the capacity constraint and proposed an improved genetic algorithm to boost the solution quality according to the optimal combination of relevant parameters of selection, crossover, and mutation operators. Kuo et al. [11] established a mathematical model aiming at maximizing the number of serviced customers and minimizing the average customer waiting time for DDVRP with uncertain service time and proposed a regular update strategy to transform the dynamic problem into a series of static problems for the routing plan. Ge et al. [12] built a mathematical model minimizing the total cost for the open DDVRP and proposed a cloud genetic algorithm with a faster convergence rate than its classic counterpart to solve the DDVRP model. At the same time, the authors also introduced the concept of the key point for updating the routes dynamically. Okulewicz and Mańdziuk [13] extended DDVRP to the dynamic demand vehicle routing problem with multidistribution centers and adopted a two-stage multiswarm particle swarm optimization (2MPSO) heuristic to solve the model. Ghannadpour et al. [14] established a multiobjective DDVRP model with capacity and fuzzy time window constraints, aiming at minimizing the number of vehicles, total driving distance, vehicle waiting time, and maximizing user satisfaction. In addition, the genetic algorithm was used to solve the model. Ge et al. [15] studied the data-driven proactive vehicle routing optimization based on the historical distribution data of logistic companies, and four models were established for predicting customer demands, customer clustering, proactive demand quotas, and replenishment strategies.Research studies onEVRP mainly focus on vehicle load, battery capacity, and time window constraints and on the battery charging or swap strategy. Conrad and Figliozzi [16] proposed an EVRP-C to solve the “range anxiety” problem. Erdoğan and Miller-Hooks [17] also proposed an EVRP-C model, which assumes that the charging time and the service time of each customer are fixed and created feasible routes by using Clark–Wright’s saving method to insert the charging stations. Schneider et al. [18] extended the green vehicle routing problem (G-VRP) to the electric vehicle routing problems with time windows (EVRP-TW). They assumed that the vehicle visited the charging station at most once during the transportation process and established a mathematical model, in which the objective is minimizing the number of vehicles and the total travel distance. Lin et al. [19] proposed a partial charging strategy for EVRP-TW, i.e., the vehicle’s charging duration and location are the decision variables. Alejandro and Christelle [20] introduced a piecewise linearization function to fit the nonlinear charging process. Compared with the traditional linearization charging function, the new function reduced the calculation error of charging time. Based on the traditional linear discharging model, Goeke and Schneider [21] proposed a nonlinear battery energy consumption model based on real travel distance, road slope, vehicle speed, and load. Yang and Sun [22] introduced a BSS-EV-LRP model based on the location of battery swap stations. It included both the locating problems and the routing optimization under the limit of battery capacity, and they designed a hybrid approach called TS-MCWS which combines Tabu search and modified the Clarke and Wright saving algorithm to jointly solve the problem. In this problem, TS is designed to search the location strategy and a modified Clarke–Wright saving method undertakes the routing decision based on this location solution. On this basis, Hof et al. [23] designed an adaptive variable neighborhood search (AVNS) algorithm for the same problem. Numerical analysis showed that the method was robust and could obtain solutions with higher solution quality than SIGALN and TS-MCWS in most cases. Shao et al. [24] studied an electric vehicle routing problem with charging time and variable travel time. For obtaining the routes, the vehicle departure time at the depot, and the charging plan, they used the genetic algorithm, and to reflect a dynamic traffic environment, they implemented fluctuations in travel time.To sum up, scholars have contributed a lot of research on dynamic vehicle routing and electric vehicle routing [25–27]. However, with the rapid development of the logistic industry, the realistic logistic scenarios, such as cash-in-transit service of banks, product oil distribution, beer supply, and garbage collecting, need to take into account the uncertainty and environmental friendliness during the distribution process, and the research studies mentioned above lack techniques and strategies to solve such real cases. Therefore, this paper studies the routing optimization problems of s with stochastic demands and dynamic remedial measures. The contributions of this paper are the following four-fold:(1)
Formulating anEVRPSD model with the opportunity constraint(2)
Establishing an effective remedy strategy to fully utilize the residual load(3)
Designing a two-stage heuristic algorithm to solve theEVRPSD model(4)
Enriching the display of the vehicle route and thus improving the decision maker’s experience


## 3. Problem Description

The problem studied in this paper can be described as follows: the customer demand is divided into stochastic demand (the location of the customer is determined, and the demand is subject to a probability distribution) and dynamic demand (the time, location, and demand of the customer can only be confirmed if and only if the customer requests have arrived), and a distribution plan consists of two cycles.The first cycle: the predicted demand of each stochastic customer is obtained according to the historical demand data in the first stage, and the initial solution is generated according to the prediction; the vehicle starts from the distribution center and visits stochastic customers in sequence according to the initial solution, while the information center will receive new requests (generated by dynamic customers) in the time window [T, T0] during the transportation process; in the second stage, the actual demand of customer i (denoted qi) is attained when the vehicle k arrives at the location along with route rk, and then the remainder capacity of k is determined as to whether qi is satisfied. If not, the vehicle k returns to the distribution center after its rest of the goods are uploaded to i (consider this customer point as a “failure point”), and the remainder demand of customer i will be assigned to the next cycle. Judging the remainder capacity of k whether it satisfies the demand of the next stochastic customer i + 1, if the probability of this estimate is less than the risk preference α, the vehicle k will go to the customer i + 1; otherwise, stochastic customers after i (consider this customer point as a “remedial point”) in route rk will be assigned into the second cycle distribution. According to the vehicle’s current capacity and dynamic customer information collected by the information center, dynamic customers would be appropriately inserted to the current route rk and the vehicle will return to the distribution center after completing its tasks. In addition, if the vehicle has rest capacity too before it successfully serves the last customer j of the route, the insert operator can be repeated until its rest capacity is close to zero, and we call j as the “key point” in this condition. If there are dynamic customers that violate the capacity constraint, they will be assigned to the second cycle.The second cycle: the stochastic customers and dynamic customers unserved in the first cycle will be served. The objective is to minimize the total cost of two-cycle delivery (distance cost and fixed cost of vehicles). Since it is necessary to ensure that the electric vehicle has a positive power in the distribution process, the battery replacement plan will change with the adjusting route plan, which will further increase the complexity of the problem. Figure1 describes the whole process of distribution, in which stochastic customers 2, 4, and 5 and dynamic customer 3 are included in the second cycle.Figure 1
The illustration of the whole distribution process. (a) First cycle phase 1. (b) First cycle phase 2. (c) Second cycle.
(a)(b)(c)To further illustrate this problem, the following assumptions are made in this paper:(1)
Stochastic customer demand is a discrete stochastic variable, and its probability distribution function is known, and their demand, always constant, is confirmed if and only if the vehicle reaches the location of the stochastic customer(2)
The expected value of the stochastic customer is regarded as the predicted value that is not greater than the maximum capacity of the distribution vehicle simultaneously(3)
The information center only receives the dynamic demand within the time window [T, T0](4)
Only one distribution center is considered in the system, and all routes are closed(5)
The vehicle is fully loaded when it departs from the distribution center(6)
The style of power supply is replacement(7)
All vehicles are fully charged when leaving the distribution center or the electrical charging station, and the vehicles can be charged for many times during the transportation process(8)
The electrical supply station is divided into two types according to the situation of congestion, i.e., first-level station (relatively congested) and secondary-level station, whose number, geographical location, and congestion distribution are known(9)
The speed of the vehicle is constant


## 4. Mathematical Model

Since the traditional mathematical model for the electric vehicle routing problems cannot fit the problem proposed above, we extend the traditionalEVRP model. This section is divided into three parts. We integrated opportunity constraints into our model to limit the lower bound of probability of successfully serving the next customer. Also, based on the results by Lin et al. [19], we proposed an improved method for calculating the battery swap station insertion cost. Furthermore, we designed a new method for calculating the dynamic customer insertion cost to evaluate the optimal insertion location.


### 4.1.EVRPSD Model

In order to better introduce the model, the symbol description is given in Table1.Table 1
Notation used in this paper.
NotationMeaningIThe set of stochastic customers with indexiWThe set of dynamic customers with indexwBThe set of battery swap stations with indexbDThe set of distribution centers with indexoVThe set of vertexesKThe set of electric vehicles with indexkdijThe distance fromi to jqiThe actual demand of customerigwThe dynamic demand of customerwEqiThe expectation of customeri’s demandrPower consumption per kilometerαThe decision maker’s attitude towards riskpvk1Remaining battery capacity of vehiclek on arrival at vpvk2Remaining battery capacity of vehiclek on departing at vujk2The remaining load when the vehicle leaves pointjUVehicle capacityQBattery capacityC1Unit distance cost (¥/km) of electric vehiclesC2Fixed cost for each electric vehicletijTravel time (hours) on arc (i, j)speedTravel speed (km/h) of electric vehiclesxijk=1,0xijk=1 if the vehicle k flows on arc (i, j); otherwise, xijk=0yik=1,0yik=1 if the customer i is serviced by vehicle k; otherwise, yik=0xok=1,0xok=1 if the vehicle k returns to the depot after completing its tasks; otherwise, xok=0The mixed integer linear programming is established as follows:(1)minC1∑k∈K∑i∈V∑j∈V,i≠jxijkdij+C2∑k∈Kxok,s.t.∑i∈V,i≠jxijk−∑i∈V,i≠jxjik=0,∀j∈V,∀k∈K,(2)Pr∑i∈I∑j∈I,j≠ixijk⋅qj≤U≥α,α∈0,1,∀k∈K,(3)∑i∈IykiEqi≤U,∀k∈K,(4)∑k∈K∑i∈V,i≠oxoik≤k,∀o∈D,(5)∑i∈I∑j∈I,j≠ixijk≤I,∀k∈K,(6)pik1=pik2,∀i∈I,∀k∈K,(7)pik2=Q,∀i∈D∪B,∀k∈K,(8)pvk1≥0,∀v∈V,∀k∈K,(9)0≤pjk1≤Q−rdijxijk,μ,(10)0≤pjk1≤pik1−rdijxijk+Q1−xijk,i∈I,j∈V,j≠i,∀k∈K,(11)uok2=U,k∈K,o∈D,(12)0≤ujk2=ujk1−qj,j∈I,k∈K.Equation (1) represents the flow balance, which ensures that vehicles enter and exit from the same point. Equation (2) is a stochastic opportunity constraint, that is, only go to point j when the probability that the current remaining load can meet the demand at point j is greater than or equal to the risk preference α. This constraint can be transformed into a deterministic opportunity constraint through equivalent transformation. The processing process of transformations is as follows.Since the probability distribution function, expectation, and variance of each customer are known in advance, the total number of stochastic demands (∑i∈Inqiyik) belonging to the route serviced by vehicle k approximately obeys the normal distribution according to the central limit theorem, and its expectation μ and standard deviation δ are as follows:(13)μ=∑i∈IEqi⋅yik,δ=u,then normalizing, we have(14)∑i∈Inqiyik−nμnδ∼N0,1.Therefore, equation (2) is rewritten as(15)Pr∑i∈Inqiyik−nμnδ≤U−nμnδ≈ϕU−nμnδ≥α.That is, when the risk preferenceα is determined, we can obtain n, i.e., the number of customers belonging to the route serviced by vehicle k.Equation (3) is the vehicle load constraint; i.e., the expectation of stochastic customer’s demand is not greater than the maximum capacity. Equation (4) indicates that the number of vehicles delivering goods is less than or equal to the total number of vehicles in the fleet. Equation (5) represents that the total number of customers served by vehicles is not greater than the total number of customers. Equation (6) guarantees that the remaining battery capacity of electric vehicles stays the same between it reaching the customer i and departing. Equation (7) indicates that the electric vehicle is fully charged when it leaves the depot or the battery swap station. Equation (8) indicates that the electric vehicle has a nonnegative battery capacity anywhere. Equation (9) represents the range of remaining battery capacity of EVs from the battery swap station or the distribution center to j. Meanwhile, equation (10) guarantees that each EV has sufficient battery power to visit customer j from i. Equation (11) ensures that each EV is fully loaded when it leaves the distribution center. Finally, equation (12) represents the remaining load of the electric vehicle when it departs from customer j.


### 4.2. Insertion Cost of Battery Swap Stations

On the basis of the study by Lin et al. [19], we considered the congestion factor of the battery swap station and established a new method to calculate the insertion cost. When a battery swap station is selected from the candidate set B(B⊂V) and inserted into the route rk, an additional cost will be obtained with equation (18):(16)minCb,kv¯=α1−gb,kv¯+α2lb,kv¯+α3eb,kv¯+α3Lb,kv¯,b∈B,v∈rk,rk∈R,where(17)α1+α2+α3+α4=1,α1,α2,α3,α4≥0,(18)gb,kv¯=∑v∈rk\0qv,k2∧−qv,k21+qk∗1+qk∗∗,(19)lb,kv¯=dvn,v¯+dv¯,vn+1−dvn,vn+1,vn,vn+1∈rk,(20)eb,kv¯=Mqv¯,k,(21)Lb,kv¯=ftv¯1.gb,kv¯ represents the feasibility improvement value of the route rk when inserting the battery swap station b into the location v¯ after node v in the route rk; qv,k∧ and qk∗∗ denote the solution feasibility of the route rk at node v after inserting b into v¯ and the worst node feasibility state in the route rk, respectively.lb,kv¯ is the cost increase caused by insertion of the battery swap station. vnandvn+1 represent the points before and after the adjacent v¯, respectively.eb,kv¯ means a penalty function and M is a maximum number. qv¯,k represents the node state of the battery swap station b at the node v¯. If the inserted position of b is a breakpoint, the battery swap station is considered unreachable, so this unfeasible insertion needs to be punished. In addition, let pv,k2 be the remaining battery capacity when the vehicle leaves point v in the route rk. For each node in the route rk, there is a set of reachable battery swap stations, denoted as Bvk⊆B, and thus Bvk=b∈Bpv,k2≥dv,bv∈rk. ∀b in Bvk, if it is selected as the battery swap station after insertion into node v, then eb,kv¯ = 0.Lb,kv¯ denotes the waiting cost generated in the process of battery replacement, which is linearly correlated with the arrival time at the battery swap station. For each station, there is a corresponding congestion function. Two types of battery swap stations are defined according to the geographical location. The first-level station represents the stations with a high probability of being visited in the city, and the second-level station represents the stations with a low probability of being visited. The congestion function of the two types is shown in Figure 2.Figure 2
Congestion functions for the two types of stations.


### 4.3. Insertion Cost of Dynamic Customers

For dynamic customer setw∈W, when inserting the dynamic customer w into the route rk, an additional cost Cw,km will occur, where m is the location of the dynamic customer, m∈M, and then we calculate the insertion cost of the dynamic customer as follows:(22)minCw,km=β1lw,km+β2ew,km+β3Cb,kv¯,w∈W,m∈rk,rk⊂R,where(23)β1+β2+β3=1,β1,β2,β≥0,(24)lw,km=dvj,v+dv,vj+1−dvj,vj+1,vj,vj+1∈rk,(25)ei,km=Si,km−SkI,Iifupk2−ui<0,in which lw,km represents the increase in distance caused by inserting the dynamic customer w into path rk in equation (24); ei,km represents a penalty, and two conditions will trigger the penalty mechanism: (a) when the dynamic customer i is too far from the current insertion position; that is, the trips to the battery swap stations increase more than twice after the dynamic customer i is inserted; (b) the vehicle overload when a dynamic customer is added to the route. And, in equation (25), I is a maximum number, Sk represents the number of trips to the battery swap stations of the original route rk, Si,km represents the number of trips to the battery swap stations when the dynamic customer i is inserted into route rk at m, upk2 represents the amount of load left when the vehicle leaves the key point or the remedial point, and ui represents the demand of the dynamic customer i. Finally, Cb,kv¯ represents the insertion cost of the battery swap station.


## 4.1.EVRPSD Model

In order to better introduce the model, the symbol description is given in Table1.Table 1
Notation used in this paper.
NotationMeaningIThe set of stochastic customers with indexiWThe set of dynamic customers with indexwBThe set of battery swap stations with indexbDThe set of distribution centers with indexoVThe set of vertexesKThe set of electric vehicles with indexkdijThe distance fromi to jqiThe actual demand of customerigwThe dynamic demand of customerwEqiThe expectation of customeri’s demandrPower consumption per kilometerαThe decision maker’s attitude towards riskpvk1Remaining battery capacity of vehiclek on arrival at vpvk2Remaining battery capacity of vehiclek on departing at vujk2The remaining load when the vehicle leaves pointjUVehicle capacityQBattery capacityC1Unit distance cost (¥/km) of electric vehiclesC2Fixed cost for each electric vehicletijTravel time (hours) on arc (i, j)speedTravel speed (km/h) of electric vehiclesxijk=1,0xijk=1 if the vehicle k flows on arc (i, j); otherwise, xijk=0yik=1,0yik=1 if the customer i is serviced by vehicle k; otherwise, yik=0xok=1,0xok=1 if the vehicle k returns to the depot after completing its tasks; otherwise, xok=0The mixed integer linear programming is established as follows:(1)minC1∑k∈K∑i∈V∑j∈V,i≠jxijkdij+C2∑k∈Kxok,s.t.∑i∈V,i≠jxijk−∑i∈V,i≠jxjik=0,∀j∈V,∀k∈K,(2)Pr∑i∈I∑j∈I,j≠ixijk⋅qj≤U≥α,α∈0,1,∀k∈K,(3)∑i∈IykiEqi≤U,∀k∈K,(4)∑k∈K∑i∈V,i≠oxoik≤k,∀o∈D,(5)∑i∈I∑j∈I,j≠ixijk≤I,∀k∈K,(6)pik1=pik2,∀i∈I,∀k∈K,(7)pik2=Q,∀i∈D∪B,∀k∈K,(8)pvk1≥0,∀v∈V,∀k∈K,(9)0≤pjk1≤Q−rdijxijk,μ,(10)0≤pjk1≤pik1−rdijxijk+Q1−xijk,i∈I,j∈V,j≠i,∀k∈K,(11)uok2=U,k∈K,o∈D,(12)0≤ujk2=ujk1−qj,j∈I,k∈K.Equation (1) represents the flow balance, which ensures that vehicles enter and exit from the same point. Equation (2) is a stochastic opportunity constraint, that is, only go to point j when the probability that the current remaining load can meet the demand at point j is greater than or equal to the risk preference α. This constraint can be transformed into a deterministic opportunity constraint through equivalent transformation. The processing process of transformations is as follows.Since the probability distribution function, expectation, and variance of each customer are known in advance, the total number of stochastic demands (∑i∈Inqiyik) belonging to the route serviced by vehicle k approximately obeys the normal distribution according to the central limit theorem, and its expectation μ and standard deviation δ are as follows:(13)μ=∑i∈IEqi⋅yik,δ=u,then normalizing, we have(14)∑i∈Inqiyik−nμnδ∼N0,1.Therefore, equation (2) is rewritten as(15)Pr∑i∈Inqiyik−nμnδ≤U−nμnδ≈ϕU−nμnδ≥α.That is, when the risk preferenceα is determined, we can obtain n, i.e., the number of customers belonging to the route serviced by vehicle k.Equation (3) is the vehicle load constraint; i.e., the expectation of stochastic customer’s demand is not greater than the maximum capacity. Equation (4) indicates that the number of vehicles delivering goods is less than or equal to the total number of vehicles in the fleet. Equation (5) represents that the total number of customers served by vehicles is not greater than the total number of customers. Equation (6) guarantees that the remaining battery capacity of electric vehicles stays the same between it reaching the customer i and departing. Equation (7) indicates that the electric vehicle is fully charged when it leaves the depot or the battery swap station. Equation (8) indicates that the electric vehicle has a nonnegative battery capacity anywhere. Equation (9) represents the range of remaining battery capacity of EVs from the battery swap station or the distribution center to j. Meanwhile, equation (10) guarantees that each EV has sufficient battery power to visit customer j from i. Equation (11) ensures that each EV is fully loaded when it leaves the distribution center. Finally, equation (12) represents the remaining load of the electric vehicle when it departs from customer j.


## 4.2. Insertion Cost of Battery Swap Stations

On the basis of the study by Lin et al. [19], we considered the congestion factor of the battery swap station and established a new method to calculate the insertion cost. When a battery swap station is selected from the candidate set B(B⊂V) and inserted into the route rk, an additional cost will be obtained with equation (18):(16)minCb,kv¯=α1−gb,kv¯+α2lb,kv¯+α3eb,kv¯+α3Lb,kv¯,b∈B,v∈rk,rk∈R,where(17)α1+α2+α3+α4=1,α1,α2,α3,α4≥0,(18)gb,kv¯=∑v∈rk\0qv,k2∧−qv,k21+qk∗1+qk∗∗,(19)lb,kv¯=dvn,v¯+dv¯,vn+1−dvn,vn+1,vn,vn+1∈rk,(20)eb,kv¯=Mqv¯,k,(21)Lb,kv¯=ftv¯1.gb,kv¯ represents the feasibility improvement value of the route rk when inserting the battery swap station b into the location v¯ after node v in the route rk; qv,k∧ and qk∗∗ denote the solution feasibility of the route rk at node v after inserting b into v¯ and the worst node feasibility state in the route rk, respectively.lb,kv¯ is the cost increase caused by insertion of the battery swap station. vnandvn+1 represent the points before and after the adjacent v¯, respectively.eb,kv¯ means a penalty function and M is a maximum number. qv¯,k represents the node state of the battery swap station b at the node v¯. If the inserted position of b is a breakpoint, the battery swap station is considered unreachable, so this unfeasible insertion needs to be punished. In addition, let pv,k2 be the remaining battery capacity when the vehicle leaves point v in the route rk. For each node in the route rk, there is a set of reachable battery swap stations, denoted as Bvk⊆B, and thus Bvk=b∈Bpv,k2≥dv,bv∈rk. ∀b in Bvk, if it is selected as the battery swap station after insertion into node v, then eb,kv¯ = 0.Lb,kv¯ denotes the waiting cost generated in the process of battery replacement, which is linearly correlated with the arrival time at the battery swap station. For each station, there is a corresponding congestion function. Two types of battery swap stations are defined according to the geographical location. The first-level station represents the stations with a high probability of being visited in the city, and the second-level station represents the stations with a low probability of being visited. The congestion function of the two types is shown in Figure 2.Figure 2
Congestion functions for the two types of stations.


## 4.3. Insertion Cost of Dynamic Customers

For dynamic customer setw∈W, when inserting the dynamic customer w into the route rk, an additional cost Cw,km will occur, where m is the location of the dynamic customer, m∈M, and then we calculate the insertion cost of the dynamic customer as follows:(22)minCw,km=β1lw,km+β2ew,km+β3Cb,kv¯,w∈W,m∈rk,rk⊂R,where(23)β1+β2+β3=1,β1,β2,β≥0,(24)lw,km=dvj,v+dv,vj+1−dvj,vj+1,vj,vj+1∈rk,(25)ei,km=Si,km−SkI,Iifupk2−ui<0,in which lw,km represents the increase in distance caused by inserting the dynamic customer w into path rk in equation (24); ei,km represents a penalty, and two conditions will trigger the penalty mechanism: (a) when the dynamic customer i is too far from the current insertion position; that is, the trips to the battery swap stations increase more than twice after the dynamic customer i is inserted; (b) the vehicle overload when a dynamic customer is added to the route. And, in equation (25), I is a maximum number, Sk represents the number of trips to the battery swap stations of the original route rk, Si,km represents the number of trips to the battery swap stations when the dynamic customer i is inserted into route rk at m, upk2 represents the amount of load left when the vehicle leaves the key point or the remedial point, and ui represents the demand of the dynamic customer i. Finally, Cb,kv¯ represents the insertion cost of the battery swap station.


## 5. TheCW-TS Heuristic

To obtain the initial routes of the model, a staticEVRPSD is solved in the first phase of the first cycle of the problem. We design a hybrid heuristic algorithm composed of the Clarke and Wright’s saving method and an improved Tabu search algorithm (CW-TS), in which the saving method algorithm obtains an initial solution and then the improved Tabu search algorithm further improves the solution quality. In the second phase of the first cycle, dynamic customers are inserted into the routes of the initial solution by means of the greedy algorithm to complete the routes’ update, and the subsequent cycles are solved by reusing the algorithm of the first cycle. The research method and the research processes are shown in Figure 3.Figure 3
Algorithm framework.


### 5.1. The Saving Method for Generating an Initial Solution

#### 5.1.1. Coding

In this paper, natural numbers are used to encode the distribution center, customers, and battery swap stations. For example, a case contains 10 stochastic customers, 3 dynamic customers, and 2 battery swap stations, and if we use 1–10 representing stochastic customers, 11-12 representing battery swap stations, 13–15 representing dynamic customers, and 0 representing the distribution center, then [0-1-3-11-5-0-1-3-11-6-0-7-8-9-0] means a feasible solution of the case’s first phase and [0-1-3-11-5-13-0-2-4-12-14-6-0-2-4-12-15-0] indicates the second phase of a feasible solution after inserting the dynamic customers and the battery swap stations.


#### 5.1.2. The Process of Generating the Initial Solution

Step 1. The saving method algorithm is utilized to obtain a giant TSP solution including all stochastic customers without considering the capacity constraint, demand constraint, and battery capacity constraint.Step 2. Relaxing the constraint of vehicle endurance, the giant TSP is divided into a VRP solution based on the capacity constraint of vehicles. At this point, the VRP solution between the distribution center and all stochastic customers is obtained.Step 3. The vehicle’s endurance constraint is used to check whether the solution is feasible or not. If it is feasible, reserve it. Otherwise, an insertion operator of the battery swap station is used to obtain an initial solution.


#### 5.1.3. Insertion Operator of Battery Swap Stations

After the steps above, we obtain an initial solution ofVRP. Since the battery capacity of the electric vehicle is limited, the solution obtained may not be feasible, so we need to insert battery swap stations into the cheapest place of the infeasible routes.For an infeasible routerk, we find a break point (node 4 in Figure 4) v∗, so the insertion position of the battery swap station can be the positions before v∗. Then, an exhaustive search is performed on candidate battery swap stations and insertion locations to select the optimal combination of a station and location with the lowest insertion cost according to formula (18). Once a battery swap station is inserted, it is necessary to determine whether the route is feasible. If not, the above operation is repeated until the route is feasible. Finally, an initial solution of the first phase of the cycle is obtained after completing the above steps.Figure 4
An illustration of inserting a battery swap station.


### 5.2. An Improved Tabu Search Algorithm

We use an improved Tabu search algorithm to optimize the initial solution obtained by the saving method. That is, a new VRP solution is obtained through neighborhood transformation; thus a new EVRP solution is generated, and the generation mode of the neighborhood candidates and the stopping criteria of the algorithm are as follows.


#### 5.2.1. The Generation of the Neighborhood Candidates

According to the current solution, an inverse operator, 1-opt exchange operator, 2-opt exchange operator, and 3-opt exchange operator are randomly used to generate neighborhood candidates. In each iteration,P candidates are obtained by randomly selecting the above operators, and the candidates with the minimum objective function value are added to the Tabu list, in which the length is the parameter Tlength.


#### 5.2.2. Stopping Criteria

When a current optimal solution keepsStpon times without changing or the number of iterations reaches the Max_iteration, the Tabu search stops and outputs the current solution for the first phase of the cycle.


### 5.3. A Greedy Insertion Algorithm for Dynamic Customer

We obtain initial routes to serve stochastic customers through theCW-TS algorithm, which can guide vehicles for delivery. As the actual demand of stochastic customers is not always consistent with the forecast, forecast deviation may occur in the distribution process. When the vehicle fails to meet the needs of the current node, a corresponding remedial strategy is required.There are two classic remedy strategies: (1) when the probability of the success of the next customeri + 1 is less than a probability α, the vehicle will directly return to the depot from node i; (2) when the vehicle arrives at node i and the remaining load is insufficient, the vehicle also directly returns to the depot.The classic strategies raise two questions: (a) when the case (1) occurs, the vehicle returns to the depot with some remaining, which is obviously wasted; (b) when the vehicle successfully serves the last stochastic customer on its route and it also has a surplus load, the vehicle has to return to the depot.Except stochastic customers, the model we proposed may also have dynamic customers whose demand and location can hardly be statistically predicted until their orders are confirmed. Therefore, we design a dynamic customer insertion strategy that allows the dynamic customers to “consume” the remaining load.


#### 5.3.1. The Insertion Strategy for Dynamic Customer

For the sake of simplicity, we only illustrate the method that inserts the dynamic customers into the routes with remedial points, and inserting behind key points is in a similar way. First, the algorithm searches the recovery pointpk of each route rk and records the remainder load and the battery capacity of each pk, removing the battery swap stations and unserved stochastic customers after remedial points. Then, for each dynamic customer i in the set W, we can assign it to the route with remedial points, which is a many-to-many matching problem, where a dynamic customer may be inserted into multiple routes and which location is the best in the route is decided according to equation (23). The detailed insertion process is shown in Figure 5.Figure 5
The insertion process of the dynamic customer.Finally, the unassigned dynamic customers and the unserved stochastic customers are stored in the set for the next delivery cycle.


#### 5.3.2. Matching Dynamic Customers and Routes

We use the greedy algorithm to rapidly search the nearest route of dynamic customers. The central coordinates of the routes are found by the K-means algorithm, and then the distance from each dynamic customer to these central coordinates is calculated. The route closest to the dynamic customer is the inserted candidate.


## 5.1. The Saving Method for Generating an Initial Solution

### 5.1.1. Coding

In this paper, natural numbers are used to encode the distribution center, customers, and battery swap stations. For example, a case contains 10 stochastic customers, 3 dynamic customers, and 2 battery swap stations, and if we use 1–10 representing stochastic customers, 11-12 representing battery swap stations, 13–15 representing dynamic customers, and 0 representing the distribution center, then [0-1-3-11-5-0-1-3-11-6-0-7-8-9-0] means a feasible solution of the case’s first phase and [0-1-3-11-5-13-0-2-4-12-14-6-0-2-4-12-15-0] indicates the second phase of a feasible solution after inserting the dynamic customers and the battery swap stations.


### 5.1.2. The Process of Generating the Initial Solution

Step 1. The saving method algorithm is utilized to obtain a giant TSP solution including all stochastic customers without considering the capacity constraint, demand constraint, and battery capacity constraint.Step 2. Relaxing the constraint of vehicle endurance, the giant TSP is divided into a VRP solution based on the capacity constraint of vehicles. At this point, the VRP solution between the distribution center and all stochastic customers is obtained.Step 3. The vehicle’s endurance constraint is used to check whether the solution is feasible or not. If it is feasible, reserve it. Otherwise, an insertion operator of the battery swap station is used to obtain an initial solution.


### 5.1.3. Insertion Operator of Battery Swap Stations

After the steps above, we obtain an initial solution ofVRP. Since the battery capacity of the electric vehicle is limited, the solution obtained may not be feasible, so we need to insert battery swap stations into the cheapest place of the infeasible routes.For an infeasible routerk, we find a break point (node 4 in Figure 4) v∗, so the insertion position of the battery swap station can be the positions before v∗. Then, an exhaustive search is performed on candidate battery swap stations and insertion locations to select the optimal combination of a station and location with the lowest insertion cost according to formula (18). Once a battery swap station is inserted, it is necessary to determine whether the route is feasible. If not, the above operation is repeated until the route is feasible. Finally, an initial solution of the first phase of the cycle is obtained after completing the above steps.Figure 4
An illustration of inserting a battery swap station.


## 5.1.1. Coding

In this paper, natural numbers are used to encode the distribution center, customers, and battery swap stations. For example, a case contains 10 stochastic customers, 3 dynamic customers, and 2 battery swap stations, and if we use 1–10 representing stochastic customers, 11-12 representing battery swap stations, 13–15 representing dynamic customers, and 0 representing the distribution center, then [0-1-3-11-5-0-1-3-11-6-0-7-8-9-0] means a feasible solution of the case’s first phase and [0-1-3-11-5-13-0-2-4-12-14-6-0-2-4-12-15-0] indicates the second phase of a feasible solution after inserting the dynamic customers and the battery swap stations.


## 5.1.2. The Process of Generating the Initial Solution

Step 1. The saving method algorithm is utilized to obtain a giant TSP solution including all stochastic customers without considering the capacity constraint, demand constraint, and battery capacity constraint.Step 2. Relaxing the constraint of vehicle endurance, the giant TSP is divided into a VRP solution based on the capacity constraint of vehicles. At this point, the VRP solution between the distribution center and all stochastic customers is obtained.Step 3. The vehicle’s endurance constraint is used to check whether the solution is feasible or not. If it is feasible, reserve it. Otherwise, an insertion operator of the battery swap station is used to obtain an initial solution.


## 5.1.3. Insertion Operator of Battery Swap Stations

After the steps above, we obtain an initial solution ofVRP. Since the battery capacity of the electric vehicle is limited, the solution obtained may not be feasible, so we need to insert battery swap stations into the cheapest place of the infeasible routes.For an infeasible routerk, we find a break point (node 4 in Figure 4) v∗, so the insertion position of the battery swap station can be the positions before v∗. Then, an exhaustive search is performed on candidate battery swap stations and insertion locations to select the optimal combination of a station and location with the lowest insertion cost according to formula (18). Once a battery swap station is inserted, it is necessary to determine whether the route is feasible. If not, the above operation is repeated until the route is feasible. Finally, an initial solution of the first phase of the cycle is obtained after completing the above steps.Figure 4
An illustration of inserting a battery swap station.


## 5.2. An Improved Tabu Search Algorithm

We use an improved Tabu search algorithm to optimize the initial solution obtained by the saving method. That is, a new VRP solution is obtained through neighborhood transformation; thus a new EVRP solution is generated, and the generation mode of the neighborhood candidates and the stopping criteria of the algorithm are as follows.


### 5.2.1. The Generation of the Neighborhood Candidates

According to the current solution, an inverse operator, 1-opt exchange operator, 2-opt exchange operator, and 3-opt exchange operator are randomly used to generate neighborhood candidates. In each iteration,P candidates are obtained by randomly selecting the above operators, and the candidates with the minimum objective function value are added to the Tabu list, in which the length is the parameter Tlength.


### 5.2.2. Stopping Criteria

When a current optimal solution keepsStpon times without changing or the number of iterations reaches the Max_iteration, the Tabu search stops and outputs the current solution for the first phase of the cycle.


## 5.2.1. The Generation of the Neighborhood Candidates

According to the current solution, an inverse operator, 1-opt exchange operator, 2-opt exchange operator, and 3-opt exchange operator are randomly used to generate neighborhood candidates. In each iteration,P candidates are obtained by randomly selecting the above operators, and the candidates with the minimum objective function value are added to the Tabu list, in which the length is the parameter Tlength.


## 5.2.2. Stopping Criteria

When a current optimal solution keepsStpon times without changing or the number of iterations reaches the Max_iteration, the Tabu search stops and outputs the current solution for the first phase of the cycle.


## 5.3. A Greedy Insertion Algorithm for Dynamic Customer

We obtain initial routes to serve stochastic customers through theCW-TS algorithm, which can guide vehicles for delivery. As the actual demand of stochastic customers is not always consistent with the forecast, forecast deviation may occur in the distribution process. When the vehicle fails to meet the needs of the current node, a corresponding remedial strategy is required.There are two classic remedy strategies: (1) when the probability of the success of the next customeri + 1 is less than a probability α, the vehicle will directly return to the depot from node i; (2) when the vehicle arrives at node i and the remaining load is insufficient, the vehicle also directly returns to the depot.The classic strategies raise two questions: (a) when the case (1) occurs, the vehicle returns to the depot with some remaining, which is obviously wasted; (b) when the vehicle successfully serves the last stochastic customer on its route and it also has a surplus load, the vehicle has to return to the depot.Except stochastic customers, the model we proposed may also have dynamic customers whose demand and location can hardly be statistically predicted until their orders are confirmed. Therefore, we design a dynamic customer insertion strategy that allows the dynamic customers to “consume” the remaining load.


### 5.3.1. The Insertion Strategy for Dynamic Customer

For the sake of simplicity, we only illustrate the method that inserts the dynamic customers into the routes with remedial points, and inserting behind key points is in a similar way. First, the algorithm searches the recovery pointpk of each route rk and records the remainder load and the battery capacity of each pk, removing the battery swap stations and unserved stochastic customers after remedial points. Then, for each dynamic customer i in the set W, we can assign it to the route with remedial points, which is a many-to-many matching problem, where a dynamic customer may be inserted into multiple routes and which location is the best in the route is decided according to equation (23). The detailed insertion process is shown in Figure 5.Figure 5
The insertion process of the dynamic customer.Finally, the unassigned dynamic customers and the unserved stochastic customers are stored in the set for the next delivery cycle.


### 5.3.2. Matching Dynamic Customers and Routes

We use the greedy algorithm to rapidly search the nearest route of dynamic customers. The central coordinates of the routes are found by the K-means algorithm, and then the distance from each dynamic customer to these central coordinates is calculated. The route closest to the dynamic customer is the inserted candidate.


## 5.3.1. The Insertion Strategy for Dynamic Customer

For the sake of simplicity, we only illustrate the method that inserts the dynamic customers into the routes with remedial points, and inserting behind key points is in a similar way. First, the algorithm searches the recovery pointpk of each route rk and records the remainder load and the battery capacity of each pk, removing the battery swap stations and unserved stochastic customers after remedial points. Then, for each dynamic customer i in the set W, we can assign it to the route with remedial points, which is a many-to-many matching problem, where a dynamic customer may be inserted into multiple routes and which location is the best in the route is decided according to equation (23). The detailed insertion process is shown in Figure 5.Figure 5
The insertion process of the dynamic customer.Finally, the unassigned dynamic customers and the unserved stochastic customers are stored in the set for the next delivery cycle.


## 5.3.2. Matching Dynamic Customers and Routes

We use the greedy algorithm to rapidly search the nearest route of dynamic customers. The central coordinates of the routes are found by the K-means algorithm, and then the distance from each dynamic customer to these central coordinates is calculated. The route closest to the dynamic customer is the inserted candidate.


## 6. Numerical Results

This section consists of two parts. The first part presents a medium-sized case to verify the effectiveness of the model, and the second part compares our algorithm with the classical Tabu search algorithm to verify the feasibility and accuracy of the algorithm by using 9 modified Solomon instances. In this paper, MATLAB R2019a software based on the Windows 10 operating system is adopted as the simulation tool, and the hardware platform is Intel Core i7-6550u CPU, 8 GB RAM and 256 GB SSD.


### 6.1. A Real Case Study

The real case investigated in this section comes from an express company in Beijing, including a distribution center (DC), 50 stochastic customers, 5 battery swap stations, and 10 dynamic customers. The specific parameters are shown in Table2, and the detailed information is shown in Tables 3–5. Figure 6 shows the geographical location of DC and customers on the Google map, where the yellow mark represents the DC, the blue mark numbered 1–50 represents the stochastic customer, the green mark numbered 51–55 represents the battery swap stations, and the pink mark numbered 56–65 represents the dynamic customer.Table 2
Parameters.
ParameterDescriptionValueUnitmaxDrivingRangeThe battery driving range120KilometerUThe loading capacity2.5Metric tonC2Vehicle fixed cost100¥/CarC1Variable costs10¥/kmrPower consumption per kilometer1kW·h/kmSpeedThe travel speed60km/hcandidateSolutionThe number of candidates inTS30—currentBestThe current best solution without improvement for a certain number10—tabuListLengthThe length of Tabu list10—αRisk preference0.6—tabuIterationMaximum number of iterations inTS50—Table 3
The detailed information of DC and dynamic customers.
No.LongitudeLatitudeForecasted demandActual demand(0)116.327832739.90135882——(1)116.6419940.0487960.256830.333879(2)116.38116339.7990750.25860.33618(3)116.11798540.0453860.2632350.210588(4)116.53184239.8294940.1895260.1326682(5)116.39915640.0452140.152410.198133(6)116.61931440.1183430.189540.246402(7)116.32212139.9545170.168730.219349(8)116.30277340.1058160.1576350.189162(9)116.50219340.0091990.00650.0078(10)116.68785839.9734750.215510.258612(11)116.52628639.7478750.004320.003456(12)116.41755939.9694150.00860.01032(13)116.49060239.9633390.1706350.136508(14)116.20536640.1451430.09240.11088(15)116.47312739.9986210.164050.213265(16)116.56415639.9192560.183760.220512(17)116.62811939.824310.281850.366405(18)116.23016840.0717810.019880.015904(19)116.37650840.0313230.120840.145008(20)116.29879739.952490.00980.00686(21)116.18828839.7858150.180680.234884(22)116.1249439.7894160.001440.001008(23)116.25768840.0993610.001380.000966(24)116.5561739.9369150.098350.07868(25)116.63592239.8321810.2784050.222724(26)116.10512639.743710.00290.00232(27)116.46641139.9886460.002490.001992(28)116.59830340.0509730.151830.182196(29)116.42464840.1275480.194230.233076(30)116.37691539.9719170.001470.001029(31)116.70764139.918890.03360.02352(32)116.19884139.7083630.058490.070188(33)116.49563939.9112950.201160.160928(34)116.2246439.7827930.1494550.1046185(35)116.5683639.8381980.0280.0336(36)116.35340739.9327410.10760.08608(37)116.27630239.7806740.175240.122668(38)116.64619939.8792430.03440.02408(39)116.32645439.8842050.2340.2808(40)116.47044539.789310.21460.25752(41)116.45459239.9192550.00630.00756(42)116.18917639.9126790.609310.792103(43)116.1664339.9918480.3612750.28902(44)116.31205439.8281310.1895180.1516144(45)116.5612139.972870.054320.070616(46)116.3921339.7477720.01460.01022(47)116.19659939.7152920.170240.221312(48)116.45281540.0404560.2496760.3245788(49)116.30385240.0544740.155060.201578(50)116.41385439.842250.435260.565838Table 4
The detailed information of battery swap stations.
BSSLongitudeLatitudeLevel(51)116.15057339.702081(52)116.52635939.7683521(53)116.15096939.9376662(54)116.4713340.0842542(55)116.39488939.891922Table 5
The detailed information of dynamic customers.
No.LongitudeLatitudeDemand(56)116.1109139.9655850.7481(57)116.2752440.1316810.4505(58)116.5366340.0016290.0838(59)116.3826139.8501520.2289(60)116.2844340.0515820.6133(61)116.3520740.018010.1523(62)116.5890840.0474720.5258(63)116.5867439.9460240.5383(64)116.3894339.976760.4961(65)116.2206239.7745720.0781Figure 6
The geographic location of the nodes.TheCW-TS algorithm is carried out 10 times, and the results are shown in Table 6. The total distance in the first phase of the first cycle is reported by column L1, the total distance in the second phase of the first cycle is reported by column L2, the total distance in the second cycle is reported by L3, and the column K, reports the number of vehicles used in this case. Finally, the total cost is presented in the last column.Table 6
The test result of the case.
L1L2L3KC(1)389.44444.76141.2866460.43(2)362.67473.73152.3766861.04(3)384.80418.75142.2566210.01(4)384.76432.69140.6966333.79(5)386.04450.38129.2766396.50(6)389.41479.08129.2766683.48(7)388.95465.88130.2666561.44(8)384.80418.75142.2566210.01(9)386.04450.38129.2766396.50(10)364.13446.60143.2266498.27The routes of the best solution are shown in Table7. C1-1 reports the routes of the first phase of the first cycle, C1-2 represents the routes of the second phase of the first cycle, and C2 indicates the routes of the second cycle.Table 7
Tour plan.
CycleRouteC1-1Route 10-2-40-4-52-11-46-37-34-21-47-32-26-22-42-20-7-36-0Route 20-27-15-9-28-6-1-10-31-38-25-17-35-16-24-45-13-33-0Route 30-43-3-14-23-18-49-8-29-54-48-5-19-30-12-41-50-0Route 40-44-39-0C1-2Route 10-2-40-4-52-11-46-37-34-21-47-32-26-22-42-20-65-0Route 20-27-15-9-28-6-1-10-31-38-25-17-35-16-24-45-0Route 30-43-3-14-23-18-49-8-29-54-48-5-19-30-12-41-64-0Route 40-44-39-59-0C2Route 10-36-7-61-60-57-56-0Route 20-13-58-62-63-33-50-0Combined with the real road network in Beijing, the actual driving paths generated by the solution are shown in Figure7.Figure 7
Actual driving path.The routes 1–6 are marked with 6 different colors, respectively, i.e., red, fuchsia, green, black, purple, and dark green.


### 6.2. Comparison of Modified Benchmark Instances

Since there are no classical benchmark instances fitting the model proposed in this paper, we modified 9 Solomon benchmark instances on different size to verify the feasibility and effectiveness of the model. The modification rules are as follows: (1) all customers of the original instances are regarded as stochastic customers, and the original demands are considered as the predicted demand of each stochastic customer. Each stochastic customer’s real demand is fluctuated by about 30% on the corresponding predicted demand, and the coordinates of the nodes remain unchanged; (2) for instances on different scales, the number of dynamic customers, the demand of dynamic customers, and the number of battery swap stations are generated randomly; (3) the vehicle capacity is 200, and the battery driving range of the vehicle is 140 km; other parameters are kept unchanged as in Table2.Then, we compare our algorithm with the classical Tabu search (TS) performed on the modified instances mentioned above. Each instance is run 10 times, and the results are shown in Table 8, where column I represents the number of stochastic customers, column W represents the number of dynamic customers, and column B means the number of battery swap stations. Finally, the column BestC reports the optimal cost, and the column Gap shows the gap between the two algorithms.Table 8
Results of comparison betweenTS and CW-TS.
TSCW-TSGap (%)InstanceIWBBestC1Time (s)BestC2Time (s)C101_2525323316.213.443231.933.145.53R101_2525325341.556.615025.436.495.92RC101_2525324896.724.914562.005.946.84C101_5050536433.9013.326040.1510.806.12R101_5050537963.5118.547426.8516.676.74RC101_5050539344.8819.108855.8318.205.23C101_100100181315344.9266.1413343.3565.6813.04R101_10010012614491.1078.3212650.1982.3512.70RC101_10010012517601.5076.1115249.0866.8013.37As can be seen from Table8, the CW-TS algorithm is superior to the classical TS algorithm, and the gap between the two solutions gradually enlarges with the increasing number of nodes. The time consumption of the two algorithms is approximate. The detailed information of the last three medium-scale instances is shown in Tables 9–11, and the tour plan of the corresponding optimal solution is shown in Tables 12–14.Table 9
Information of dynamic customers and battery swap stations in C101_100.
No.CoordinateDemandBSSCoordinate(1)[40, 40]181[40, 33](2)[75, 65]162[31, 80](3)[40, 80]183[18, 58](4)[9, 68]114[15, 30](5)[22, 40]125[45, 20](6)[22, 45]176[60, 45](7)[28, 65]127[80, 48](8)[29, 68]108[55, 71](9)[49, 10]149[70, 70](10)[48, 13]1910[20, 40](11)[70, 33]1211[25, 20](12)[72, 18]1512[40, 60](13)[80, 25]1313[85, 20](14)[80, 29]10(15)[90, 48]15(16)[20, 20]18(17)[88, 60]12(18)[85, 70]10Table 10
Information of dynamic customers and battery swap stations in R101_100.
No.CoordinateDemandBSSCoordinate(1)[10, 10]181[15, 40](2)[15, 31]22[20, 13](3)[30, 21]123[50, 25](4)[15, 12]124[35, 60](5)[38, 41]145[13, 70](6)[50, 37]186[50, 65](7)[60, 45]2(8)[60, 33]16(9)[52, 20]10(10)[55, 60]14(11)[13, 70]4(12)[35, 10]10Table 11
Information of dynamic customers and battery swap stations in RC101_100.
No.CoordinateDemandBSSCoordinate(1)[30, 80]181[20, 70](2)[31, 78]122[20, 40](3)[20, 60]163[25, 15](4)[22, 61]184[45, 75](5)[22, 58]145[80, 50](6)[15, 35]18(7)[16, 41]2(8)[23, 30]16(9)[50, 20]10(10)[53, 25]14(11)[70, 40]4(12)[65, 70]10Table 12
C101_100 optimal solution.
CycleRouteC1-1Route 10-75-1-2-4-6-9-11-10-8-7-3-5-0Route 20-13-17-18-19-15-16-14-12-0Route 30-99-100-97-93-92-94-95-96-98-0Route 40-90-87-86-83-82-84-85-88-89-91-0Route 50-66-61-72-80-79-77-73-70-71-76-78-81-0Route 60-67-65-63-74-62-64-55-68-69-0Route 70-59-60-58-56-53-54-57-0Route 80-43-42-41-40-44-46-45-48-50-51-31-52-49-47-0Route 90-29-34-36-39-38-37-35-33-32-0Route 100-21-22-23-26-28-30-27-25-24-20- 0C1-2Route 10-75-1-2-4-6-9-11-10-8-7-3-5-120-0Route 20-13-17-18-19-15-16-14-12-117-0Route 30-99-100-97-93-92-94-95-96-0Route 40-90-87-86-83-82-84-85-88-89-91-115-0Route 50-66-61-72-80-79-77-73-70-71-76-78-81-127-124-0Route 60-67-65-63-74-62-64-55-68-69-114-0Route 70-59-60-58-56-53-54-123-0Route 80-43-42-41-40-44-46-45-48-50-51-31-52-49-47-118-0Route 90-29-4-36-39-38-37-35-33-32-0Route 100-21-22-3-26-28-30-27-25-24-20-119-0C2Route 10-121-116-108-98-109-131-130-128-126-125-122-57-129-0Table 13
R101_100 optimal solution.
CycleRouteC1-1Route 10-53-58-40-21-73-74-72-4-54-12-26-28-0Route 20-80-68-3-77-76-50-51-30-70-10-62-0Route 30-5-84-17-45-46-36-49-47-48-82-8-83-60-0Route 40-19-11-105-64-63-90-32-20-66-65-71-35-34-78-9-0Route 50-33-81-79-29-24-55-25-39-67-23-56-75-0Route 60-22-41-15-43-14-44-38-86-16-61-91-100-98-0Route 70-99-59-93-85-37-92-97-87-42-57-2-13-0Route 80-94-95-96-6-89-18-52-7-88-31-1-69-27-0C1-2Route 10-53-58-40-21-73-74-72-4-54-12-26-28-112-0Route 20-80-68-3-77-76-50-51-30-70-10-62-0Route 30-5-84-17-45-46-36-49-47-48-82-117-105-0Route 40-19-11-105-64-63-90-32-20-66-65-71-35-34-78-9-116-0Route 50-33-81-79-29-24-55-25-39-67-23-56-0Route 60-22-41-15-43-14-44-38-86-16-61-91-100-98-102-110-108-0Route 70-99-59-93-85-37-92-97-87-42-57-2-0Route 80-94-95-96-6-89-18-52-7-88-31-1-69-27-0C2Route 10-111-113-114-103-115-75-118-109-13-60-83-107-8-0Table 14
RC101_100 optimal solution.
CycleRouteC1-1Route 10-68-61-70-100-55-88-53-10-11-9-0Route 20-13-15-16-17-47-14-12-78-73-79-60-0Route 30-1-3-5-45-4-46-8-7-6-2-0Route 40-42-44-43-40-36-35-37-38-39-0Route 50-41-72-71-105-32-33-89-76-63-85-51-84-56-0Route 60-93-67-34-31-29-27-26-28-30-50-62-95-0Route 70-65-66-83-24-64-91-92-94-96-54-81-80-0Route 80-20-49-19-18-48-21-23-25-22-0Route 90-57-103-77-58-75-97-59-87-74-86-52-99-0Route 100-69-98-82-90-0C1-2Route 10-68-61-70-100-55-88-53-10-11- 9-111-112-0Route 20-13-15-16-17-47-14-12-78-73-79-60-110-0Route 30-1-3-5-45-4-46-8-7-6-0Route 40-42-44-43-40-36-35-37-38-39-117-0Route 50-41-72-71-105-32-33-89-76-63-85-51-84-56-103-114-0Route 60-93-67-34-31-29-27-26-28-30-50-62-95-116-105-0Route 70-65-66-83-24-64-91-92-94-96-54-81-0Route 80-20-49-19-18-48-21-23-25-22-0Route 90-57-103-77-58-75-97-59-87-74-86-52-0Route 100-69-98-82-90-0C2Route 10-107-106-2-101-109-108-99-113-115-80-0


## 6.1. A Real Case Study

The real case investigated in this section comes from an express company in Beijing, including a distribution center (DC), 50 stochastic customers, 5 battery swap stations, and 10 dynamic customers. The specific parameters are shown in Table2, and the detailed information is shown in Tables 3–5. Figure 6 shows the geographical location of DC and customers on the Google map, where the yellow mark represents the DC, the blue mark numbered 1–50 represents the stochastic customer, the green mark numbered 51–55 represents the battery swap stations, and the pink mark numbered 56–65 represents the dynamic customer.Table 2
Parameters.
ParameterDescriptionValueUnitmaxDrivingRangeThe battery driving range120KilometerUThe loading capacity2.5Metric tonC2Vehicle fixed cost100¥/CarC1Variable costs10¥/kmrPower consumption per kilometer1kW·h/kmSpeedThe travel speed60km/hcandidateSolutionThe number of candidates inTS30—currentBestThe current best solution without improvement for a certain number10—tabuListLengthThe length of Tabu list10—αRisk preference0.6—tabuIterationMaximum number of iterations inTS50—Table 3
The detailed information of DC and dynamic customers.
No.LongitudeLatitudeForecasted demandActual demand(0)116.327832739.90135882——(1)116.6419940.0487960.256830.333879(2)116.38116339.7990750.25860.33618(3)116.11798540.0453860.2632350.210588(4)116.53184239.8294940.1895260.1326682(5)116.39915640.0452140.152410.198133(6)116.61931440.1183430.189540.246402(7)116.32212139.9545170.168730.219349(8)116.30277340.1058160.1576350.189162(9)116.50219340.0091990.00650.0078(10)116.68785839.9734750.215510.258612(11)116.52628639.7478750.004320.003456(12)116.41755939.9694150.00860.01032(13)116.49060239.9633390.1706350.136508(14)116.20536640.1451430.09240.11088(15)116.47312739.9986210.164050.213265(16)116.56415639.9192560.183760.220512(17)116.62811939.824310.281850.366405(18)116.23016840.0717810.019880.015904(19)116.37650840.0313230.120840.145008(20)116.29879739.952490.00980.00686(21)116.18828839.7858150.180680.234884(22)116.1249439.7894160.001440.001008(23)116.25768840.0993610.001380.000966(24)116.5561739.9369150.098350.07868(25)116.63592239.8321810.2784050.222724(26)116.10512639.743710.00290.00232(27)116.46641139.9886460.002490.001992(28)116.59830340.0509730.151830.182196(29)116.42464840.1275480.194230.233076(30)116.37691539.9719170.001470.001029(31)116.70764139.918890.03360.02352(32)116.19884139.7083630.058490.070188(33)116.49563939.9112950.201160.160928(34)116.2246439.7827930.1494550.1046185(35)116.5683639.8381980.0280.0336(36)116.35340739.9327410.10760.08608(37)116.27630239.7806740.175240.122668(38)116.64619939.8792430.03440.02408(39)116.32645439.8842050.2340.2808(40)116.47044539.789310.21460.25752(41)116.45459239.9192550.00630.00756(42)116.18917639.9126790.609310.792103(43)116.1664339.9918480.3612750.28902(44)116.31205439.8281310.1895180.1516144(45)116.5612139.972870.054320.070616(46)116.3921339.7477720.01460.01022(47)116.19659939.7152920.170240.221312(48)116.45281540.0404560.2496760.3245788(49)116.30385240.0544740.155060.201578(50)116.41385439.842250.435260.565838Table 4
The detailed information of battery swap stations.
BSSLongitudeLatitudeLevel(51)116.15057339.702081(52)116.52635939.7683521(53)116.15096939.9376662(54)116.4713340.0842542(55)116.39488939.891922Table 5
The detailed information of dynamic customers.
No.LongitudeLatitudeDemand(56)116.1109139.9655850.7481(57)116.2752440.1316810.4505(58)116.5366340.0016290.0838(59)116.3826139.8501520.2289(60)116.2844340.0515820.6133(61)116.3520740.018010.1523(62)116.5890840.0474720.5258(63)116.5867439.9460240.5383(64)116.3894339.976760.4961(65)116.2206239.7745720.0781Figure 6
The geographic location of the nodes.TheCW-TS algorithm is carried out 10 times, and the results are shown in Table 6. The total distance in the first phase of the first cycle is reported by column L1, the total distance in the second phase of the first cycle is reported by column L2, the total distance in the second cycle is reported by L3, and the column K, reports the number of vehicles used in this case. Finally, the total cost is presented in the last column.Table 6
The test result of the case.
L1L2L3KC(1)389.44444.76141.2866460.43(2)362.67473.73152.3766861.04(3)384.80418.75142.2566210.01(4)384.76432.69140.6966333.79(5)386.04450.38129.2766396.50(6)389.41479.08129.2766683.48(7)388.95465.88130.2666561.44(8)384.80418.75142.2566210.01(9)386.04450.38129.2766396.50(10)364.13446.60143.2266498.27The routes of the best solution are shown in Table7. C1-1 reports the routes of the first phase of the first cycle, C1-2 represents the routes of the second phase of the first cycle, and C2 indicates the routes of the second cycle.Table 7
Tour plan.
CycleRouteC1-1Route 10-2-40-4-52-11-46-37-34-21-47-32-26-22-42-20-7-36-0Route 20-27-15-9-28-6-1-10-31-38-25-17-35-16-24-45-13-33-0Route 30-43-3-14-23-18-49-8-29-54-48-5-19-30-12-41-50-0Route 40-44-39-0C1-2Route 10-2-40-4-52-11-46-37-34-21-47-32-26-22-42-20-65-0Route 20-27-15-9-28-6-1-10-31-38-25-17-35-16-24-45-0Route 30-43-3-14-23-18-49-8-29-54-48-5-19-30-12-41-64-0Route 40-44-39-59-0C2Route 10-36-7-61-60-57-56-0Route 20-13-58-62-63-33-50-0Combined with the real road network in Beijing, the actual driving paths generated by the solution are shown in Figure7.Figure 7
Actual driving path.The routes 1–6 are marked with 6 different colors, respectively, i.e., red, fuchsia, green, black, purple, and dark green.


## 6.2. Comparison of Modified Benchmark Instances

Since there are no classical benchmark instances fitting the model proposed in this paper, we modified 9 Solomon benchmark instances on different size to verify the feasibility and effectiveness of the model. The modification rules are as follows: (1) all customers of the original instances are regarded as stochastic customers, and the original demands are considered as the predicted demand of each stochastic customer. Each stochastic customer’s real demand is fluctuated by about 30% on the corresponding predicted demand, and the coordinates of the nodes remain unchanged; (2) for instances on different scales, the number of dynamic customers, the demand of dynamic customers, and the number of battery swap stations are generated randomly; (3) the vehicle capacity is 200, and the battery driving range of the vehicle is 140 km; other parameters are kept unchanged as in Table2.Then, we compare our algorithm with the classical Tabu search (TS) performed on the modified instances mentioned above. Each instance is run 10 times, and the results are shown in Table 8, where column I represents the number of stochastic customers, column W represents the number of dynamic customers, and column B means the number of battery swap stations. Finally, the column BestC reports the optimal cost, and the column Gap shows the gap between the two algorithms.Table 8
Results of comparison betweenTS and CW-TS.
TSCW-TSGap (%)InstanceIWBBestC1Time (s)BestC2Time (s)C101_2525323316.213.443231.933.145.53R101_2525325341.556.615025.436.495.92RC101_2525324896.724.914562.005.946.84C101_5050536433.9013.326040.1510.806.12R101_5050537963.5118.547426.8516.676.74RC101_5050539344.8819.108855.8318.205.23C101_100100181315344.9266.1413343.3565.6813.04R101_10010012614491.1078.3212650.1982.3512.70RC101_10010012517601.5076.1115249.0866.8013.37As can be seen from Table8, the CW-TS algorithm is superior to the classical TS algorithm, and the gap between the two solutions gradually enlarges with the increasing number of nodes. The time consumption of the two algorithms is approximate. The detailed information of the last three medium-scale instances is shown in Tables 9–11, and the tour plan of the corresponding optimal solution is shown in Tables 12–14.Table 9
Information of dynamic customers and battery swap stations in C101_100.
No.CoordinateDemandBSSCoordinate(1)[40, 40]181[40, 33](2)[75, 65]162[31, 80](3)[40, 80]183[18, 58](4)[9, 68]114[15, 30](5)[22, 40]125[45, 20](6)[22, 45]176[60, 45](7)[28, 65]127[80, 48](8)[29, 68]108[55, 71](9)[49, 10]149[70, 70](10)[48, 13]1910[20, 40](11)[70, 33]1211[25, 20](12)[72, 18]1512[40, 60](13)[80, 25]1313[85, 20](14)[80, 29]10(15)[90, 48]15(16)[20, 20]18(17)[88, 60]12(18)[85, 70]10Table 10
Information of dynamic customers and battery swap stations in R101_100.
No.CoordinateDemandBSSCoordinate(1)[10, 10]181[15, 40](2)[15, 31]22[20, 13](3)[30, 21]123[50, 25](4)[15, 12]124[35, 60](5)[38, 41]145[13, 70](6)[50, 37]186[50, 65](7)[60, 45]2(8)[60, 33]16(9)[52, 20]10(10)[55, 60]14(11)[13, 70]4(12)[35, 10]10Table 11
Information of dynamic customers and battery swap stations in RC101_100.
No.CoordinateDemandBSSCoordinate(1)[30, 80]181[20, 70](2)[31, 78]122[20, 40](3)[20, 60]163[25, 15](4)[22, 61]184[45, 75](5)[22, 58]145[80, 50](6)[15, 35]18(7)[16, 41]2(8)[23, 30]16(9)[50, 20]10(10)[53, 25]14(11)[70, 40]4(12)[65, 70]10Table 12
C101_100 optimal solution.
CycleRouteC1-1Route 10-75-1-2-4-6-9-11-10-8-7-3-5-0Route 20-13-17-18-19-15-16-14-12-0Route 30-99-100-97-93-92-94-95-96-98-0Route 40-90-87-86-83-82-84-85-88-89-91-0Route 50-66-61-72-80-79-77-73-70-71-76-78-81-0Route 60-67-65-63-74-62-64-55-68-69-0Route 70-59-60-58-56-53-54-57-0Route 80-43-42-41-40-44-46-45-48-50-51-31-52-49-47-0Route 90-29-34-36-39-38-37-35-33-32-0Route 100-21-22-23-26-28-30-27-25-24-20- 0C1-2Route 10-75-1-2-4-6-9-11-10-8-7-3-5-120-0Route 20-13-17-18-19-15-16-14-12-117-0Route 30-99-100-97-93-92-94-95-96-0Route 40-90-87-86-83-82-84-85-88-89-91-115-0Route 50-66-61-72-80-79-77-73-70-71-76-78-81-127-124-0Route 60-67-65-63-74-62-64-55-68-69-114-0Route 70-59-60-58-56-53-54-123-0Route 80-43-42-41-40-44-46-45-48-50-51-31-52-49-47-118-0Route 90-29-4-36-39-38-37-35-33-32-0Route 100-21-22-3-26-28-30-27-25-24-20-119-0C2Route 10-121-116-108-98-109-131-130-128-126-125-122-57-129-0Table 13
R101_100 optimal solution.
CycleRouteC1-1Route 10-53-58-40-21-73-74-72-4-54-12-26-28-0Route 20-80-68-3-77-76-50-51-30-70-10-62-0Route 30-5-84-17-45-46-36-49-47-48-82-8-83-60-0Route 40-19-11-105-64-63-90-32-20-66-65-71-35-34-78-9-0Route 50-33-81-79-29-24-55-25-39-67-23-56-75-0Route 60-22-41-15-43-14-44-38-86-16-61-91-100-98-0Route 70-99-59-93-85-37-92-97-87-42-57-2-13-0Route 80-94-95-96-6-89-18-52-7-88-31-1-69-27-0C1-2Route 10-53-58-40-21-73-74-72-4-54-12-26-28-112-0Route 20-80-68-3-77-76-50-51-30-70-10-62-0Route 30-5-84-17-45-46-36-49-47-48-82-117-105-0Route 40-19-11-105-64-63-90-32-20-66-65-71-35-34-78-9-116-0Route 50-33-81-79-29-24-55-25-39-67-23-56-0Route 60-22-41-15-43-14-44-38-86-16-61-91-100-98-102-110-108-0Route 70-99-59-93-85-37-92-97-87-42-57-2-0Route 80-94-95-96-6-89-18-52-7-88-31-1-69-27-0C2Route 10-111-113-114-103-115-75-118-109-13-60-83-107-8-0Table 14
RC101_100 optimal solution.
CycleRouteC1-1Route 10-68-61-70-100-55-88-53-10-11-9-0Route 20-13-15-16-17-47-14-12-78-73-79-60-0Route 30-1-3-5-45-4-46-8-7-6-2-0Route 40-42-44-43-40-36-35-37-38-39-0Route 50-41-72-71-105-32-33-89-76-63-85-51-84-56-0Route 60-93-67-34-31-29-27-26-28-30-50-62-95-0Route 70-65-66-83-24-64-91-92-94-96-54-81-80-0Route 80-20-49-19-18-48-21-23-25-22-0Route 90-57-103-77-58-75-97-59-87-74-86-52-99-0Route 100-69-98-82-90-0C1-2Route 10-68-61-70-100-55-88-53-10-11- 9-111-112-0Route 20-13-15-16-17-47-14-12-78-73-79-60-110-0Route 30-1-3-5-45-4-46-8-7-6-0Route 40-42-44-43-40-36-35-37-38-39-117-0Route 50-41-72-71-105-32-33-89-76-63-85-51-84-56-103-114-0Route 60-93-67-34-31-29-27-26-28-30-50-62-95-116-105-0Route 70-65-66-83-24-64-91-92-94-96-54-81-0Route 80-20-49-19-18-48-21-23-25-22-0Route 90-57-103-77-58-75-97-59-87-74-86-52-0Route 100-69-98-82-90-0C2Route 10-107-106-2-101-109-108-99-113-115-80-0


## 7. Conclusions

Aiming at minimizing the total cost of electric vehicle routing problems with stochastic demands and dynamic remedial measures, an opportunity-constrained EVRPSD model and two insert cost models are established, and theCW-TS heuristic algorithm and a greedy algorithm are designed to solve them. Combined with a case of a logistic company in Beijing, the effectiveness of the algorithm is verified. At the same time, the accuracy and feasibility of the algorithm are verified by comparing the results of the 9 modified Solomon instances. In a theoretical sense, first of all, we studied a new variant of EVRP that considered both the randomness and dynamics of the demand; in addition, an EVRPSD model is formulated with the opportunity constraint; and finally, a two-stage heuristic algorithm is designed to solve the EVRPSD model. In a practical sense, first of all, when choosing a battery swap station, not only its location attributes but also its congestion attributes are considered; secondly, an effective remedy strategy is established to fully utilize the residual load; finally, the display of the vehicle route is enriched, thus improving the decision maker’s experience.For future research, the time windows constraint of customers can be added into the model to improve their satisfaction. In addition, the discharging process of the battery is nonlinear, so how to describe the process more accurately is also a very challenging direction.



---

*Source: 8795284-2020-08-25.xml*