Mathematical Problems in Engineering
Volume 2020 (2020), Article ID 8795284, 15 pages
https://doi.org/10.1155/2020/8795284

Research Article

Electric Vehicle Routing Problems with Stochastic Demands and Dynamic Remedial Measures

Xianlong Ge,1,2 Ziqiang Zhu,¹ and Yuanzhi Jin1,3

¹School of Economics and Management, Chongqing Jiaotong University, Chongqing 400074, China
²Key Laboratory of Intelligent Logistics Network, Chongqing Jiaotong University, Chongqing 400074, China
³Department of Computer Technology and Information Engineering, Sanmenxia Polytechnic, Sanmenxia 472000, China

Correspondence should be addressed to Xianlong Ge; gexianlong@cqjtu.edu.cn

Received 29 January 2020; Accepted 23 June 2020; Published 25 August 2020

Academic Editor: Kalyana C. Veluvolu

Copyright © 2020 Xianlong Ge et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 the capacitated 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 Section 2, 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 on CVRPSD. 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 the CVRPSD 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 on DDVRP 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 on EVRP 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 an EVRPSD 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 the EVRPSD 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. Figure 1 describes the whole process of distribution, in which stochastic customers 2, 4, and 5 and dynamic customer 3 are included in the second cycle.

(a)

(b)

(c)

Figure 1: The illustration of the whole distribution process. (a) First cycle phase 1. (b) First cycle phase 2. (c) Second cycle.

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 traditional EVRP 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 Table 1.

Table 1: Notation used in this paper.

Notation
Meaning

The set of stochastic customers with index i

The set of dynamic customers with index

The set of battery swap stations with index b

The set of distribution centers with index o

The set of vertexes

The set of electric vehicles with index k

The distance from i to j

The actual demand of customer i

The dynamic demand of customer

The expectation of customer i’s demand

Power consumption per kilometer

The decision maker’s attitude towards risk

Remaining battery capacity of vehicle k on arrival at

Remaining battery capacity of vehicle k on departing at

The remaining load when the vehicle leaves point j

Vehicle capacity

Battery capacity

Unit distance cost (¥/km) of electric vehicles

Fixed cost for each electric vehicle

Travel time (hours) on arc (i, j)

Travel speed (km/h) of electric vehicles

if the vehicle k flows on arc (i, j); otherwise,

if the customer i is serviced by vehicle k; otherwise,

if the vehicle k returns to the depot after completing its tasks; otherwise,

The mixed integer linear programming is established as follows:

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 () 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:

then normalizing, we have

Therefore, equation (2) is rewritten as

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

where

represents the feasibility improvement value of the route when inserting the battery swap station b into the location after node in the route ; and denote the solution feasibility of the route at node after inserting b into and the worst node feasibility state in the route , respectively.

is the cost increase caused by insertion of the battery swap station. represent the points before and after the adjacent , respectively.

means a penalty function and M is a maximum number. represents the node state of the battery swap station b at the node . 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 be the remaining battery capacity when the vehicle leaves point in the route . For each node in the route , there is a set of reachable battery swap stations, denoted as , and thus . in , if it is selected as the battery swap station after insertion into node , then  = 0.

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 set , when inserting the dynamic customer into the route rk, an additional cost will occur, where m is the location of the dynamic customer, , and then we calculate the insertion cost of the dynamic customer as follows:

where

in which represents the increase in distance caused by inserting the dynamic customer into path rk in equation (24); 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, represents the number of trips to the battery swap stations of the original route rk, represents the number of trips to the battery swap stations when the dynamic customer i is inserted into route rk at m, represents the amount of load left when the vehicle leaves the key point or the remedial point, and represents the demand of the dynamic customer i. Finally, represents the insertion cost of the battery swap station.

5. The CW-TS Heuristic

To obtain the initial routes of the model, a static EVRPSD 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 of VRP. 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 route rk, we find a break point (node 4 in Figure 4) , so the insertion position of the battery swap station can be the positions before . 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 keeps Stpon 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 the CW-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 customer i + 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 point pk 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 Table 2, 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.

