SVND Enhanced Metaheuristic for Plug-In Hybrid Electric Vehicle Routing Problem

: Plug-in Hybrid Electric Vehicles (PHEVs), as a new type of environmental-friendly low cost transportation, have attracted growing interests for logistics. The path-planning optimization for PHEV has become a major challenge. In fact, PHEV-based routing optimization is a type of hybrid vehicle routing problem (HVRP). Compared with the traditional Traveling Salesman Problem (TSP) and Vehicle Routing Problem (VRP), the PHEV routing problem should consider more constraints, such as time limits, capacity constraints (including fuel tank capacity and battery capacity), electric stations, fuel stations and so forth. In this paper, a Mixed Integer Linear Programming formulation is presented and a novel hybrid metaheuristic approach (HMA_SVND) is proposed. Our method is a combination of memetic algorithm (MA), sequential variable neighborhood descent (SVND) and a revised 2_opt method. Comparative studies show that our proposed method outperformed previous works.


Introduction
Transportation plays an important role in the logistics system. Actually, transportation for logistic mainly contained two phases: sending raw materials from suppliers to manufacturers and the delivery of finished products from manufacturers to customers [1]. However, current transportation mostly employs fuel powered vehicles, which emit a lot of greenhouse gases (GHG) and harmful substances, such as carbon dioxide and sulfide, causing various environmental problems. Therefore, energy saving and emission reduction are more than necessary for logistic transportation. Recently, the conception of green logistics has been proposed. Green logistics is more environment-friendly. It has quickly become a hot research subject. By optimizing the usage of logistic resources, it focuses on restraining the environmental harm while reducing the whole logistic cost.
In order to build green logistics, the first objective is to make a reasonable layout and path planning. This is usually regarded as a Vehicle Routing Problem (VRP). The VRP is a well-known NP-hard combinatorial optimization problem. It was first introduced in 1959 [2]. The VRP problem is different from traditional routing problems such as Traveling Salesman Problem (TSP), it can be more complicated by taking consideration of multiple vehicles and various constraints. The classical VRP problem can be summarized as follows [3]-there is a certain number of customers with different cargo needs. Distribution centers provide cargo to customers. A fleet is responsible for distributing cargo and organizing appropriate driving routes. The aim is to satisfy the customer's needs and to achieve goals like minimizing the travel distance, reducing costs or optimizing travel time under certain constraints. Lin et al. [4] provided a review study on VRP and classified them into two categories-traditional VRP and Green Vehicle Routing Problem (GVRP) which is firstly proposed by Erdogan and Miller-Hooks [5].

•
The PHEV routing problem is formulated as a Mixed Integer Linear Programming model, different from previous research, various hybrid vehicles are considered.

•
This study proposes a new hybrid metaheuristic which combines the memetic algorithm with a powerful local search technique SVND; the experimental results show that our algorithm obtains better results than previous ones.
The rest of this article is organized as follows. Section 2 reviews the relevant literature. Section 3 defines the problem. Section 4 describes the proposed memetic algorithm and SVND neighborhood search method. Section 5 reports the experimental results and the conclusion is given.

