Abstract
With energy and environmental issues becoming increasingly prominent, electric vehicles (EVs) have become the important transportation means in the logistics distribution. In the real-world urban road network, there often exist multiple paths between any two locations (depot, customer, and charging station) since the time-dependent travel times. That is, the travel speed of an EV on each path may be different during different time periods, and thus, this paper explicitly considers path selection between two locations in the time-dependent electric vehicle routing problem with time windows, denoted as path flexibility. Therefore, the integrated decision-making should include not only the routing plan but also the path selection, and the interested problem of this paper is a time-dependent electric vehicle routing problem with time windows and path flexibility (TDEVRP-PF). In order to determine the optimal path between any two locations, an optimization model is established with the goal of minimizing the distance and the battery energy consumption associated with travel speed and cargo load. On the basis of the optimal path model, a 0-1 mixed-integer programming model is then formulated to minimize the total travel distance. Hereinafter, an improved version of the variable neighborhood search (VNS) algorithm is utilized to solve the proposed models, in which multithreading technique is adopted to improve the solution efficiency significantly. Ultimately, several numerical experiments are carried out to test the performance of VNS with a view to the conclusion that the improved VNS is effective in finding high-quality distribution schemes consisted of the distribution routes, traveling paths, and charging plans, which are of practical significance to select and arrange EVs for logistics enterprises.
1. Introduction
1.1. Background
In the past two decades, with the rapid development of economy and e-commerce in China, the annual parcel deliveries have been growing continuously, bringing opportunities and challenges for the development of the logistics industry. For example, the annual business volume of express delivery service enterprises reached 50.71 billion pieces in 2018 with a year-on-year growth of 26.6% in China [1].
Logistics vehicles are characterized by high energy consumption, emissions, and pollution. The Organization for Economic Cooperation and Development points out that the freight traffic in major cities of developed countries accounts for 10%–15% of the total urban traffic volume, while the environmental pollution caused by freight vehicles accounts for 40%–60% of the total urban traffic volume [2]. Therefore, as important participants in emission activities, logistic enterprises should fulfill their social responsibilities.
With increasingly prominent energy and environmental problems, electric vehicles (EVs) stand out by virtue of their characteristics of environmental protection, energy-saving, and low cost, and they gradually become an important transportation means in the logistics industry. Hence, electric vehicle routing problem (EVRP) has become a popular research field, which attracts the attention of many scholars. Nevertheless, EVs always have some limitations such as limited battery capacity and long recharging times. With the rapid development of advanced technology, the shorter charging time and the improvement of battery endurance bring opportunities for the wide application of electric vehicles. Furthermore, this paper is in line with the goal of the Ministry of Transport of China to achieve the successful completion of the task of tackling the key problems of transportation pollution prevention and control in 2020, and it also has responded to the call of promoting scientific and technological innovation in transportation and the development of intelligent transportation, as well as intensive logistics.
1.2. Literature Review
To the best of our knowledge, Dantzig and Ramser [3] first proposed the vehicle routing problem (VRP). Clarke and Wright [4] then improved the method proposed by Dantzig and Ramser and presented the C-W saving algorithm. With continuous development in recent years, researchers have considered different influencing constraints to form many branches.
Vehicle routing problem with time window (VRPTW) is a classical VRP problem with the time window constraint. Solomon [5] first considered the time window constraint to vehicle routing and proposed an insertion heuristic algorithm with satisfactory results. After that, a great deal of meta-heuristic algorithms is proposed to solve the VRP. Garcia et al. [6] used the tabu search (TS) algorithm to solve VRPTW and constructed the initial solution with the insertion heuristic algorithm. Gambardella et al. [7] designed the ant colony algorithm to solve VRPTW. Ombuki et al. [8] proposed a hybrid heuristic algorithm based on the genetic algorithm and TS, which is used to optimize the number of vehicles and travel distance. Bräysy [9] proposed a four-stage heuristic algorithm to solve VRPTW based on variable neighborhood search (VNS). Pellegrini et al. [10] studied VRP constrained by different vehicles and multiple time windows and used the ant colony algorithm to solve the proposed model. Recently, Alzaqebah et al. [11] presented artificial bee colony solutions, which are used by the scout bees to memorize the abandoned solutions. Binart et al. [12] addressed the routing problem with stochastic travel and service times by giving a two-stage solution approach. Hiermann et al. [13] considered time windows at customer points, which is a common and important constraint in real-world routing planning problems and solved this problem by means of branch-and-price algorithm. Huang et al. [14] employed route-path approximation and took the Beijing road network as an example to calculate the average speed of distinctive districts during different periods. He and Li [15] investigated the routing optimization problem with the goal of minimizing fuel consumption under the time window and heterogeneous fleet, which is solved by using the tabu search algorithm.
Since the 21st century, the problem of environmental pollution and energy shortage has attracted people’s attention. When people began to look for more energy-saving and clean vehicles, green vehicles emerged at the historic moment. For the past few years, with people’s further attention, the continuous innovation of electric vehicle technology has made it become popular in recent years and become a mainstream of green vehicles. Consequently, a majority of large logistics enterprises have also begun to develop electric vehicle distribution. Therefore, EVRP, which is the combination of EVs and VRP, become a new research direction.
Lin et al. [16] considered that EVRP has comparable traveling time and distance when it is compared with the diesel truck VRP. The authors emphasized the influence of vehicle load on energy consumption and proposed a method of minimum cost including travel cost and energy cost. Goeke and Schneider [17] studied a rich fleet mixing problem, which is traditional fuel oil vehicle and electric vehicle hybrid fleet, in which traveling speed, slope, and load are considered to establish a more practical energy consumption model. They developed an adaptive large neighborhood search (ALNS) approach with an embedded local search procedure which is using a surrogate function to evaluate changes efficiently. Bruglieri et al. [18] proposed a mixed-integer linear programming formulation to assume that the battery recharging level reached at each station is a decision variable in order to guarantee more flexible routes. The model aims to minimize the total travel, waiting and recharging time, and the number of the employed EVs. Chen and Murata [19] presented a novelty model for solving an electric vehicle routing problem with the objective of minimizing the total electrical consumption. As for the energy consumption of electric vehicles, Zheng [20] proposed an energy consumption model related to the rolling resistance coefficient of tires, the air resistance coefficient and the windward area, and a constant of the automobile.
In the study of the EVRP, due to the great influence of congestion and service demand problems, there are a stream of studies on the electric vehicle problem with time widow, which we refer to as E-VRPTW. Penna et al. [21] proposed an iterative local search algorithm coupled with a set partitioning model to solve the electric fleet size and mix vehicle routing problem with time windows and recharging stations. Several types of electrical vehicles with varying driving range, capacity, and fixed cost should serve a set of customers within their time limits, and during their tours, each vehicle can be recharged in stations. Liu et al. [22] introduced a green VRPTW-PF model, which takes the fuel consumption of acceleration and waiting on traffic lights into consideration. Keskin and Çatay [23] formulated a mixed-integer linear programming model by allowing partial recharges with three recharging configurations which can be referred to as normal, fast, and superfast recharges. To solve the larger problems, this paper developed a matheuristic approach which couples the adaptive large neighborhood search approach with an exact method.
There is also a great deal of research studies on the charging of EVs. Schneider et al. [24] considered the limited vehicles and battery capacity; the vehicles consume battery energy according to a constant consumption factor when traveling along the arc. In addition, vehicles are possible to visit charging stations along the route, and the recharging time depends on the current battery charge on arrival at the station. The authors proposed a hybrid algorithm combined with VNS and TS and tested its performance on the instances for the multidepot VRP with interdepot routes. Conrad and Figliozzi [25] proposed the recharge vehicle routing problem, which consists of routing a fleet of EVs, limited by vehicle capacity and time window and limited driving range. Felipe et al. [26] studied an extension of the green VRP, which considered different types of charging stations with different costs and charging speeds. In this paper, a heuristic is constructed, and heuristic algorithm and simulated annealing algorithm are proposed for several local searches. Keskin and Çatay [27] allowed all EVs to recharge partially (not fully) at the charging stations in transit. The paper presented a mathematical model and designed the ALNS algorithm to solve the proposed model. Shao et al. [28] presented the hybrid genetic algorithm to get the routes and charging plan. Koç et al. [29] introduced the electric vehicle routing problem with shared charging stations to minimize the sum of the fixed opening cost of charging stations and the drivers’ cost.
1.3. Contributions of This Paper
In this paper, we are committed to planning reasonable distribution routes to save logistics costs and improve economic efficiency for enterprises. Concretely, we consider the EVRP with time windows and flexible paths. Besides, vehicles in this research are assumed to be fully charged in charging stations.
The main contributions of the paper can be summarized as follows:(1)This paper first considers the path flexibility between any two locations in the EVRP. In the process of distribution, the path selection mainly depends on the departure time from the locations and the congestion levels in the current road network. In this case, the electric vehicle may travel along different paths since the time-dependent speed in the real-world road condition, which is defined as path flexibility (PF) in this paper. In addition, with the consideration of path flexibility, it is able to save battery energy and thus plays an important role in cost savings.(2)This paper provides two mathematical models to formulate the time-dependent electric vehicle routing problem with time windows and path flexibility (TDEVRPTW-PF). The first model is used to select the optimal path between any two locations including consumers and charging stations to achieve the shortest path and the least energy consumption. The second model gives the general travel route of multiple EVs’ distribution, taking path flexibility and time windows into consideration. The total goal is to minimize the total distance during the travel, to save costs, and to increase profits for enterprises.(3)This paper proposes an effective VNS algorithm to solve TDEVRPTW-PF. By inserting virtual distribution centers into the traveling sequence to represent the routes of multiple vehicles, the results can be seen more intuitively. In addition, considering the time and space complexity of the algorithm, the traveling sequence and charging record are separated, and multithreading is adopted, which improves the efficiency significantly. When testing the charging conditions, two locations in advance are used to determine the charging conditions in order to ensure that the vehicle has enough electricity to continue serving and to provide protection for the vehicle’s endurance.
The remainder of this paper is organized as follows: Section 2 gives a statement on problem description and presents two mathematical models. Section 3 proposes an improved VNS algorithm for the solution. Numerical experiments are given in Section 4. Finally, we summarize our work with concluding remarks in Section 5.
2. Problem Statement
This paper focuses on the routing optimization of EVs with time windows and time-dependent travel times, and it also involves charging and battery energy consumption. Before formulating the multiple EVs’ routing optimization model with time windows and path flexibility, we first make the following assumptions:
Assumption 1. The electric vehicles are assumed to be a homogeneous fleet. All vehicles hence have the same load and volume. Furthermore, all EVs are associated with the same maximum battery capacity and recharging rate.
Assumption 2. The EV is always fully charged, and it can be recharged multiple times. Besides, a charging station can hold many EVs simultaneously, and we do not need to consider the queue time for charging.
Assumption 3. According to the actual distribution situations of logistics enterprises, we divide the working time into several time periods, and therefore, the speed on each path between any two locations is constant in a certain time period, so it is expressed in terms of average speed.
Based on the above assumptions, the problem studied in this paper can be formulated as a fleet of EVs fully charged, departing from a distribution center, can serve a series of customers with different demands. Each vehicle delivers its goods to the customers within predetermined time windows. Moreover, in order to prevent an EV that needs to charge from not reaching the charging stations, the energy consumption analysis should be carried out in advance.
Next, we will give an example of how to analyze the energy consumption so as to ensure the vehicle that needs to charge can reach the charging station. Figure 1 illustrates an example involving three EVs (EV1, EV2, and EV3), eight customers (C1–C8), three charging stations (S1–S3), and a distribution center (D). The arrow lines represent the optimal paths between any two locations, and the percentage values on the arrow lines represent the energy consumption. Besides, the percentage values in the box show the battery state when the vehicle arrives at a customer. From Figure 1, we can see that EV1 first serves customers C1 and C2 and returns to the distribution center without any charging; EV2 needs to charge after serving customers C3 and C4; however, it cannot arrive currently at nearest charging station S2. Therefore, EV2 should be recharged at S1 after serving C3, and then EV2 continues to serve C4 and C5. Accordingly, EV3 needs to charge at S3 after servicing C6 and visits S3 again after servicing C7. Finally, it returns to the distribution center after servicing C8. From Figure 1, we can also observe that an EV can be charged many times (such as EV3), and each charging station is not necessarily visited (such as S2).
For the sake of convenience, relevant notations used in the models are listed in Table 1. Let denote customers and denote charging stations, and we set . For simplicity, notations and both denote the distribution center, and each route starts at and ends at . The set of nodes in the distribution network can be denoted by , and we set , . Hence, the actual road network is described by the graph , where A represents the set of arcs. Each arc can correspond to a set of multiple paths in the graph, denoted as . A path , connects nodes i and j in . Each path is associated with distance and travel time . The EV arrives at customer i, , at which has a demand with weight and volume , service time , and time window . According to the actual conditions of delivery, we divide the working time of a day into R time periods, and each period is represented by notation r, . The average speed and travel time in period r are assumed to be and , respectively, when a vehicle travels on path. Therefore, our goal is to optimize the distribution scheme by minimizing the total traveling distance and energy consumption with the consideration of time windows and path flexibility. In the following, we propose two 0-1 mixed-integer programming models of the problem.
2.1. Analysis for Energy Consumption of Electric Vehicles
Referring to Zheng [20], this paper analyzes the energy consumption of electric vehicles (see Figure 2). The energy consumption in unit distance mainly includes two aspects, i.e., driving energy and vehicle accessory energy. In addition, some electric vehicles can recover part of the electrical energy through kinetic energy when braking. The energy consumption during driving is mainly affected by three elements. (1) The working environment, such as temperature and humidity, will affect the power conversion efficiency of the battery. (2) The driving style will also affect the energy consumption. For example, a driver who likes to brake suddenly tends to consume more energy. (3) The vehicles work by overcoming resistance, including rolling resistance and air resistance. The energy consumption of vehicle accessories is mainly affected by two elements, including the average speed of the vehicle in various working environments and the electrical energy consumption of comfort accessories, safety accessories, and functional accessories.
Based on the above discussion, the energy consumed by the driving style and vehicle accessories is difficult to calculate. Therefore, the energy consumption of working by overcoming resistance is firstly calculated, and then the proportion of this energy consumption to total energy consumption is given based on experience and prediction. Finally, the total energy consumption can be obtained according to this proportion.
The calculated process of energy consumption from node i to j on path p in time period r is given in equation (4).
Hereinafter, the notions in equation (4) are defined as in Table 2.
For the charging process, we assume that the charging time is constant, and the EV is always fully charged. Besides, it can be recharged multiple times. The dynamic performance of electric vehicles is reflected on the influence of dynamic speed and different path choices on the distribution process.
2.2. Optimal Path Selection Model between Any Two Locations
In order to determine the optimal path between any two locations including customers, charging stations, and distribution center, we will establish a 0-1 mixed-integer programming model to minimize the energy consumption and the distance.
Firstly, the decision variables can be defined as
Since a vehicle traveling from a location to another one may cover multiple time periods, its total time is the sum of travel time in these time periods. For example, we assume that the distance between two locations is 10 km, and the duration of the time period is one hour, i.e., 60 min, and thus, the step speed function is shown in Figure 3. When an EV arrives at customer i at 8:50, it cannot reach the next customer j at the current speed 40 km/h in this time period, and it have to change the speed to 60 km/h and keep traveling in the second time period from 9:00 to 10:00. Therefore, the distance of the path between locations i and j can also be represented as
Correspondingly, the total energy consumption of an EV is the sum of the consumption in all of the time periods. Thus, we can get the equation as
To select the optimal path on arc (i, j) with time-dependent travel speeds, we can formulate the model as follows:
The objective function (5) is to minimize the energy consumption and the distance, each of which has weights and (). Constraint (6) guarantees that the vehicle selects exactly one path if it travels on the arc . Constraints (7) and (9) enforce the vehicle may go through multiple time periods. Constraint (8) represents the flow conservation. Finally, constraints (10)–(12) define the binary decision variables.
2.3. Multiple Electric Vehicles’ Routing Programming Model with Time Windows
The speeds on different paths between any two locations are not always the same during different time periods, which leads to different energy consumption. We can get the optimal path with least energy consumption and shortest distance between any two locations at any departure time from the above model. Accordingly, we can use the variables on the optimal path between i and j to define the variables on the arc , such as
On the basis of the optimal path selection model, the time-dependent electric vehicle routing problem with time windows and path flexibility (TDEVRP-PF) can be formulated as
The objective function (14) is to minimize the total traveling distance. Constraint (15) handles the connectivity of customers and ensures that each customer is visited by only one vehicle. Constraint (16) imposes that there are no more than K vehicles. Vehicle flow conservation constraint (17) presents the number of vehicles that arrive at a node must equal the number of vehicles that depart for all nodes including customers and charging stations. Constraints (18) and (19) ensure the time feasibility of arcs leaving the customers and the charging stations, respectively. Constraint (20) enforces the time windows of the customers. Constraints (21) and (23) guarantee that weight and volume of cargoes satisfy all customers’ demand. Constraints (22) and (24) ensure that the weight and volume of each vehicle should not exceed the vehicle’s capacity. Constraints (25) and (26) keep track of the battery state of charge and make sure that it is never negative. Finally, constraint (27) defines the binary decision variables.
3. A Variable Neighborhood Search Algorithm for TDEVRPTW-PF
The proposed model is an NP-hard problem, and it involves many factors, such as multivehicle simultaneous delivery, hard time window, path flexibility, and charging on the route. Through comprehensively considering both advantages and disadvantages of various heuristic algorithms (tabu search algorithm, genetic algorithm, simulated annealing algorithm, variable neighborhood search algorithm, etc.), we finally choose the variable neighborhood search (VNS) algorithm which is more inclusive to solve the proposed model. Based on the idea of the VNS algorithm, we improve and perfect it according to the actual needs, so as to solve the proposed model effectively and efficiently. The objective of the algorithm is to minimize the total traveling distance while satisfying the constraints.
The rest of this section can be divided into four parts. The first part shows a way to solve the PF problem. According to the idea of solving the multiobjective decision model, we collect the data of the dynamic road network and select the optimal path between any two locations to generate the database for TDEVRPTW-PF. In the second part, we present a form of the initial solution for simultaneous distribution of multiple vehicles, which is an important part of the algorithm and a prerequisite for data experiments. The third part is to analyze the characteristics of the target problem (that is, the constraints in the presented model) and to filter the solutions. Finally, the whole VNS framework and the termination condition of the algorithm are introduced in the fourth part.
3.1. A Method to Solve the PF Problem
As shown in Figure 4, between any two nodes i and j, the traveling speed of each segment in different time periods may fluctuate under the influence of traffic conditions, namely, the speed of vehicles traveling in different paths between two nodes at the same time may be different, and the speed of vehicles traveling in the same path between two nodes at different times may also be different. With the change of traveling speed, the optimal path between the two nodes may also change.
In order to select the dynamic path, we construct an optimal path selection model. In the algorithm for solving it, we store the average speed of multiple optional paths between any two locations corresponding to different time periods in the database. In other words, the path flexibility (PF) problem influenced by the time factor has been commendably solved through the optimal path selection model.
3.2. Generation of the Initial Solution
In the conventional VRP, vehicle routes are usually represented in the form of customer point sequences. Nevertheless, the solution form of traditional route planning cannot tackle the multiple electric vehicles’ problem appropriately for the VNS algorithm. In TDEVRPTW-PF, we aim at proposing a route planning scheme for multiple electric vehicles to distribute simultaneously, which can both meet the constraints and give a charging scheme. Therefore, we propose a concept of virtual distribution center and fit the drop-shaped distribution route of multiple vehicles into a circular distribution path, which is explained concretely in Figure 5.
Representation of traditional VRP solutions is given in Figure 6.
Representation of the solution after generating virtual distribution centers is given in Figure 7.
The sites in this study that constitute the distribution route include customer sites and charging station sites. The traditional representation of the solution to solve the VRP (as shown in Figure 6) is given as a sequence with the distribution center as the two ends and the customer sites in the middle, which is difficult to satisfy the traveling sequence of multiple vehicles in this study. On the basis of the solution form of the traditional VRP, we use the sequence of mixed sorting of customer sites and virtual DCs (distribution centers) to represent the vehicle route (see Figure 7). Considering the time and space complexity of the algorithm, the charging station sites are not included in the solution, but the charging record, i.e., a vehicle chooses to go to a charging station for charging at a certain customer site, based on the traveling route and charging record stored separately. Instead of performing neighborhood transformation on the entire sequence including the charging station sites, we can transform the sequence containing only customer points and then generate a charging history based on the customer point sequence. In this way, the transformation operation of the solution is simpler, and it is made more convenient to carry out large-scale data experiments, and meanwhile, the time and space efficiency of the algorithm are improved. Ultimately, the form of the solution can be organized by two parts: traveling sequence and charging record.
When generating the initial traveling sequence, the total number of virtual distribution centers in the sequence is equal to the number of vehicles minus 1 (see Figure 7). Meanwhile, the total number of distribution centers in the sequence is equal to the number of vehicles plus 1. First, both ends of the sequence must be 0 and represent the actual distribution center. Then, the virtual distribution centers are inserted into the randomly scrambled customer site sequence while ensuring that they are not adjacent.
In order to obtain an initial feasible solution, a feasibility test is carried out on the generated initial driving route (see Section 3.3), and the charging record in the initial feasible solution is obtained. Hereinafter, a specific example will be taken to illustrate the generation of the initial solution. Suppose that there are 10 customers, numbered from 1 to 10, who are served by a fleet of 3 homogeneous vehicles which are marked by A, B, and C. The vehicles depart from 2 virtual distribution centers marked by 0 and charge in 5 charging stations which are marked from 11 to 15.
The process of generating the initial solution can be introduced as follows. Firstly, the customers are randomly arranged to 5-2-4-6-9-7-8-10-3-1. Additionally, at both ends, two 0’s are inserted to represent the actual distribution center. Then, according to the principle of nonadjacency, two 0’s which represent virtual distribution centers are inserted into the sequence separately. Finally, the initial solution includes charging record which is obtained through feasibility examination.
The initial solution consists of two parts: one is the traveling route and the other is charging record. Now that 0 represents the virtual distribution center, the sequence between each two 0’s denotes a traveling route of a single vehicle. In this case, the sequence (0-5-2-4-0-6-9-7-0-8-10-3-1-0) means the first vehicle A serves customers 5, 2, and 4 in sequence. Likewise, vehicle B serves customers 6, 9, and 7, and vehicle C serves customers 8, 10, 3, and 1 in order. Furthermore, according to the charging history of vehicles, after serving customer 9, vehicle B has to go to charging station 12 to recharge, and vehicle C ought to get energy in charging station 15 after serving customer 3.
3.3. Discriminant of Feasible Solutions
In order to examine the feasibility of the generated solution, a checking procedure is applied to distinguish whether the initial sequence is feasible or not, which can be divided into the following four parts (see the left part of Figure 8).
3.3.1. Decision Number of Vehicles
For an initial sequence, if two 0’s are adjacent to each other in a traveling route sequence, it indicates that the actual number of vehicles is reduced, and the number of vehicles in such an initial sequence does not satisfy the constraints, namely, it is an unfeasible sequence.
3.3.2. Decision Volume and Load of Vehicles
The total weight of freight to be delivered at all customer sites between the two adjacent 0’s of a traveling route sequences is compared with the maximum load of the vehicle. If the load is exceeding the maximum load, the sequence is not feasible. The measurement of volume is similar to that of load.
3.3.3. Generation and Decision of Charging Records
We refer to Zheng [20] to calculate the energy consumption during traveling between two locations, which have been introduced in Section 2.1. Considering the factors such as the driving speed changing with the time period and path, the actual load, rolling resistance, and energy conversion efficiency, the calculated energy consumption is closer to the reality.
As shown in Figure 9, examination of electricity is divided into the following two aspects:(1)According to the energy consumption during traveling to determine whether the vehicle needs to charge. As is shown in Figure 9, suppose the current node of the vehicle is i, and its current electric quantity is SOC. Let the energy consumption from the previous node i − 1 to the current node i be E0. The energy consumption from node i to node i + 1 is E2. Customers i + 1 to i + 2 consume electricity E4. The energy consumption from node i to the optimal charging station (the sum of the distance between the charging station and nodes i and the charge station and i + 1 is the minimum, which is to make the minimum distance to deviate the vehicle from the original route and return to the original route) is E0. If SOC < (E2 + E4) and SOC < (E2 + E3), go to charge station 1 and record the number of charge stations and the customer site number.(2)Determine whether the route sequence is feasible according to the current electric quantity during traveling. When the residual battery capacity of the current vehicle is unable to reach the charge station nearest to the current point, that is, when SOC < E1, it is judged as an unfeasible solution.
3.3.4. Decision of Time Windows
We update the current time every time when vehicles reach a customer site through the time advance mechanism. Comparing the current time with the time window of the customer, if it exceeds the time window range, it is determined as an unfeasible sequence.
3.4. An Improved VNS Component
Mladenovic and Hansen [30] first proposed the variable neighborhood search (VNS) algorithm, which subsequently became a research hotspot in the field of operations research and was applied to the NP-hard problem with good performance. The basic idea of VNS is to obtain the approximate optimal solution by a series of neighborhood operations on the feasible solution. After meeting the termination condition, it is considered that the current approximate optimal solution is the globally approximate optimal solution.
The traditional VNS algorithm has a considerable adaptability to the TSP and other single-vehicle problems, but it is unable to solve multivehicle problems. Therefore, by generating virtual distribution centers, we propose a solution suitable for multivehicle distribution, which has been described in detail in Section 3.2. Next, we will introduce the VNS algorithm structure in detail on the basis of this solution form.
3.4.1. Neighborhood Structure
The neighborhood structure is crucial in the VNS algorithm. The purpose of the neighborhood structure is to constantly transform the current feasible solution so that a large number of other solutions can be derived from the feasible solution in the neighborhood. In the neighborhood structure, the core of the algorithm to transform the feasible solution is the neighborhood operations. Thus, in our algorithm structure, two complementary and effective neighborhood operations are selected after screening so that the operations can be fully traversed in the neighborhood.
(1) Neighborhood Operation 1. Above all, two positions must be determined in the feasible solution sequence S (except the two ends), and then the sequence is reversed in the middle of the two ends to obtain a new traveling sequence , as shown in Figure 10.
(2) Neighborhood Operation 2. First, two positions in the feasible solution sequence ought to be determined (except the two ends), and then the sequence integrally in the middle of the two positions is moved to the beginning part of the overall sequence. It should be noted that the initial point represents the distribution center, so the position of the initial point remains unchanged. The moving part of the sequence is actually inserted after the first point, as shown in Figure 11.
The following describes the procedure of the neighborhood structure algorithm by taking neighborhood structure 1 as an example: Step 1: enter a feasible solution S0, and set the neighborhood approximate optimal solution as S = S0. Step 2: the position in the sequence in which the initial operation is performed is changed to i = 2, j = i + 1. Step 3: set the maximum number M of neighborhood operations. The initial record is m = 0. Step 4: M++; if m ≥ M, stop. Step 5: (1) if the position of i is in the penultimate of the sequence, i = 2, j = i + 1, return to Step 4. (2) If the position of j is at the end of the sequence, i++, j = i + 1. (3) Perform neighborhood operation 1 on the two positions i and j in the feasible solution sequence to generate S′. Step 6: determine whether S′ is feasible. If not, j++, return to Step 5. Step 7: determine whether S′ is better; if better, S = S′. Step 8: j++; return to Step 5.
3.4.2. VND Component
VND is a major component of the VNS algorithm. Its main function is to perform neighborhood operation on the feasible solutions of the input, generate the approximate optimal solutions within the range of multiple neighborhoods of the current solution by turn, and select the new approximate optimal solution through comparison. Two types of neighborhood structures and the corresponding operations have been introduced in Section 3.4.1, which will not be repeated here.
The specific operation steps of VND are as follows: Step 1: enter a feasible solution and assign it to the current solution. Step 2: generate neighborhood 1 by performing neighborhood structure 1 on the current solution, and screen out the local approximate optimal solution in neighborhood 1. Step 3: compare the local approximate optimal solution in neighborhood 1 with the objective function value of the current solution. If the objective function of the local approximate optimal solution is better, replace the current solution with the local approximate optimal solution of neighborhood 1 and return to Step 2; otherwise, proceed to Step 4. Step 4: generate neighborhood 2 by performing neighborhood structure 2 on the current solution, and screen out the local approximate optimal solution in neighborhood 2. Step 5: compare the local approximate optimal solution in neighborhood 2 with the objective function value of the current solution. If the objective function of the local approximate optimal solution value is better, replace the current solution with the local approximate optimal solution of neighborhood 2, and return to Step 2; otherwise, perform Step 6. Step 6: current approximate optimal solution = current solution. Output the current solution.
3.4.3. Shaking Procedure
After VND operation is performed on the current approximate optimal solution, the VND output needs to be compared with the global approximate optimal solution. If the value of the objective function of the current approximate optimal solution in many iterations is greater than that of the global approximate optimal solution, we may get in the situation of the local approximate optimal solution. To skip from local optimization, large transformations are required for the current approximate optimal solution, which is the main significance of the shaking procedure.
For a simple example, as shown in Figure 12, remove “0” of two ends from the sequence 0-5-2-4-0-6-7-9-0-8-10-3-1-0. Thus, it turns to 5-2-4-0-6-7-9-0-8-10-3-1. Subsequently, we divide it into four segments: the first segment is 5-2-4, the second is 0-6-7, the third is 9-0-8, and the fourth is 10-3-1. The original third, fourth, second, and first segments are arranged to form a new current solution 0-9-0-8-10-3-1-0-6-7-5-2-4-0.
3.4.4. An Improved VNS Algorithm for TDEVRPTW-PF
The components of the VNS algorithm have been briefly introduced above. Figure 8 describes the whole procedure of VNS, and the left part of it has been introduced in Section 3.3, while the right part of it illustrates the rest part of the VNS algorithm. The steps of the VNS algorithm are introduced in the following: Step 1: generate the initial route sequence according to the principle. Step 2: check the sum of vehicles; if the conditions are mismatched, return to Step 1. Step 3: check whether the maximum load and volume of the vehicle meet the requirements. If not, return to Step 1. Step 4: generate charging history and check the electric quantity. If the condition is not met, return to Step 1. Step 5: check the time window. If the condition is not met, return to Step 1; if met, take it as a feasible solution. Step 6: set an iteration counter m = 0; set the maximum number of iterations of VND M. Step 7: M++; perform VND operations on the feasible solution to find approximate optimal solutions in multiple neighborhoods. Step 8: determine whether the new feasible solution is better than the current feasible solution. If better, current feasible solution = new feasible solution. Step 9: perform shaking to avoid getting in the local approximate optimal solution. Step 10: if m > M, stop; if not, return to Step 7.
3.5. Double-Optimization Method
In order to improve the optimization degree and stability of the global approximate optimal solution obtained by the VNS algorithm, we generate multiple initial sequences according to the above principles and operate VNS at the same time. By comparing the global approximate optimal solution generated by each group, the best one is taken as the final optimal solution. For the purpose of improving the time efficiency of the whole process, we use multithreading as the main framework of double-optimization when programming and achieved a good result in subsequent experiments and data analysis.
4. Numerical Experiments
In this section, a large number of numerical experiments are conducted to evaluate the performance of our modified VNS algorithm. In the first part, we present the development tools, the experimental environment, and the parameter setting of the numerical experiments. In the second part, a set of solutions are generated under the initial conditions. Then, we evaluate the performance of the solutions and explain the traveling routes in detail. In addition, we also analyze the efficiency of the improved algorithm and the degree of solution optimization. Finally, the sensitivity analyses of battery capacity, load, volume, and the number of EVs are demonstrated in the third part.
4.1. Experimental Environment and Parameter Setting
Our modified VNS is implemented as a code in C++ in Microsoft Visual Studio 2017 software on the Windows 10.0 platform of a personal computer with Intel(R) Core (TM)i7-6500U CPU @ 2.50 GHz and 4 gigabyte memory.
Our data are randomly generated to simulate real distribution scenarios. We randomly select 30 customers, 12 charging stations, and 5 EVs as our experimental instances. Table 3 shows the basic parameters on TDEVRPTW-PF instances. According to the actual conditions, the parameters about energy consumption of EVs are set in Table 4. Parameters about EVs themselves and customers are shown, respectively, in Tables 5 and 6.
The geographical distribution of customers, charge stations, and the distribution center is shown in Figure 13, and Table 6 gives the parameters of customers.
4.2. Experiments on TDEVRPTW-PF Instances
4.2.1. Generation of TDEVRPTW-PF Instances
The basic conditions for conducting numerical experiments are given in the previous part. In this part, the speed is made to fluctuate with time and path, which is limited to the range from 50 km/h to 60 km/h. After that, the VNS algorithm program is run with the maximum number of iterations set, 300 times, and the approximate optimal solution is obtained in Table 7. In the approximate optimal solution, the total traveling distance of 30 customers served by 5 EVs is 1512.4 km. The specific traveling routes and charging records are as follows.
4.2.2. Performance of VNS on TDEVRPTW-PF
We use the generated TDEVRPTW-PF test instances to analyze the performance of our VNS algorithm. Figure 14 shows the specific traveling routes, where it can be seen that the traveling routes of EVs are closed-loop routes with the shape of water drops. It shows obvious division of the regional distribution situation. This is consistent with the distribution arrangement of logistics enterprises in reality, which proves the correctness of our solution.
However, we also find some seemingly unreasonable routes at the same time, such as EV 5 chooses to deliver from customer 13 to customer 9 rather than to customer 14. Instead, EV 2 is needed to return from customer 2 and then head to customer 14. By exploring the reasons, we arrange customer 14 to be behind customer 13. Then, after checking the new sequences, it is found that EV 2 reaches customer 16 at 11:54 which is not in its service time window (12:00–22:00). Similarly, we try to arrange customer 16 to be behind customer 18 and customer 26 to be behind customer 29, however, all of which do not accord with its service time window when arriving at customer 16, and the underlying reason is that all vehicles are set to start to deliver at the same time as stipulated in our preconditions.
On the basis of a series of experiments, we draw a conclusion that the factors affecting the performance of the solutions are mainly divided into two categories. One is the constraint conditions, which include the hard time windows of the customers, the vehicle load constraint, and the vehicle volume constraint. The other is the VNS algorithm itself, which is based on the maximum number of iterations of the algorithm and the maximum number of iterations in the neighborhood. The partial irrationality of the solution affected by the factor of the first type is immutable, but that influenced by the second type can be changed, which will be explained in detail in the next part.
4.2.3. Analysis on the Efficiency of VNS
It can be seen from experiments on TDEVRPTW-PF instances that our improved VNS algorithm has been able to generate the solution we want to solve TDEVRPTW-PF. However, the problem we need to solve is to make the process more efficient and to optimize the distribution scheme. By adjusting the important parameters of the improved algorithm, the further experiments are shown in Table 8.
From the analysis of the experimental results in Table 8, it can be found that the setting of multithreading can significantly improve the optimization degree and robustness of the approximate optimal solution at the expense of certain time efficiency under other constant conditions. Similarly, the maximum number of iterations in the neighborhood has no significant influence on the approximate optimal solution. Finally, when the maximum numbers of iterations in the thread and neighborhood are constant, the maximum number of iterations of the algorithm can improve the optimization of the solution to a certain extent. Under the given experimental environment and acceptable conditions, we finally adopt a variable neighborhood search algorithm with the maximum number of iterations in the neighborhood equal to 30, the maximum number of iterations in the algorithm equal to 300, and the number of threads equal to 10 to improve time efficiency and optimization degree of solving TDEVRPTW-PF.
4.3. Sensitivity Analysis
4.3.1. Sensitivity Analysis of Electricity
In the case of experimental environment and parameter setting of the above examples, the sensitivity analysis is carried out by changing the battery capacity (electricity when an EV is fully charged) of EVs. Table 9 shows the corresponding experimental results of ten groups of different battery capacities.
Experimental results drawn in a two-dimensional coordinate system indicate that the objective function value of the solution decreases as the battery capacity of EVs increases in the range from 120 kwh to 260 kwh. There is a strong correlation between the objective function value of the solution and the battery capacity of EVs, and the correlation coefficient is 0.9683. The reason is that when the number of vehicles is less, the average number of customers that each vehicle needs to serve will increase, which will increase the charging times and the total traveling distance of the scheme. When the battery capacity of an EV is greater than or equal to 260 kwh, the objective function value of the solution fluctuates steadily within a certain range, and an EV can even complete its distribution task without charging. The decreasing trend of the objective function value of the approximate optimal solution also disappears. In terms of the robustness of the solution, the volatility of each approximate optimal solution is about 5%, which is within the acceptable range when the battery capacity is constant.
Through the sensitivity analysis of the battery capacity of electric vehicles, we get a suggestion that according to the actual distribution requirements, logistics enterprises should select EVs with the battery capacity at the inflection point of the C-E curve in Figure 15, which can effectively reduce costs and improve economic efficiency.
4.3.2. Sensitivity Analysis of the Volume
In addition to electricity, the maximum load and maximum volume of EVs also affect the performance of the solution. Because of the similar properties, we select the maximum volume to perform the sensitivity analysis. The experimental results in Table 10 show that when the maximum volume is greater than or equal to 20 , with the increase of the maximum volume, the objective function value of the solution presents a downward trend with volatility. When the maximum volume reaches 29 , the objective function value of the solution no longer decreases. The reason is that the number of customers an EV can serve is directly proportional to the maximum volume of the vehicle, and when the number of vehicles is fixed, the larger the maximum volume of the vehicle is, the more feasible driving sequences are available.
This is of practical significance in vehicle procurement, and we can get the best vehicle volume level by analyzing the inflection point of the C-V curve in Figure 16. Therefore, the type and model of the EVs can be selected on the basis of that, so as to further improve the economic efficiency of the distribution.
4.3.3. Sensitivity Analysis of the Number of EVs
The instances given in Section 4.2.1 show how to deliver and charge in the process of distribution under a given number of EVs. However, in fact, too many or too few EVs serving for certain customers will affect the costs or the quality of service. Therefore, we carry out sensitivity analysis on the number of EVs. Table 11 shows our experimental results. In Figure 17, the trend line C-N is fitted according to the results. In the case of 30 customers, in the process of increasing the number of vehicles from 3, the objective function value of the solution decreases first, then increases, and reaches the minimum value when the number of vehicles is equal to 5. It shows that, under the current distribution condition, the costs of 5 vehicles can be reduced to the lowest level. In addition, the robustness of the objective function value of the solution is the best at this time, and the average absolute deviation of the objective function value of all approximate optimal solutions is about 1.27%. Hence, we should choose the best number of EVs and make a tradeoff according to the robustness required when arranging the distribution plan in practice.
5. Conclusions
In this paper, we present TDEVRPTW-PF with the consideration of time windows, charging, energy consumption, and path flexibility. In order to determine the optimal path between any two locations, a 0-1 mixed-integer programming model is established, minimizing the energy consumption and the distance. Based on the optimal path selection model, the other model is established to solve multiple electric vehicles’ routing problems with time windows.
Then, we choose the VNS algorithm with good effect on the NP-hard problem in the heuristic algorithm as our algorithm framework. We modify and optimize it on the basis of TDEVRPTW-PF. After that, programming and debugging of the program are completed based on the improved VNS.
In the end, a large number of numerical experiments are carried out by running the program. We explain the experimental instances and analyze the efficiency of the algorithm. In addition, the sensitivity analyses are also carried out. By adjusting some parameters such as battery capacity, maximum load, maximum volume, and the number of planned EVs, we draw some conclusions that have practical significance on purchasing operation and vehicle arrangement for logistics enterprises.
However, as a vehicle routing problem, the optimization of TDEVRPTW-PF is never-ending. In the future, further studies may concentrate on energy consumption analysis and the remaining electricity estimation of EVs. Furthermore, the diverse types of the vehicles need to be considered. Further research on this topic may also address that detailed work needs to be executed such as the different departure times of EVs and the insertion of the random distribution task. All of the studies are aimed at making the model established closer to the reality, and the results of scientific research are more capable of serving the society.
Data Availability
The data used to support the findings of this study are included within the paper. Further details may be available to readers from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research work was supported by the National Natural Science Foundation of China (no. 71801018) and the Fundamental Research Funds for the Central Universities (no. 2018RC39).