Parameter
Description
Value
Unit

maxDrivingRange
The battery driving range
120
Kilometer

U
The loading capacity
2.5
Metric ton

C₂
Vehicle fixed cost
100
¥/Car

C₁
Variable costs
10
¥/km

r
Power consumption per kilometer
1
kW·h/km

Speed
The travel speed
60
km/h

candidateSolution
The number of candidates in TS
30
—

currentBest
The current best solution without improvement for a certain number
10
—

tabuListLength
The length of Tabu list
10
—

α
Risk preference
0.6
—

tabuIteration
Maximum number of iterations in TS
50
—

Table 3: The detailed information of DC and dynamic customers.

No.
Longitude
Latitude
Forecasted demand
Actual demand

(0)
116.3278327
39.90135882
—
—

(1)
116.64199
40.048796
0.25683
0.333879

(2)
116.381163
39.799075
0.2586
0.33618

(3)
116.117985
40.045386
0.263235
0.210588

(4)
116.531842
39.829494
0.189526
0.1326682

(5)
116.399156
40.045214
0.15241
0.198133

(6)
116.619314
40.118343
0.18954
0.246402

(7)
116.322121
39.954517
0.16873
0.219349

(8)
116.302773
40.105816
0.157635
0.189162

(9)
116.502193
40.009199
0.0065
0.0078

(10)
116.687858
39.973475
0.21551
0.258612

(11)
116.526286
39.747875
0.00432
0.003456

(12)
116.417559
39.969415
0.0086
0.01032

(13)
116.490602
39.963339
0.170635
0.136508

(14)
116.205366
40.145143
0.0924
0.11088

(15)
116.473127
39.998621
0.16405
0.213265

(16)
116.564156
39.919256
0.18376
0.220512

(17)
116.628119
39.82431
0.28185
0.366405

(18)
116.230168
40.071781
0.01988
0.015904

(19)
116.376508
40.031323
0.12084
0.145008

(20)
116.298797
39.95249
0.0098
0.00686

(21)
116.188288
39.785815
0.18068
0.234884

(22)
116.12494
39.789416
0.00144
0.001008

(23)
116.257688
40.099361
0.00138
0.000966

(24)
116.55617
39.936915
0.09835
0.07868

(25)
116.635922
39.832181
0.278405
0.222724

(26)
116.105126
39.74371
0.0029
0.00232

(27)
116.466411
39.988646
0.00249
0.001992

(28)
116.598303
40.050973
0.15183
0.182196

(29)
116.424648
40.127548
0.19423
0.233076

(30)
116.376915
39.971917
0.00147
0.001029

(31)
116.707641
39.91889
0.0336
0.02352

(32)
116.198841
39.708363
0.05849
0.070188

(33)
116.495639
39.911295
0.20116
0.160928

(34)
116.22464
39.782793
0.149455
0.1046185

(35)
116.56836
39.838198
0.028
0.0336

(36)
116.353407
39.932741
0.1076
0.08608

(37)
116.276302
39.780674
0.17524
0.122668

(38)
116.646199
39.879243
0.0344
0.02408

(39)
116.326454
39.884205
0.234
0.2808

(40)
116.470445
39.78931
0.2146
0.25752

(41)
116.454592
39.919255
0.0063
0.00756

(42)
116.189176
39.912679
0.60931
0.792103

(43)
116.16643
39.991848
0.361275
0.28902

(44)
116.312054
39.828131
0.189518
0.1516144

(45)
116.56121
39.97287
0.05432
0.070616

(46)
116.39213
39.747772
0.0146
0.01022

(47)
116.196599
39.715292
0.17024
0.221312

(48)
116.452815
40.040456
0.249676
0.3245788

(49)
116.303852
40.054474
0.15506
0.201578

(50)
116.413854
39.84225
0.43526
0.565838

Table 4: The detailed information of battery swap stations.

BSS
Longitude
Latitude
Level

(51)
116.150573
39.70208
1

(52)
116.526359
39.768352
1

(53)
116.150969
39.937666
2