Literature Review
The Vehicle Routing Problem (VRP) was first introduced in 1959. Since then, many variants of the VRP problem have been studied. Recently, the environmental impacts of logistics transportation have grabbed more attention. Green logistics has therefore become a hot research topic. It was first proposed by Dekker et al. [12] and Sbihi et al. [13], and it takes environmental problems into combinatorial optimization problems. Indeed, green logistics has attracted increasing research attention. For example, the GVRP was proposed by Erdogan and Miller Hooks et al. [5], and focuses on reducing emissions and eliminating environmental impacts. In their work, they considered the limited fuel capacity of vehicles and the possibility of refueling at fuel stations along the way. In recent years, there have been many studies on the variants of GVRP and there are some studies combining carbon emissions or energy consumption with different constraints and optimization objectives. For example, Messaoud et al. [14] proposed the Dynamic Green Vehicle Routing Problem (DGVRP), which combines the Dynamic Vehicle Routing Problem (DVRP) with GVRP. Leu et al. [15] took into account the distribution efficiency and carbon emissions and developed a cargo flow modeling method to analyze logistics planning. Ubeda et al. [16] conducted a study aimed at minimization of both traveling distances and pollutant emissions. For GVRP, there are two kinds of solutions based on exact and heuristic algorithms. Zhou et al. [17] proposed a Lagrangian Relaxation-Based Solution for GVRP. Zhang et al. [18] proposed a new tabu search named the RS-TS algorithm to solve the Vehicle Routing Problem with fuel consumption and carbon emissions. Demir et al. [19] used an adaptive large neighborhood search algorithm (ALNS) to minimize driving time and fuel consumption, respectively. Lately, many works have focused on the path planning problem of new energy vehicles. The path planning problem for electric vehicles is called the electric vehicle routing problem (EVRP). Schiffer et al. [20] proposed a location routing method for charging stations, which considered both the path of electric vehicles and the location of charging stations. On the basis of EVRP, Shao et al. [21] considered the impact of an electric vehicle's driving speed and cargo load on energy consumption. They adopted a hybrid genetic algorithm to solve the problem. For EVRP, latter works also considered time window constraints. Roberti and Wen [22] investigated the Electric Traveling Salesman Problem with Time Windows. Keskin and Çatay [23] discussed the electric vehicle routing problem with time-windows and partial recharging (EVRP-TWPR). In recent works on EVRP, the objectives have been to reduce energy consumption, cost, distance and so on. In terms of constraints, the layout of charging stations, partial charging, time window and variable speed are considered. For example, Mavrovouniotis et al. [24] applied the ant colony algorithm (ACO) for EVRP which estimated whether the electric vehicle has a charging station in its working range. Schneider et al. [25] introduced a hybrid heuristic for EVRP with time windows and recharging stations. Hannan et al. [26] provided a comprehensive survey of HEVs' source combination, models, energy management systems and so forth. Wen et al. [27] discussed the plug-in electric vehicles (PEVs) charging selection problem.
Plug-in hybrid vehicles are regarded as one of the most environment-friendly solutions for logistic transportation. Actually, there are more related works that try to provide a solution for HVRP. The vehicle routing problem for new hybrid energy vehicles is termed HVRP, which was proposed by Mancini [6]. The author defined HVRP as an extension of GVRP, in which vehicles can use both electricity and fuel, the unitary travel cost is much lower for distances covered in electric mode, a mixed integer linear programming formulation and a large neighborhood search based on matheuristic are proposed. Yu et al. [8] proposed a simulated annealing (SA) algorithm based on restart strategy for HVRP. Hiermann et al. [9] presented a complex meta-heuristic method, which combines a genetic algorithm with a local and large-scale neighborhood search. Macrina et al. [28] proposed an iterative local search heuristic algorithm to optimize the routing of hybrid vehicles with electric engines and internal combustion engines. Sun and Zhou [29] proposed a Cost Optimization Algorithm (COA) based on the fact that plug-in hybrid electric vehicles have two power sources, so even along the same route, different power choices will lead to different energy consumption. In Tarek Abdallah's thesis [30], the vehicle routing problem for PHEV with time window constraint is studied. The objective of optimization is to minimize the fuel powered driving time within the time window. A Lagrange relaxation method is used to solve this problem, also a new tabu search algorithm is proposed. To our knowledge, Tarek Abdallah was the first to put forward the vehicle routing problem for PHEV.
Generally, HVRP is a combinatorial optimization problem. The solutions can be divided into three categories-neighborhood-centered, population-based and hybrid methods [31]. The first type of approach is generally proceeded by iteratively exploring neighborhoods of a single incumbent solution. Such as Simulated Annealing (SA) [8], Adaptive Large Neighborhood Search (ALNS) [32], Large Neighborhood Search (LNS) [5]. The second type is usually based on natural laws, such as Genetic Algorithms (GA), Evolutionary Algorithms (EA) and Memetic Algorithm (MA) [33]. The third combines the previous methods, it is a hybridization of the metaheuristic and a powerful local search technique. Vidal (2012) et al. [34] proposed a new approach based on a hybrid genetic algorithm for Multidepot and periodic vehicle routing problem. The algorithm combined population evolutionary search, meta-heuristic neighborhood search and adjusted population diversity management to improve the solution's general performance.
As a matter of fact, there are many studies about the vehicle routing problem of electric vehicles, but few articles about the vehicle routing problem of PHEV. To some extent, this paper presents a Mixed Integer Linear Programming model and provides a new algorithm for the vehicle routing problem of PHEV. In our previous work [35], we employed a memetic algorithm for solving the PHEV based HVRP problem by considering both electric and fuel stations. In order to obtain a better solution for a large-scale HVRP with PHEV, we have made some improvements on the previous algorithm. We proposed a novel hybrid memetic algorithm. It mainly combines memetic algorithm with sequential variable neighborhood descent (HMA_SVND). We tested our new approach on larger test instances and detailed comparative studies on open source benchmarks are performed.

