The electric vehicle routing problem with partial recharge and vehicle recycling

In this paper, a new variant of the electric vehicle (EV) routing problem, which considers heterogeneous EVs, partial recharge, and vehicle recycling, is investigated based on logistic companies' practical operation. A mixed integer linear programming (MILP) model is proposed to formulate the problem. For small-scale scenarios, commercial solver, e.g., CPLEX, is leveraged. For large-scale instances faced by practical applications, a hybrid metaheuristic is designed through integrating a modified Greedy Algorithm with the Variable Neighborhood Search (VNS). The proposed algorithm was tested by real-world instances from JD, an e-commerce enterprise in China. Computational results indicate that partial recharge and vehicle recycling can save costs effectively. It also shows that the number of charging stations is an important factor for the application of EVs.


Introduction
Greenhouse gas (mainly CO 2 ) emissions have spurred a major concern as the past 9 years witnessed the intense growth of daily CO 2 of roughly 5.15% from 388.72 parts per million in 2010 to 408.58 parts per million in 2019 [1. Data from the European Commission (2016) [1] indicate that the transportation system is one of the major drivers for greenhouse gas emissions, accounting for 25% of the CO 2 emissions, and the ratio is expected to be double by 2050. Not surprisingly, regulations promulgated by the European Commission are beginning to promote alternative fuel vehicles (AFV) that adopt greener fuel sources (e.g., bio-diesel, liquid/compressed natural gas, and electricity) to replace traditional internal combustion commercial vehicles (ICCVs) that utilize fossil fuels. Among AFVs, the choice of electric vehicles remains the most attractive.
Studies of NASA [2] show that the use of EVs entails the following merits: (1) EVs exhaust fewer carbon emissions, which remarkably improve the air quality; (2) EVs require lower operational and maintenance costs. (3) EVs can offer quicker, quieter, and smoother acceleration, and therefore can be heavily used in densely populated areas with noise restrictions. These advantages of EVs and government policy support have led to a dramatic growth of the EV industry. EVs are beginning to be widely used in logistics companies' delivery fleets [3] such as TNT Express, JD, UPS, etc. However, a major challenge faced by EVs is their limited driving range, which requires additional charging facilities and longer charging time. As a ramification of these restrictions, a new variant of the vehicle routing problem (VRP)-electric vehicle routing problem (E-VRP) arises.
Nowadays, along with the development of e-commerce, logistics has been recognized as one of the key factors in our daily lives. The primary concern of companies is how to save the cost effectively under many practical constraints. Generally, EV's fixed acquisition cost is very high, and the recycling use of EVs can effectively save costs by increasing vehicle utilization rate and reducing necessary labor costs. Most of the previous works do not consider the vehicle recycling and ignore the volume limit, a practical constraint tantamount to the weight constraint. However, in the practical operation of logistic companies, e.g., JD logistics, both the weight and volume of the goods have to be considered, and the EVs are allowed to return to the warehouse and restart again after a short period of time.
Focused on vehicle recycling and volume constraint, a new variant of E-VRP (H-EVRP-PR&VR) is investigated, 1 3 which comprehensively considers mixed fleet, time windows, and vehicle capacities (weight and volume), partial recharge (PR), and vehicle recycling (VR). These constraints are faced by many logistic companies employing EV, complicating the problem devastatingly. For large logistic companies like JD, the EV routes are needed to be optimized for daily package distribution in hundreds of cities, which is an impetus to developing efficient approaches to such kind of problems.
The main contributions of this study are summarized as follows. • We introduce the H-EVRP-PR&VR, which incorporates a fleet of heterogeneous EVs, partial recharge, and vehicle recycling schemes. The objective is to minimize the total cost. • The MILP model for the problem is developed and solved by CPLEX. • We propose a hybrid metaheuristic combining a construction heuristic with a VNS algorithm. New perturbation operators and neighborhood structures are designed and incorporated in the VNS metaheuristic.
The remainder of the paper is organized as follows. The section "Literature review" presents a review of related works. In the section "Model development", an MILP formulation of H-EVRP-PR&VR is presented. We describe the VNS algorithm in the section "The hybrid metaheuristic algorithm" and reports the experimental results in the section "Numerical experiments". The section "Conclusions and future works" concludes the paper and presents a discussion on future research directions.

Literature review
Recently, there is a considerable amount of literature on E-VRP. Erdelić et al. [4] presented a comprehensive review on existing studies for E-VRP. The literatures related with the investigated H-EVRP-PR&VR are reviewed in this section for the convenience of readers.
The Recharging VRP (RVRP), where vehicles with limited range are allowed to recharge at the customer, was firstly introduced by Conrad and Figliozzi [5]. This is a problem that not only the traditional routing decisions are optimized with respect to the load capacity of the vehicle, but also the selection of charging stations should be considered with respect to the limited traveling range of EV. The number of vehicles and the total cost are minimized. Erdogan and Miller-Hooks [6] proposed GVRP by considering charging stations and optimized the problem by minimizing the total distance. Two heuristics named modified Clarke and Wright Savings (MCWS) and Density-Based Clustering Algorithm (DBCA) were presented in the paper. Koc et al. [7] developed a branch-and-cut algorithm which combines a simulated annealing heuristic algorithm and several valid inequalities. Still, considering the possibility of two consecutive charging stations, Bahrami et al. [8] introduced two MILP formulations and consolidated them by some dominant criteria and valid inequalities. Montoya et al. [9] proposed a multi-sampling heuristic, the method includes two phases: sampling and assembly. In the sampling phase, a pool of route is built via a "route first and cluster second" heuristic algorithm, while in the assembly phase, the solution is obtained by solving a set partitioning formulation over the pool. Andelmin et al. [10] solved the problem by a multi-start local search heuristic.
Furthermore, Schneider et al. [11] extended the model to E-VRPTW by adding time windows on the customer and assumed that each vehicle can visits at most one charging station per route and becomes fully recharged with a nonfixed recharging time, and solved the problem by a hybrid heuristic algorithm combining the VNS with tabu search (TS). Keskin et al. [12] developed a two-phase metaheuristic combining the Adaptive Large Neighborhood Search (ALNS) with an exact approach. Felipe et al. [13] studied the GVRP with partial recharge (GVRP-PR) by considering partial recharge, in which three charging technologies with varied charging speed were also considered. Local Search (LS) combined with non-deterministic Simulated Annealing (SA) was utilized to solve it. Based on their work, Ding et al. [14] considered the capacity limitation of the recharging station by restricting the number of EVs to be recharged at a recharging station, and a hybrid heuristic method combining VNS and Tabu Search (TS) was developed. Keskin et al. [15] investigated the E-VRPTW with partial recharge (E-VRPTWPR) and solved it by ALNS. Results show that partial recharge strategies can greatly improve routing decisions. Desaulniers et al. [16] presented branch price and cut algorithms for four variants of E-VRPTW which were different in recharge numbers per route or whether allowed partial recharge. Schiffer et al. [17] extended the problem by considering the location of charging stations. Mao et al. [18] studied E-VRPTW with multiple recharging options and solved it by an improved ant colony optimization algorithm.
Goeke et al. [19] extended the work to E-VRPTWMF by considering the mixed fleet of battery electric vehicles (BEVs) and ICCVs, in which mixed fleet combining a single type of BEVs and ICCVs is considered, and the loaddependent energy consumption is minimized. An ALNS method that allows infeasible solutions during the search is introduced. Hiermann et al. [20] extended the problem by focusing on heterogeneous BEVs with the difference in transport capacity, battery size, and acquisition cost. Besides, in [21], they considered E-VRP which combined ICCVs, BEVs, and plug-in hybrid vehicles, and solved the problem with a metaheuristic approach combining a genetic algorithm (GA) and LNS. Lebeau et al. [22] studied the combinations of different vehicle class and vehicle technology. However, the recharging was only allowed at the depot. Kopfer et al. [23] proposed the mixed use of ICCVs and BEVs, and illustrated the advantaged of mixed fleet by comparing range, payload, and efficiency. Sassi et al. [24] extended Hiermann's work by allowing multiple charging and partial recharge technologies and considering timedependent charging costs. Macrina et al. [25] investigated the mixed use of EVs and ICCVs with partial recharge and figured out the problem with an LNS algorithm.
Additionally, most of the literatures consider using linear functions to approximate charging process. Montoya et al. [26] established a mixed integer programming named E-VRP with nonlinear functions (E-VRPNL) and solved it by a hybrid heuristic method. They used a piecewise linear function to capture the nonlinear behavior of the charging process. Koç et al. [27] extended the work by considering a shared charging station which aims to minimize the aggregate driver cost and construction cost of charging stations. In [28], Froger et al. considered a piecewise linear approximation by tracking the time and the state of charge by an arc instead of a node. An overview of different variants of E-VRP is summarized in Table 1.
To sum up, most of the prior works do not take the volume constraint and vehicle recycling into account, which would be useful for the practical application in logistic companies. Therefore, a vehicle routing problem with E-VRPTW constraints, volume constraints, heterogeneous fleet, partial recharge, as well as vehicle recycling is presented in this paper.

Model development
In this section, the problem is defined and an MILP model that considers vehicle recycling and volume constraints along with partial recharge is consequently developed.

Problem analysis
H-EVRP-PR&VR is an extension of the E-VRP by considering a combination of mixed fleet, time windows, vehicle GVRP Traveling distance Koc et al. [7] GVRP Traveling distance Bahrami et al. [8] GVRP √ Traveling distance Montoya et al. [9] GVRP Traveling distance Andelmin et al. [10] GVRP √ Traveling distance Schneider et al. [11] E-VRPTW √ √ Traveling distance Keskin et al. [12] E-VRPTWPR √ √ √ Traveling distance Felipe et al. [13] GVRP √ √ Total recharging costs Ding et al. [14] GVRP √ √ √ Traveling distance Keskin et al. [15] E-VRPTW-PR √ √ √ Total energy cost Desaulniers et al. [16] E-VRPTW √ √ √ √ Total routing costs Schiffer et al. [17] E-VRPTW-PR √ √ √ Total costs Mao et al. [18] E-VRPTW √ √ √ Total costs Goeke et al. [19] E-VRPTW √ √ √ √ Distance and costs Hiermann et al. [20] E-VRPTW √ √ √ Total costs Hiermann et al. [21] E-VRPTW √ √ √ Total costs Lebeau et al. [22] E-VRPTW √ √ √ Total costs Kopfer et al. [23] E-VRPTW √ √ √ Total energy Sassi et al. [24] HEVRP-TDMF √ √ √ Total costs Macrina et al. [25] GMFVRPPRTW √ √ √ √ Total costs Montoya et al. [26] E-VRP-NL √ Total time Koç et al. [27] E-VRP-NL √ Total costs Froger et al. [28] E-VRP-NL-C √ Total time capacities (weight and volume), partial charging, and vehicle recycling. The problem is defined by a directed graph G = (N, A) with a set of nodes N = C ∪ F ∪ O and a set of arcs A = {(i, j)|i, j ∈ N, i ≠ j . L e t C = {1, 2, … , n} a n d F = {1, 2, … , m} be customers, and the set of recharging stations and their copies, respectively. Vertex O denotes the depot. Each customer i ∈ C entails a capacity p i , a volume q i , a service time s i , and a hard time window [e i , l i ] . Each arc (i, j) ∈ A has a relevant distance d ij and time t ij . There are K types of vehicles available for transport. Each vehicle type k ∈ K encompasses a maximum load capacity Q k , a maximum load volume V k , a maximum energy capacity Y k , a travel cost per kilometer r k , and a fixed cost f k . The maximum number of each type of vehicle is V.
The objective is to find a route plan departing from the depot satisfying customers' demands and arriving at the depot. The vehicles can be recharged at charging stations or depot with non-fixed charging time at the charging station and fixed charging time at the depot. The aim is to minimize the total cost comprised of fixed costs, charging costs, travel costs, and waiting costs. We need to decide vehicle type, visiting sequence, departure time, recharging quantity, and return time for each vehicle. Figure 1 illustrates an example which includes 9 customers (C1-C9), 4 stations (S1-S4), and the depot (D). This case contains two types of vehicle routes marked in different colors, and the percentage on each path represents battery state. Vehicle 1 visits S1 after servicing C1, recharges its battery to 90% before visiting C2 and C3, and consequently goes back to the depot with the battery fully consumed after servicing C4. On the other hand, Vehicle 1 is recharged twice at S1, indicating that a charging station can be visited more than once. Vehicle 2 returns to the depot after serving C5 and C6 and sets off again, via C7-C9. Obviously, the recycling of the vehicle saves a vehicle.
As the proposition proposed by Keskin et al. in [12], we assume that the vehicle departs from the depot with full energy and arrives at the depot with its energy fully consumed if it has been recharged along the route. For convenience, we note the proposition as follows.
Proposition 1 If an optimal solution exists, such that an EV leaves the depot with its battery partially charged, then the same fully charged EV departing from the depot is also optimal.

Model formulation
In this section, a mixed linear model of the problem is formulated with the notations of the model presented at the outset. Then, a mixed linear programming (MLP) is formulated and a proposition, which proves that the model can be transformed into an MILP model with the path consistent with the non-integer programming model, is given.

Notations
The parameters and variables involved in the article are shown in Table 2

Model formulation
In this model, we use a four-dimension binary variable x kl ij to define the type of used vehicle and the number. We also define binary variable x kl iOj to denote whether the vehicle is recycled. kl i represents the cost of time recharged at i by lth vehicle of type k. In this problem, the objective to minimize is given as: The objective function consists of four parts: (1) the total cost of the vehicle, (2) the total travel cost, (3) the total  Model constraints are listed as follows: (1)

3
Equation (1) ensures that each customer should be served by one vehicle and exactly once, while Eq. (2) guarantees flow conservation, i.e., all the inflows to a node must be equal to the outflows from that node, this constraint also ensures that all routes start and end at the depot when j = O . Equation (3) describes that for each type k , the number of vehicles departing again from the depot is less than the number of vehicles departing for the first time. The relationship between x kl iOj and x kl ij is covered by Eqs. (4) and (5). The timing constraints for the vehicle are covered by Eqs. (6)- (13). Equation (6) describes that the starting time θ i at node i is later than arrival time i . Equation (7) ensures that the starting time i is inside the time window [e i , l i ] . Equation (8) ensures the arrival time j is equal to the starting time j when the vehicle departs for the first time, i.e., the service can begin as soon as the vehicle arrives, which differs from Eqs. (9) and (10), in which we consider the vehicle recycling. Equations (9) and (10) consider that i is a customer and charging station, respectively. Equation (9) guarantees that if the vehicle restarts from the depot, its arrival time at j is equal to starting time i plus service time s i , travel time t iO and t Oj , and adjustment time at the depot s O . Equation (10) considers the same case for charging station by replace service time with recharge time kl j . Equation (11) ensures that the arrival time of next node j depends on the starting time at the previous node j plus service time s i and travel time t ij . Equation (12) considers the case when the previous node i is a charging station. Equation (13) restricts the starting time i to be equal to arrival time i when i is either charging station or depot.
The constraints of weight and volume are covered by Eqs. (14)- (19). Equation (14) guarantees that the load volume of the vehicle at first node q kl j is equal to the volume of the node q j . Equation (15) indicates that the load volume of next node q kl j is equal to the load volume of previous node q kl i plus demand q j . Equation (16) guarantees that the load volume q kl i is a positive number and never exceeds the maximum volume of vehicle type k . Equations (17)- (19) are corresponding constraints for load.
Equations (20)- (25) show the restrictions on current available energy y kl ij . Equation (20) guarantees the vehicle can reach next node j from current node i . Equation (21) ensures that if j is a customer or depot node, current available energy y kl ij at j is equal to the previous available energy (32) x kl iOi = 0, ∀k ∈ K, ∀l ≤ V, ∀i ∈ N, y kl si at i minus the distance d ij . Equation (22) considers when j is a charging station, current available energy y kl ij depends on the previous available energy y kl si at i plus the amount of energy the vehicle recharged at station i and minus the distance d ij . Equation (23) guarantees that each time the vehicle leaves the depot with full energy. Equation (24) enforces the current available energy y kl j never exceeds the maximum capacity Y k for each vehicle k . Equation (25) guarantees that the recharge time does not exceed a large enough number.
Equations (26) and (27) place a once-only recycle restriction on each vehicle. Equation (28) makes sure that the vehicle is used in order. Equations (29)-(33) define the domains of variables x kl ij and x kl iOj . By requiring the charging time kl i to be integer, the above model can be easily converted to an integer programming model.

Proposition 2 If is an optimal solution to the MLP, then ⌈ ⌉ is the optimal solution to the corresponding MILP.
Proof If is optimal, then vehicles shall arrive at the depot with its battery fully consumed (if it is charged among the route). ⌈ ⌉ is the smallest integer that not only meets the charging requirements but also minimizes the remaining battery when the vehicle returns to the depot. As the two sides of the time window [e j , l j ] are integers, assume the vehicle visits node j after accessing charging station i , and then, the arrival time of node j satisfies ⌈ j ⌉ − j = ⌈ i ⌉ − i < 1 . If j < j , ⌈ j ⌉ = j . If j = j , then ⌈ j ⌉ − j < 1 . Since j ≤ l j , therefore ⌈ j ⌉ ≤ l j , i.e., the new solution also satisfies the time window, ⌈ ⌉ must be optimal for the corresponding integer linear programming.

The hybrid metaheuristic algorithm
To solve the H-EVRP-PR&VR problem, a hybrid metaheuristic combining a construction heuristic with the VNS algorithm is developed. The former incorporates the idea of the greedy method to construct the initial solution, while the latter improves the solution in each iteration, VNS randomly perturbs the current solution according to the defined neighborhood structure. Then, the departure time is adjusted by a departure time operator which aims to save the waiting cost. The pseudo-code of the hybrid metaheuristic algorithm is shown in Algorithm 1.

Variable neighborhood search
Our VNS structure is similar to the classic VNS that consists of a variable neighborhood descent and a shaking procedure; however, we define new neighborhood structures and a perturbation operator to expand the search space. New labeling algorithm and recycling operator are also developed in the heuristic.

Labeling algorithm
A labeling algorithm is proposed to find out the optimal charging station. For each route, we determine whether it needs to be plugged into a charging station, record the inaccessible location, and find the closest non-customer point. The interval where the charging station can be inserted is found and the location with the lowest insertion cost is chosen. The recharging time is determined by dividing the distance from the current charging station to the next charging station by the power consumption. The algorithm is as follows:

Generation of initial solution
The initial solution for the H-EVRP-PR&VR problem is constructed by a modified Greedy Algorithm (MGA). The first point on each route is selected among the remaining customers according to the urgency of the customer, which is defined as follows: This formula shows that the urgency of the customer U i is determined by the difference between the latest service time l i and the time required to reach the customer from the depot t Oi .
The vehicle type with the lowest cost is selected. After the first point is confirmed, the point with the least travel costs and waiting costs is chosen, until any constraints are violated. While the vehicle cannot reach the next client, the nearest charging station is selected to recharge the vehicle.
exchanged randomly with the cost saving recorded and ultimately, the assignment algorithm is used to find the optimal exchange. The algorithm is as follows:

Perturbation algorithm
A new perturbation operator is proposed with each route divided into three segments according to a randomly given time period. The middle parts of the two routes are 1 3 The saving value defaults to a large number M , and thus, the diagonal entries of the saving matrix are all M. The perturbation algorithm changes one value in the matrix at a time.

Neighborhood structures
As we know, VNS is an improved local search algorithm that makes use of the neighborhood structure formed by different actions to perform alternate searches and diversify the search space. Our neighborhood structures are defined on the basis of the exchange operator. We proposed k-neighborhood structures by the number of exchange nodes with each k-neighborhood structures representing the exchange of k * 100 customers. The exchange cost is calculated for each customer, and the exchange with the greatest saving is selected. Neighborhoods are defined in terms of latitude and longitude as well as the time window coincidence.
Classical operators' Inter-route 2-opt and Relocate are also used in the local search to improve the solution. Only feasible solutions are accepted throughout the local search, but during the insertion process of the Inter-route 2-opt and Relocate, the mileage constraint is allowed to be violated, i.e., after the insert, the remaining mileage of the vehicle cannot reach a certain customer point. Then, the infeasible path is repaired by the Labeling algorithm; if it cannot be repaired, the insert is abandoned.

Recycling operator
A vehicle recycling operator is proposed to save fixed costs by reducing the number of vehicles and improve vehicle utilization. For two routes, if the arrival time of one route is earlier than the departure time of the other, reconstruct the route based on the customers contains in the two routes. The above procedure is carried out until no further improvements can be done.

Departure time adjustment heuristic
As the departure time affects waiting time, an algorithm is designed to adjust the departure time. The algorithm uses the difference between the arrival time and time window of each point to calculate the adjustment time.
For route r , we calculate the waiting time and the delay time of each customer n i , i.e., w i = max{0, e i − i } and sl i = max{0, l i − i } , and figure out the maximum waiting time and minimum delay time for all customers. Afterward, we compute the maximum waiting time and minimum delay time for all customers, and the smallest value among them is chosen as the adjustment time.

Numerical experiments
In this section, we first present an illustration for the data set provided by JD Logistics and solve six small instances by the commercial software, CPLEX. The performance of the  Fig. 2 The best-known solution for the instance with eight customers

3
improved VNS algorithm is tested. The impact of vehicle recycling, partial recharge, and the number of charging stations is analyzed. The VNS algorithm is coded in Python.
All the experiments are implemented on a laptop with Intel(R) Core(TM) i7 processor (2.2 GHz) and 8 GB RAM.

Experiment description
The algorithm is tested on the real-world data of Beijing provided by JD. The logistics company provides urban distribution services to 1000 customer nodes distributed in the whole city. Each customer has a time window. Two types of EVs are employed, and each EV has fixed capacities on both weight and volume. There are 100 charging stations distributed in the whole region, which can provide charging services for all EVs. The EV departs from the logistics center after 8 o'clock and returns before 24 o'clock. When an EV returns earlier, it can depart again after a period of 1 h ( Table 3).
The table below lists out the main information about the two types of EVs.
Some notes on other parameters are as follows: 1. Vehicles departing from the distribution center with battery fully recharged. The EVs are allowed to return to the depot and restart again after 1 h and calculating waiting costs. 2. The unloading time is 0.5 h, while the loading time is not taken into account. 3. The charging cost is 100 yuan/h. 4. There are no restrictions on the charging station. The charging rate g is 4000 km per minute. 5. The waiting cost w is 24 yuan/h.

Computational results of small instances
The correctness of the model is demonstrated by combining the MATLAB toolbox YALMIP with the commercial solver CPLEX for solve small-scale cases. As this problem is raised for the first time, existing benchmark results are not available. A small dataset which includes 11 points, 8 customers, 2 charging stations, and one depot is constructed based on  the scale of the actual data. There are two types of vehicles available, with different vehicle loads and maximum mileage. The distribution of customer points and the best-known solution are shown in Fig. 2. The plane coordinates are converted from latitude and longitude using the depot as the origin. Lines of different colors represent different vehicle types and green dots represent charging stations. In Table 4, the results of commercial solver and VNS are compared for different sizes of customers. The second and third columns record the solution and computation time of CPLEX and the last two columns are the counterparts of the VNS. The results show that for customer points 3-7, VNS returns the same solution as CPLEX, but with a significantly shorter computation time; for 8-customer scenario, VNS can get a better solution than CPLEX in less than 1 s. Figure 3 shows the change in computational time as the customer increases. From the above results, we can see that the calculation time increases exponentially as the number of customers increases. As the number of customers increases to 8, it is almost infeasible to find the optimal solution within an allotted time. Therefore, CPLEX is impossible to find the optimal solution for the large-scale problem, and the corresponding heuristic algorithm needs to be designed.

Computational results of the real-world cases
In this section, the performance of the algorithm on different scales is presented, and the effects of partial recharge, vehicle cycling, and changes in the number of charging stations are analyzed.

Performance of the algorithm
To verify the performance of the algorithm, we modify the MCW algorithm proposed by Schneider et al. [11] as a comparison. The gap between each solution, considering the definition proposed by Hiermann et al. in [21], is calculated by gap = solution−VNS solution VNS solution .
In Table 5, we present the results of MCW, MGA, and VNS, respectively. The second column shows the number of routes with each vehicle type shown in brackets. The total costs are displayed in the third column, while the gaps are shown in the fourth column. The results suggested that the solution obtained by VNS is optimal, and the gap between VNS and MCW is 9.6%, while the gap between MGA and the final solution is 18.9%.
We construct new instances by randomly sampling 500-950 customers from the example. It is more critical for companies to obtain a feasible solution quickly than the optimal solution in practical applications. Therefore, we fix the number of iterations K max as 5 and limit the maximum running time to 10 min. The results are shown in Table 6.
To avoid accidental errors, all experiments are performed ten times with the best of them taken as a result. The distribution of the solutions is presented in Fig. 4. It can be seen from the graph that the box widths are narrow and has fewer abnormalities, reflecting that the solutions are stable and  Fig. 5 The number of utilized vehicles with/without EV recycling Fig. 6 Comparison of the charging cost for partial and full recharge strategies evenly distributed, while the cost reduction trend is steady at around 6% as the number of customers decreases. In Table 6, it can be observed that VNS improves the solution by roughly 9.91%. For each 50 point difference, the cost will drop by about 11,900. The widest gap is more than 15,900, while the smallest gap is 7278.

Effect of partial recharge and vehicle recycling
In Table 7, we show the effects of vehicle recycling and partial recharge on the large-scale case. The second column is the result without considering partial recharge and vehicle recycling, while the third column is the solution considering partial recharge and vehicle recycling. These results indicate that partial recharge and recycling could significantly save the total cost, reducing the initial and final solution by 5.6 and 4%, respectively.
In Fig. 5, the impact of the vehicle recycling on instances with different number of customers is shown. The number of vehicles between recycling and non-recycling is compared. The results demonstrate that vehicle recycling can save roughly 4.0% of the fixed cost, and in fact, the recycling of vehicles also reduces the vehicle purchase requirements and driver costs.
In Fig. 6, we display the comparison between partial recharge and fully recharge. For each instance, partial recharge can save about 1700 costs. The results demonstrate that partial recharge can save about 57.9% of the charging fee.

Effect of charging stations
In E-VRP, another key factor that affects the costs is the number of charging stations. We study the impact of the charging station by varying its number from 100 to 50.
In Table 8, CS stands for the number of charging stations, while Cus represents customers' number. It can be drawn from Table 8 that as the number of charging stations decreases, the optimal solution cannot be obtained for largescale points, and the overall solution projects an increasing trend. For scenarios with less than 700 customers, the optimal solutions are obtained given 70 or 60 charging stations. This means that for smaller scale, 60 charging stations are already sufficient, and due to the reduction in the number of customers, it is easier to get the optimal solution. However, when the number of charging stations is further reduced, the costs increase.

Conclusions and future works
This study introduces the heterogeneous electric vehicle routing problem with time windows, partial recharge, and vehicle recycling. The objective is to minimize the total cost comprised of fixed costs, driving costs, waiting costs, and charging costs. A mathematical model is established. Experiments indicate that commercial solvers, e.g., CPLEX, can only solve small instances of the problem. For largescale instances, a two-stage algorithm is designed. An initial feasible solution is constructed by the MGA heuristic and then is improved by VNS. A new shaking operator and a departure time adjustment heuristic are also included in the VNS to improve the search efficiency. Computational results show that the VNS algorithm can greatly improve the initial solution. Besides, the effectiveness of partial recharge and vehicle recycling is demonstrated, while the impact of the number of charging stations is analyzed.
There are several directions for future research. As the proposed problem is a new variant of E-VRP and integrates many practical constraints faced by logistic companies in the daily operation, the problem turns out to be complicated and computationally expensive. Thus, it remains an essential impetus to study more efficient algorithms for the problem. The problem presented in this paper can also be extended by considering the multi-depot situation and different charging technologies.