(54)
116.47133
40.084254
2

(55)
116.394889
39.89192
2

Table 5: The detailed information of dynamic customers.

No.
Longitude
Latitude
Demand

(56)
116.11091
39.965585
0.7481

(57)
116.27524
40.131681
0.4505

(58)
116.53663
40.001629
0.0838

(59)
116.38261
39.850152
0.2289

(60)
116.28443
40.051582
0.6133

(61)
116.35207
40.01801
0.1523

(62)
116.58908
40.047472
0.5258

(63)
116.58674
39.946024
0.5383

(64)
116.38943
39.97676
0.4961

(65)
116.22062
39.774572
0.0781

Figure 6: The geographic location of the nodes.

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

L1
L2
L3
K
C

(1)
389.44
444.76
141.28
6
6460.43

(2)
362.67
473.73
152.37
6
6861.04

(3)
384.80
418.75
142.25
6
6210.01

(4)
384.76
432.69
140.69
6
6333.79

(5)
386.04
450.38
129.27
6
6396.50

(6)
389.41
479.08
129.27
6
6683.48

(7)
388.95
465.88
130.26
6
6561.44

(8)
384.80
418.75
142.25
6
6210.01

(9)
386.04
450.38
129.27
6
6396.50

(10)
364.13
446.60
143.22
6
6498.27

The routes of the best solution are shown in Table 7. 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.

Cycle
Route

C1-1
Route 1
0-2-40-4-52-11-46-37-34-21-47-32-26-22-42-20-7-36-0

Route 2
0-27-15-9-28-6-1-10-31-38-25-17-35-16-24-45-13-33-0

Route 3
0-43-3-14-23-18-49-8-29-54-48-5-19-30-12-41-50-0

Route 4
0-44-39-0

C1-2
Route 1
0-2-40-4-52-11-46-37-34-21-47-32-26-22-42-20-65-0

Route 2
0-27-15-9-28-6-1-10-31-38-25-17-35-16-24-45-0

Route 3
0-43-3-14-23-18-49-8-29-54-48-5-19-30-12-41-64-0

Route 4
0-44-39-59-0

C2
Route 1
0-36-7-61-60-57-56-0

Route 2
0-13-58-62-63-33-50-0

Combined with the real road network in Beijing, the actual driving paths generated by the solution are shown in Figure 7.

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

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 between TS and CW-TS.

TS
CW-TS
Gap (%)

Instance
I
W
B
BestC₁
Time (s)
BestC₂
Time (s)

C101_25
25
3
2
3316.21
3.44
3231.93
3.14
5.53

R101_25
25
3
2
5341.55
6.61
5025.43
6.49
5.92

RC101_25
25
3
2
4896.72
4.91
4562.00
5.94
6.84

C101_50
50
5
3
6433.90
13.32
6040.15
10.80
6.12

R101_50
50
5
3
7963.51
18.54
7426.85
16.67
6.74

RC101_50
50
5
3
9344.88
19.10
8855.83
18.20
5.23

C101_100
100
18
13
15344.92
66.14
13343.35
65.68
13.04

R101_100
100
12
6
14491.10
78.32
12650.19
82.35
12.70

RC101_100
100
12
5
17601.50
76.11
15249.08
66.80
13.37

As can be seen from Table 8, 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.
Coordinate
Demand
BSS
Coordinate

(1)
[40, 40]
18
1
[40, 33]

(2)
[75, 65]
16
2
[31, 80]

(3)
[40, 80]
18
3
[18, 58]

(4)
[9, 68]
11
4
[15, 30]

(5)
[22, 40]
12
5
[45, 20]

(6)
[22, 45]
17
6
[60, 45]

(7)
[28, 65]
12
7
[80, 48]

(8)
[29, 68]
10
8
[55, 71]

(9)
[49, 10]
14
9
[70, 70]

(10)
[48, 13]
19
10
[20, 40]

(11)
[70, 33]
12
11
[25, 20]

(12)
[72, 18]
15
12
[40, 60]

(13)
[80, 25]
13
13
[85, 20]