Problem Description and Mathematical Formulation
Inspired by the works of Yu et al. [8], Erdogan [5] and Mancini [6], we describe the problem as follows-the PHEV routing problem can be defined as an undirected complete graph G = (V, E), where vertex set V is a combination of the customer set I = {v 1 , v 2 , . . . , v n }, the depot v 0 , a set of fuel stations (S f ) and a set of electric stations (S e ). The fuel station or electric station can be visited several times by the vehicle. Each node j is characterized by a service time p j . Each edge (v i , v j ) is associated with a non-negative travel time t ij and distance d ij . This problem seeks to find at most m tours, visiting a subset of vertices including S e and S f when needed such that the total transportation cost is minimized. The studied problem considers two kinds of constraints: time constraint and capacity constraint (including fuel tank capacity and electric capacity). The time constraint means that the duration of the route in a feasible solution should not exceed the T max . Admittedly, in certain cases, there is a possibility that all customers are visited by one vehicle. The following assumption is considered in this work.

1.
All customers need to be visited and each customer can only be visited once.

2.
There are at most m vehicles available and each vehicle must start and finish in the depot.

3.
Vehicles can enter S e or S f when needed. The station type Z(S z ) can be visited more than once or not visited at all.

4.
There are no restrictions on the number of vehicles accessing electric or fuel stations.

5.
When the vehicle visits the power station, the corresponding fuel or electric energy will be filled up. 6.
Once the vehicle returns to the depot, both its fuel and battery will be filled to its maximum. 7.
Each route may be covered using only the electric engine, or a part of it can be covered using traditional fuel propulsion, we give energy priority to electricity power. 8.
Travel speeds are assumed to be constant over a link. 9.
The travel time of each vehicle used in this problem should not exceed the predetermined maximum time T max .
The mathematical formulation of PHEV routing problem is as follows: Objective function: Subject to: f ijkz ≤ Q kz ∀i, j ∈ V 0 , ∀k ∈ M, ∀z ∈ Z, and i = j (11) Objective (1) is to minimize the total cost of travel using PHEVs in a given day. Constraint (2) and (3) ensure that each customer is visited exactly once by one vehicle and each station may be visited at most once, respectively. Route continuity is guaranteed by constraint (4) while a maximum number of routes is imposed by constraints (5) and (6).
The arrival time at each vertex by one vehicle is tracked through constraint (7). Constraint (8) ensures that each vehicle returns to the depot before the time limit T max . Constraint (9) implies the lower and upper bound on arrival time at customer and station vertices. Constraint (10) tracks a vehicle's energy level based on vertex sequence and type. If vertex j is visited right after vertex i by vehicle k (X ijk = 1), then the vehicle's remaining energy level should minus the fuel and electric consumed on this route. Constraints (11) and (12) implies the upper bound of energy consumption of type z from vertex i to vertex j by vehicle k. Constraint (13) implies that vehicle k's energy consumption from vertex i to vertex j is based on the distance and vehicle k's energy consumption rate. Constraint (14) resets the power level of type z to Q kz after visiting the depot or the station vertex of type z. Finally, the domain of the variables are specified by the constraints (15) to (18).

Methodology
A memetic algorithm is a hybridization of a genetic algorithm with a local search, so this work proposes a memetic algorithm using an SVND.

Memetic Algorithm
Memetic Algorithm (MA), a population-based metaheuristic, is an improvement of the evolutionary algorithm (EA). It is defined by the collaboration between cultural evolution and natural evolution to explore and exploit the search [36]. The concept of MA was inspired by Dawkin's notion of memes (a unit of cultural evolution), which consists of self-improvement [37]. The classical memetic algorithm is a combination of EA and local search procedure (LSP). The pseudo code of standard MA is given in Algorithm 1.
Firstly, an initial population is generated randomly and a local search procedure is applied to optimize it. Then, the fitness values of all the individuals are calculated by using an evaluation process. Afterwards, during the whole loop, the selection, crossover and mutate operators are employed to generate new offsprings. These offsprings are further optimized through the local search procedure, then these solutions are evaluated through the evaluation process. After that, the SelectNewPop operator tries to choose the best individuals for the population of the next generation. These steps are repeated until the stopping criterion is satisfied.

Representation and Evaluation
The solution consists of n customers and q(q ≥ 0) selected fuel or electric stations, where the power stations come from sets S e and S f . Each customer must be visited and can only visited once, but each electric station and fuel station may be visited more than once or not at all. Therefore, the number of visited power stations in solutions may be different. Index 0 represents the depot, the i th number in 1, 2, ..., n denotes the i th customer to be served. The following number of n + 1, n + 2, ..., n + m then denotes the index of power stations. Figure 1 shows two feasible solutions.
An example in Figure 1 shows the representation for an instance where 15 customers need to be visited and 4 power stations were considered, indexes 1 to 15 represent the customers, 16 and 17 electric stations, 18 and 19 fuel stations. As can be seen from Figure 1a,b, we adopt a chromosome coding method in which −1 is used as the delimiter between routes. In Figure 1a only 15 customers were visited without access to the energy station and 4 vehicles were used to visit all customers. In Figure  1b, all customers were visited by 3 vehicles, and the second vehicle and the third vehicle visit a electric station and a fuel station in their tour, respectively. Comparing two representations, the latter uses fewer vehicles because there are vehicles visiting the power station on the way, their driving distance capacity was extended. Of course, their driving time still should not exceed T max . As we considered more constrains than GVRP [5] and HVRP [6] in this work, there would be more infeasible solutions during the computation of the objective function; however, the elimination of all infeasible solutions would harm the diversity of the population and the solution would highly likely fall into local optimum. Therefore, we introduced a penalty mechanism to maintain the diversity of population. The evaluation of the solution can be described as follows-the delimiter is a sign to distinguish different routes. For each route, we need to determine whether time constraint and vehicle energy capacity constraint are satisfied. We set the time constraint as a soft time constraint, and the solution beyond T max is also acceptable. Only a certain penalty should be added to reduce the probability of the solution selected, the objective function of the solution with penalty P will become worse and be eliminated in the process of iterating to produce the new solution. The reason for introducing penalty P is to ensure the diversity of population. However, when the residual electric and fuel energy of the vehicle are not enough to support its return to the depot or to the nearest energy station, the solution is not feasible. The feasibility checking of the solution refers to the checking of the vehicle energy capacity constraint, which is a constraint condition. If the constraint is violated, the solution is not feasible. We added an infinite number to the objective function to cancel the probability of the solution being selected.the pseudo-code for evaluation procedure is shown in Algorithm 2, where the symbols vn, EPCost i , FPCost i , time_spent i represent the number of the route, the cost of electric, the cost of fuel and the time spent in route i, respectively.

Algorithm 2: Evaluation of the individual
Input: An individual Output: The value of objective function Begin Cost←0

Individual Initialization
Generally, the initial solutions are generated randomly but the convergence of the solution can be accelerated by using a heuristic. In this work, in order to speed up the convergence, we generate an initial population including a heuristic solutions and random solutions to explore the different regions of solution space. A single individual can be generated by a novel approach which is revised from Nearest Neighborhood Heuristic (NNH) [38] and the other individuals are generated randomly. The detailed procedure of this method is shown in Algorithm 3.

Select Individual for the Recombine Operator and The Mutate Operator
Inspired by Singh and Banghel [39], Singh and Gupta [40], to avoid premature convergence of programs, we have used a binary tournament individual selection based on the probability p better . The candidate with better fitness is selected with this probability. It is a pre-set value, and this value is suggested to 0.8 [41]. The binary tournament selection takes two individuals randomly from the population, it means that in each selection process, we first generate a random number. If this random number is smaller than the p better , we choose the better one in the two individuals, otherwise we choose the worst one. This process will be repeated twice to select two different individual as the parents for the recombine operator and the mutate operator. The pseudo code of the binary tournament selection is given in Algorithm 4.

Algorithm 3: Method to generate a signal individual
Input: The data of problem Output: An initial solution R Begin