(14)
[80, 29]
10

(15)
[90, 48]
15

(16)
[20, 20]
18

(17)
[88, 60]
12

(18)
[85, 70]
10

Table 10: Information of dynamic customers and battery swap stations in R101_100.

No.
Coordinate
Demand
BSS
Coordinate

(1)
[10, 10]
18
1
[15, 40]

(2)
[15, 31]
2
2
[20, 13]

(3)
[30, 21]
12
3
[50, 25]

(4)
[15, 12]
12
4
[35, 60]

(5)
[38, 41]
14
5
[13, 70]

(6)
[50, 37]
18
6
[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]
10

Table 11: Information of dynamic customers and battery swap stations in RC101_100.

No.
Coordinate
Demand
BSS
Coordinate

(1)
[30, 80]
18
1
[20, 70]

(2)
[31, 78]
12
2
[20, 40]

(3)
[20, 60]
16
3
[25, 15]

(4)
[22, 61]
18
4
[45, 75]

(5)
[22, 58]
14
5
[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]
10

Table 12: C101_100 optimal solution.

Cycle
Route

C1-1
Route 1
0-75-1-2-4-6-9-11-10-8-7-3-5-0

Route 2
0-13-17-18-19-15-16-14-12-0

Route 3
0-99-100-97-93-92-94-95-96-98-0

Route 4
0-90-87-86-83-82-84-85-88-89-91-0

Route 5
0-66-61-72-80-79-77-73-70-71-76-78-81-0

Route 6
0-67-65-63-74-62-64-55-68-69-0

Route 7
0-59-60-58-56-53-54-57-0

Route 8
0-43-42-41-40-44-46-45-48-50-51-31-52-49-47-0

Route 9
0-29-34-36-39-38-37-35-33-32-0

Route 10
0-21-22-23-26-28-30-27-25-24-20- 0

C1-2
Route 1
0-75-1-2-4-6-9-11-10-8-7-3-5-120-0

Route 2
0-13-17-18-19-15-16-14-12-117-0

Route 3
0-99-100-97-93-92-94-95-96-0

Route 4
0-90-87-86-83-82-84-85-88-89-91-115-0

Route 5
0-66-61-72-80-79-77-73-70-71-76-78-81-127-124-0

Route 6
0-67-65-63-74-62-64-55-68-69-114-0

Route 7
0-59-60-58-56-53-54-123-0

Route 8
0-43-42-41-40-44-46-45-48-50-51-31-52-49-47-118-0

Route 9
0-29-4-36-39-38-37-35-33-32-0

Route 10
0-21-22-3-26-28-30-27-25-24-20-119-0

C2
Route 1
0-121-116-108-98-109-131-130-128-126-125-122-57-129-0

Table 13: R101_100 optimal solution.

Cycle
Route

C1-1
Route 1
0-53-58-40-21-73-74-72-4-54-12-26-28-0

Route 2
0-80-68-3-77-76-50-51-30-70-10-62-0

Route 3
0-5-84-17-45-46-36-49-47-48-82-8-83-60-0

Route 4
0-19-11-105-64-63-90-32-20-66-65-71-35-34-78-9-0

Route 5
0-33-81-79-29-24-55-25-39-67-23-56-75-0

Route 6
0-22-41-15-43-14-44-38-86-16-61-91-100-98-0

Route 7
0-99-59-93-85-37-92-97-87-42-57-2-13-0

Route 8
0-94-95-96-6-89-18-52-7-88-31-1-69-27-0

C1-2
Route 1
0-53-58-40-21-73-74-72-4-54-12-26-28-112-0

Route 2
0-80-68-3-77-76-50-51-30-70-10-62-0

Route 3
0-5-84-17-45-46-36-49-47-48-82-117-105-0

Route 4
0-19-11-105-64-63-90-32-20-66-65-71-35-34-78-9-116-0

Route 5
0-33-81-79-29-24-55-25-39-67-23-56-0

Route 6
0-22-41-15-43-14-44-38-86-16-61-91-100-98-102-110-108-0