Recombine Operator
Before the recombine operator, we obtained two parent individuals by using binary tournament selection method. In this paper we attempt to exchange memes between two parents to obtain new individuals. The recombine operator was proposed firstly by Singh and Baghel [39] to solve a multiple traveling salesperson problem (MTSP) where the number of vehicles is fixed. In our studied problem, due to the uncertainty of vehicle number and the choice of the power station, the individual may have different tour numbers and different individual lengths. Thus, a new revised recombine operator is proposed to solve the studied problem, the details of the reorganization operation are as follows: Let x 1 and x 2 be the two selected parents to recombine operator and let m 1 and m 2 be the number of tours present in x 1 and x 2 , and mR is defined as the minimum tour number mR = min(m 1 , m 2 ). The recombine operator can be divided into the following steps. The first step, sorting all tours by their distances in an ascending way within x 1 and x 2 , and then it constructs a offspring individual x in an iterative way. In each iteration, it selects one of the two parent individuals randomly and moves the i th tour to x . This process is repeated in mR times in our case and the parents are selected with equal probability. Clearly, at the end of first stage, some customers still remain unassigned and they will be moved into an unassigned list.
In the second stage, we adopt the method proposed by Singh and Baghel [39] to generate good solutions. During each iteration, the customers in the unassigned list are reassigned one by one to x . Each time it assigns a customer to a particular tour whose distance increases the least by this reassignment. To achieve this, we have to check all possible insertion positions in every tour. Simultaneously, in order to avoid having only the depot in the route of the generated offspring x , a repair operation is performed. In this repair operation, we choose a tour randomly that has more than two customers, and move one customer to the route which has only one depot. The pseudo code of the recombine algorithm is given in Algorithm 5.

Mutate Operator
In the mutate operator, all tours in x are sorted in ascending order. Each customer from a tour i is moved to x with varied probability P copy , otherwise the customers are moved to the unassigned list. This process is repeated until all customers are moved from x. Then, the customers in the unassigned list will be reassigned to x by the same process in the recombine operator.
Venkatesh and Singh [41] explained why they did not use the constant probability, but the probability related to the number of iterations and they defined P copy in detail, it can be determined by the following equation: Where iter is the running time or the number of iterations, iter max is the maximum running time or the maximum number of iterations, P maxcopy and P mincopy represent respectively the maximum value and the minimum value of P copy . The algorithm for mutate is given in Algorithm 6.

Sequential Variable Neighborhood Descent
Variable Neighborhood Search (VNS) is a metaheuristic based on a systematic change of the neighborhood structures within a search introduced by Mladenović and Hansen [42]. It is based on the idea of a systematic change of neighborhood both in a descent phase to find a local optimum and in a perturbation phase to get out of the corresponding valley [43]. It has been proved that it could solve various combinatorial optimization problems. The Sequential Variable Neighborhood Descent (SVND) is a sample variant of the basic VNS where the systematic changes of the neighborhood structures are performed in a deterministic way. The pseudo code of SVND is shown in Algorithm 7. It should be noted that there are two ways to explore local search, one is called first improvement strategy (a move made when an improvement in the neighborhood is once found) and the other is called the best improvement strategy (a move to the best solution in the neighborhood) [44]. In this work, we adopted the former.
The SVND in this paper is applied as a local search algorithm to improve the performance of a global search algorithm, and five neighborhood structures are proposed in this study (l max = 5). Among all the new individuals produced, the SVND is used to improve the best one among them. The SVND method starts from a given solution x, and the lth neighborhood structure is applied to generate a new solution x , then the fitness value f (x ) is compared with the current value f (x). If the new value is better, the solution x moves to x and continue the search with premier neighborhood (l ← 1), otherwise the next neighborhood structure is applied with l ← l + 1. If the search operation runs l max times without finding better values. The process will be terminated. Specifically, during the search process, an extra search technique 2_opt is applied when the new solution x is better than the incumbent, it generates x and update the solution x only when f (x ) < f (x). The 2_opt method and the neighborhood structure in the SVND algorithm are described in the following.

2_Opt Method
The 2_opt method used in this paper was proposed by Lin et al. [45]. It is a simplified form of the Lin-Kernighan algorithm [46], which is known as the k-opt method. In this paper, we can get a better solution by constantly comparing the distances of customers in the cycle. If the following equation is satisfied, then the indices of targets i + 1 and j should be swapped. Such comparisons can be repeated many times over a route, number of repetitions equals to the length of route-2. The same operation finally updates all routes. The pseudo-code for the 2_opt method is given in Algorithm 8.

Neighborhood Search Structure
The performance of the local search procedure depends on the choice of neighborhood structures. The size of the neighborhood is a key factor; the larger the neighborhood, the more likely it could contain better solutions [46]. We utilize five neighborhood structures in this work, for example, Swap, Insertion, Reverse, Insertion stations and Delete stations. The first three neighborhood structures are common changes in other VRP papers, and the latter two are unique to the hybrid vehicle routing problem in this paper. These operators are briefly explained below and illustrated in Figure 2.

1.
Swap: The swap move is carried out by exchanging the positions of randomly selected i and j nodes of solution x. It should be noted that i and j may come from the same route or may come from different routes. That is to say, the swap operation may be an Intra-route movement or an Inter-route movement.

2.
Insertion: The insertion is implemented by randomly selecting a node from index i and inserting it immediately following a randomly selected j node of x. Similar to swap, insertion operation may be an Intra-route movement or an Inter-route movement.

3.
Reverse: A reverse is performed by randomly selecting i and j nodes from the same route of solution x and then reversing the substring in-between them. 4.
Insertion Stations: An insert station is executed by choosing a station each from S e and S f , and then inserting them into the preceding random i index on solution x. Fuel stations, electric stations and insertion locations are randomly selected. 5.
Delete Stations: A delete station means randomly taking out a power station from solution x Before performing this operation, it is necessary to determine whether there is an power stations in the solution, if there are power stations in the solution, select a deletion from the existing power station, and if not, the operation does not make any changes to the solution.

Hybrid Memetic Algorithm Based on SVND for Solving PHEV Routing Problem
In conjunction with the above sections, the overall scheme of the proposed algorithm in this paper is shown in Figure 3. To make our solution more applicable to the PHEV routing problem, we made some customized improvement on the Memetic algorithm. First, for population initialization, we proposed a new algorithm based on the nearest neighborhood heuristic, it takes the factors of time, vehicle capacitance and fuel capacity as well as the access to the energy station into account. Secondly, we introduced penalty mechanism to maintain the diversity of population. Last, we improved the traditional 2_opt method with a cyclic comparison technique. The algorithm can be summarized as the following steps.

1.
Initialization operation is used to generate a set of solutions, an initial solution is generated by rules and others are generated randomly. The total number of initial solutions is defined as P s .

2.
The 2_opt method is used to improve the generated solutions.

3.
The objective function evaluation is carried out on these solutions.

4.
The binary tournament selection operation is applied to choose the parents, then recombine operator or mutate operator is executed to generate new offspring. This step is repeated P s times, in order to generate a new population with P s individuals.

5.
Evaluating the new population by using objective evaluation. 6.
The SVND procedure is used to further improve the best individual at each generation, which searches all its neighborhoods according to a sequence and returns the best individual among all the searched neighbors. 7.
Select the best individuals from the current population and the offsprings, they are considered as the population for the next generation, and the population size are always kept to P s . 8.
If the stopping criterion is satisfied, the algorithm will be terminated and export the result, otherwise go to step 4.
Steps 4 to 7 are executed in the iteration loop, in which the better solutions have a larger probability to be selected for the recombine operator or the mutate operator. It is noted that the recombine and mutation operators are used in a mutually exclusive manner, each offspring is generated by a recombine operator or a mutate operator but not both. The former is employed with probability P c , otherwise the latter operator is used. The number of new individuals generated in step 4 of each iteration is P s . In step 7, we arrange the total solutions in ascending order and select the non-repetitive P s individuals in turn, if the objective values of the two solutions are identical but the compositions of the solutions are different, then both solutions should be retained, provided that they will be arranged in P s . Finally, the termination of the algorithm needs to satisfy the stopping condition.