Route 7
0-99-59-93-85-37-92-97-87-42-57-2-0

Route 8
0-94-95-96-6-89-18-52-7-88-31-1-69-27-0

C2
Route 1
0-111-113-114-103-115-75-118-109-13-60-83-107-8-0

Table 14: RC101_100 optimal solution.

Cycle
Route

C1-1
Route 1
0-68-61-70-100-55-88-53-10-11-9-0

Route 2
0-13-15-16-17-47-14-12-78-73-79-60-0

Route 3
0-1-3-5-45-4-46-8-7-6-2-0

Route 4
0-42-44-43-40-36-35-37-38-39-0

Route 5
0-41-72-71-105-32-33-89-76-63-85-51-84-56-0

Route 6
0-93-67-34-31-29-27-26-28-30-50-62-95-0

Route 7
0-65-66-83-24-64-91-92-94-96-54-81-80-0

Route 8
0-20-49-19-18-48-21-23-25-22-0

Route 9
0-57-103-77-58-75-97-59-87-74-86-52-99-0

Route 10
0-69-98-82-90-0

C1-2
Route 1
0-68-61-70-100-55-88-53-10-11- 9-111-112-0

Route 2
0-13-15-16-17-47-14-12-78-73-79-60-110-0

Route 3
0-1-3-5-45-4-46-8-7-6-0

Route 4
0-42-44-43-40-36-35-37-38-39-117-0

Route 5
0-41-72-71-105-32-33-89-76-63-85-51-84-56-103-114-0

Route 6
0-93-67-34-31-29-27-26-28-30-50-62-95-116-105-0

Route 7
0-65-66-83-24-64-91-92-94-96-54-81-0

Route 8
0-20-49-19-18-48-21-23-25-22-0

Route 9
0-57-103-77-58-75-97-59-87-74-86-52-0

Route 10
0-69-98-82-90-0

C2
Route 1
0-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 the CW-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.

Data Availability

The relevant parameters of the simulation experiment are shown in Table 2. The customer’s coordinate information is shown in Tables 3–5. The optimal solution in Table 6 can be verified by the path order in Table 7.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the National Social Science Foundation of China (19CGL041) and the Cultural and Social Science Research Project of the Chongqing Education Commission (19JD033).

References

W. R. Stewart and B. L. Golden, “Stochastic vehicle routing: a comprehensive approach,” European Journal of Operational Research, vol. 14, no. 4, pp. 371–385, 1983.

H. N. Psaraftis, “A dynamic programming solution to the single vehicle many-to-many immediate request dial-a-ride problem,” Transportation Science, vol. 14, no. 2, pp. 130–154, 1980.

F. Gonçalves, S. R. Cardoso, and S. Relvas, “Optimization of a distribution network using electric vehicles,” in Proceedings of the IO2011-15 Congresso da Associação Portuguesa de Investigação Operacional, pp. 18–20, Coimbra, Portugal, 2011.

M. Dror, G. Laporte, and P. Trudeau, “Vehicle routing with stochastic demands: properties and solution frameworks,” Transportation Science, vol. 23, no. 3, pp. 166–176, 1989.

C. Gauvin, G. Desaulniers, and M. Gendreau, “A branch-cut-and-price algorithm for the vehicle routing problem with stochastic demands,” Computers & Operations Research, vol. 50, pp. 141–153, 2014.

M. Gendreau, G. Laporte, and R. Séguin, “An exact algorithm for the vehicle routing problem with stochastic demands and customers,” Transportation Science, vol. 29, no. 2, pp. 143–155, 1995.

H. Lei, G. Laporte, and B. Guo, “The capacitated Vehicle routing problem with stochastic demands and time windows,” Computers & Operations Research, vol. 38, no. 12, pp. 1775–1783, 2011.

L. Bertazzi and N. Secomandi, “Faster rollout search for the vehicle routing problem with stochastic demands and restocking,” European Journal of Operational Research, vol. 270, no. 2, pp. 487–497, 2018.