Computational Experiments
The experimental results and analysis are described in detail in this section. A comparative study is also conducted with the work of Yu et al. [8]. The proposed method was coded in C++ language and compiled on Microsoft Visual Studio 2017. the computing platform configuration are as following: Intel (R) Core (TM) i7 at 4.20 GHz, 16 GB of RAM with 64-bit platform under Windows 10 Operating System.
The vehicle data used in this paper are the same as the work of Yu et al. [8]. The relevant data are shown in Table 1. The test data we used for testing is a classic example of open source vehicle routing problem. It is an extension of the classical benchmark of VRP. The data of these instances are derived from Capacitated Vehicle Routing Problem (CVRP). This instances consists of datasets for CVRP proposed by Augerat et al.(datasets A, B, and P), by Christofides and Eilon (for dataset E), by Fisher(for dataset F). All instances are available in http://neo.lcc.uma.es/vrp/vrp-instances/. In order to better study the impact of electric stations and fuel stations on PHEV routing problem research, we divided the test into three situations, taking different number of electric stations and fuel stations for testing. The first case consists of 0 electric station and 0 fuel station (E0F0). The second scenario considers 2 electric stations and 2 fuel stations (E2F2) and the third case has 4 electric stations and 4 fuel stations (E4F4). In order to facilitate the comparison with existing methods, we use the same station locations as in the work of Yu et al. [8]. Running tests were conducted on 14 kinds of test data under these three conditions, and their costs and distances were recorded. Specifically, for CVRP, it can be understood as a specific case of our E0F0 situation in this work where the fuel and electric stations are not considered. By replacing the time constraint with load constraint, some modifications to the proposed algorithm can be used to solve CVRP. Compared with HVRP defined by Simona Mancini, our optimization goal is to minimize the total cost and our article covers all the constraints of their articles and adds a few additional constraints. It is shown that the proposed algorithm can solve the problems with energy stations and general vehicle routing problems without energy stations.
Parameter sensitivity analysis has been done in the work of Venkatesh and Singh, Singh and Gupta. Therefore, in this work, our parameter settings are consistent with theirs which are P c = 0.5, P better = 0.8, P mincopy = 0.15, P maxcopy = 0.9, iter max = 1000. When E2F2 and E4F4 are used, l max is set to 5 and has five neighborhood structures. When E0F0 is used, l max is set to 3. At this time, insert and delete station are not performed. For the layout of the fuel station and electric station, the methods we used are as follows-E2F2, 1 and 2 for charging stations, 3 and 4 for gas stations, then the customer sequence will start from 5, and −1 for depot. Similarly, we can get the distribution of E4F4.
Yu et al. [8] compared Plug-in Hybrid Electric Vehicle (PHEV) with Mild Hybrid Vehicles and Conventional Vehicles. The results suggested that PHEV has more advantages than the others. Thus, we considered the PHEV in our work. This paper focuses on the PHEV routing problem, which considers a series of constraints in the VRP problem such as time limits, fuel capacity, electric capacity, fuel station and electric station. Table  2 presents the experimental results of the proposed algorithm for this problem with different power station arrangement (e.g., E0F0, E2F2 and E4F4). We have 14 instances tested under these cases, and each configuration is tested 10 times. The size of the population is 1000, the number of iterations is set to 100. With the increasing number of customers in test cases, the average CPU time of these instances are between 3 and 30 s. The best value of cost, the mean value of cost, the related mile value of the best obtained solution and the mean value of mile for each configuration are presented in Table 2. Since the objective of this paper is to minimize the total cost, the related mile in the table refer to the distances corresponding to the solution with best cost value, instead of the shortest distance in the 10 runs. In most cases, the related mile is less than the average mile in Table 2, but there are few exceptions, such as instance B-35-K5 under E4F4 where the related mile is 408 miles and the average mile is 401.2 m. In this paper, electricity or fuel can be selected for each route, and the HVRP is an extension of GVRP, so we give priority to using electric energy. Therefore, as long as the electric stations are visited enough times along a route, there may be a long distance but low costs for the route.
In order to verify the performance of the proposed algorithm, we considered simultaneously the precision and the robustness. Table 3 shows the comparison results of HMA_SVND and SA_RS_CF to solve the 14 instances with different power station arrangements. The SA_RS_CF is the method proposed by Yu et al. [8]; it means that simulated annealing with a restart strategy and using the Cauchy Function to determine the acceptance probability of poor solutions. The best obtained solution with the lowest cost and related mile of the different methods are presented in this table. We proved that HMA_SVND stability analysis presented in Table 4 where the standard deviation and error value have been shown. In Table 3, we have shown the comparison results of HMA_SVND and SA_RS_CF. The first line of Table 3 is explained. This problem is defined with 33 nodes, in the situation of E0F0, there are 1 depot and 32 customers to be visited, and this problem is considered as a classical CVRP problem. In the case of E2F2, there is 1 depot, 2 electric stations, 2 fuel stations and 28 customers, their indexes are arranged in the above order, in the case of E4F4, there is also a similar circumstance with 4 electric stations and 4 fuel stations. The cost value and mile value of HMA_SVND and SA_RS_CF correspond their best found solution. In the first line, the HMA_SVND has obtained better solutions than SA_RS_CF under both 3 situations. For most cases in Table 3, HMA_SVND has found better solution than the other (39 cases were better of the 42 configurations), the comparison results proved the effectiveness of the HMA_SVND algorithm.
In order to present the comparison results more clearly, we draw the bar chart shown in Figure 4, where Figure 4a-c are the results of the comparison under E0F0, E2F2 and E4F4, respectively. According to Figure 4, we can clearly observe the effectiveness of our proposed algorithm HMA_SVND, which outperforms the SA_RS_CF algorithm to a great extent. According to these graphs, it is not difficult to find that in some benchmark problems, our proposed algorithm can greatly improve the solution, some are slightly improved for most cases.
In order to prove the robustness of the algorithm, the standard deviation (SD) and error of HMA_SVND for tested configurations are shown in Table 4, the latter is defined by the following equation where Aver is the mean value of the results for 10 runs, Opt is the best value (related value for mile) for the 10 runs. As shown in Table 4, in the cases of E0F0, the maximum cost Error is 4.2% and the minimum cost Error is 0.8%, the SD values of cost are between 0.35 and 3.46, in the cases of E2F2, the cost SD values and cost error values are between 0.51 and 0.38, 0.8% and 4.8% respectively, in the cases of E4F4, the cost SD values are from 0.2 to 4.39, the cost error values are from 0.3% to 5.4%. In this study, our aim is to minimize the total cost, the value of mile is the total distance of the solution which has the minimum cost value. The total distance is not the considered point of our problem; therefore, the stability of mile values for the tested instances is less than that of the cost value. In a special case of the B-n35-K5 in E4F4, the mile error is −1.6%, this is because the related mile is greater than the average value in this case. The standard deviation and the error could be used to express the degree of data dispersion. The smaller the two values are, the more centralized distribution of experimental results is, on the contrary, the more decentralized the distribution is.  Figure 4. The comparative results of the proposed HMA_SVND method and SA_RS_CF method of Yu et al. [8] In terms of the robustness, Figure 5a,b graphically illustrate the box plot of the standard deviation and error of cost for different power station configurations. It can be observed that, in the case of E0F0 and E2F2, the stability of the algorithm is better when the number of energy stations increases, the robustness of the algorithm will become slightly worse.
In summary, compared to SA_RS_CF algorithm, our proposed HMA_SVND algorithm obtained a better performance in terms of precision and the robustness.

Conclusions
In this paper, we focused on solving the routing problem for PHEV, a new hybrid optimization algorithm which merges the memetic algorithm and the sequential variable neighborhood descent is proposed. The proposed approach exhibits a competitive performance in terms of precision and robustness, achieving superior results to prevalent algorithms. However, our studied problem is still based on some simplified assumptions. For example, the charging time of the charging station was not considered and the electric station is set to be always available. In real life, charging time and mode should be considered, and there may be queues at stations; vehicles may need to wait for charging and refilling in energy stations, so the cost of waiting should also be considered. At the same time, each customer may have a time window for their services. Therefore, in future research, we will concentrate on more comprehensive problems for real-life PHEV routing.
Author Contributions: The research work was conducted in collaboration with all authors. X.L., X.S. and H.L. defined the research theme and designed the proposed method. X.L., X.S. and Y.Z. implemented the prototype and wrote the article. X.S. and Y.D. discussed and evaluated the data test. Both authors have thoroughly contributed to reading and validating the manuscript. All authors have read and agreed to the published version of the manuscript. Time of arrival at vertex j from vertex i Q kz Vehicle k's tank capacity of power type z r kz Vehicle k's consumption rate (gallons per mile) on power type z c z Cost rate for power source type z, which includes c e and c f T max Pre-determined maximum time for travel from depot until return to depot C A large positive constant value Decision variables x ijk Binary variable equal to 1 if vehicle k travels from vertex i to vertex j; otherwise, 0 y jkz The remaining tank Level of vehicle k's power type z upon arrival to vertex j τ j Time variable specifying the time of arrival of a vehicle at vertex j, initialized to zero upon departure from the depot f ijkz Vehicle k's consumption from vertex i to vertex j using power type z (electric or fuel)