A. G. Qureshi, E. Taniguchi, and T. Yamada, “A microsimulation based analysis of exact solution of dynamic vehicle routing with soft time windows,” Procedia-Social and Behavioral Sciences, vol. 39, pp. 205–216, 2012.

A. M. F. M. Abdallah, D. L. Essam, and R. A. Sarker, “On solving periodic re-optimization dynamic vehicle routing problems,” Applied Soft Computing, vol. 55, pp. 1–12, 2017.

R. J. Kuo, B. S. Wibowo, and F. E. Zulvia, “Application of a fuzzy ant colony system to solve the dynamic vehicle routing problem with uncertain service time,” Applied Mathematical Modelling, vol. 40, no. 23-24, pp. 9990–10001, 2016.

X. Ge, X. Wang, and L. Deng, “Research on open and dynamic vehicle routing problems based on joint distribution,” Journal of Industrial Engineering and Engineering Management, vol. 27, no. 3, pp. 60–68, 2013.

M. Okulewicz and J. Mańdziuk, “The impact of particular components of the PSO-based algorithm solving the dynamic vehicle routing problem,” Applied Soft Computing, vol. 58, pp. 586–604, 2017.

S. F. Ghannadpour, S. Noori, R. Tavakkoli-Moghaddam, and K. Ghoseiri, “A multi-objective dynamic vehicle routing problem with fuzzy time windows: model, solution and application,” Applied Soft Computing, vol. 14, pp. 504–527, 2014.

X. Ge, G. Xue, and P. Wen, “Proactive two-level dynamic distribution routing optimization based on historical data,” Mathematical Problems in Engineering, vol. 2018, Article ID 5191637, 15 pages, 2018.

R. G. Conrad and M. A. Figliozzi, “The recharging vehicle routing problem,” in Proceedings of the 61st Annual Conference and Expo of the Institute of Industrial Engineers, pp. 21–25, Reno, NV, USA, May 2011.

S. Erdoğan and E. Miller-Hooks, “A green vehicle routing problem,” Transportation Research Part E: Logistics and Transportation Review, vol. 48, no. 1, pp. 100–114, 2012.

M. Schneider, A. Stenger, and D. Goeke, “The electric vehicle-routing problem with time windows and recharging stations,” Transportation Science, vol. 48, no. 4, pp. 500–520, 2014.

J. Lin, W. Zhou, and O. Wolfson, “Electric vehicle routing problem,” Transportation Research Procedia, vol. 12, pp. 508–521, 2016.

M. Alejandro and G. Christelle, “The electric vehicle routing problem with nonlinear charging function,” Transportation Research Part B, vol. 103, pp. 87–110, 2017.

D. Goeke and M. Schneider, “Routing a mixed fleet of electric and conventional vehicles,” European Journal of Operational Research, vol. 245, no. 1, pp. 81–99, 2015.

J. Yang and H. Sun, “Battery swap station location-routing problem with capacitated electric vehicles,” Computers & Operations Research, vol. 55, pp. 217–232, 2015.

J. Hof, M. Schneider, and D. Goeke, “Solving the battery swap station location-routing problem with capacitated electric vehicles using an AVNS algorithm for vehicle-routing problems with intermediate stops,” Transportation Research Part B: Methodological, vol. 97, pp. 102–112, 2017.

S. Shao, W. Guan, and B. Ran, “Electric vehicle routing problem with charging time and variable travel time,” Mathematical Problems in Engineering, vol. 2017, Article ID 5098183, 13 pages, 2017.

J. Su, Y. Yu, and X. Zhang, “Knowledge transfer efficiency measurement with applications for open innovation networks,” International Journal of Technology Management, vol. 81, no. 1, pp. 118–142, 2019.

W. Wang, X. Ge, L. Li, and J. Su, “Proactive and reactive multi-project scheduling in uncertain environment,” IEEE Access, vol. 7, pp. 88986–88997, 2019.

X. Ge, X. B. Ge, and W. Wang, “A path-based selection solution approach for the low carbon vehicle routing problem with a time-window constraint,” Applied Science, vol. 10, no. 4, 2020.