A modified clustering search based genetic algorithm for the proactive electric vehicle routing problem

In this paper, an electric vehicle routing problem with time windows and under travel time uncertainty (U-EVRW) is addressed. The U-EVRW aims to find the optimal proactive routing plan of the electric vehicles under the travel time uncertainty during the route of the vehicles which is rarely studied in the literature. Furthermore, customer time windows, limited loading capacities and limited battery capacities constraints are also incorporated. A new mixed integer programming (MIP) model is formulated for the proposed U-EVRW. In addition to the commercial CPLEX Optimizer version 20.1.0, a modified Clustering Search based Genetic algorithm (MCSGA) is developed as a solution method. Numerical tests are conducted on the one hand to validate the effectiveness of the proposed MCSGA and on the other hand to analyze the impact of travel time uncertainty of the electric vehicle on the solutions quality.


Introduction
Nowadays, to manage climate change and reduce greenhouse gases emission, the transportation system is forced to operate with the safest fuel sources instead of the usage of petroleum-based fuel.One of the cleanest and most achievable fuel sources is electricity.Indeed, the most important delivery companies have incorporated electric vehicles (EVs) into their fleet.Operating with EVs have many advantages, for example, lower C02 emissions, less maintenance and minimal noise.However, due to the limited battery power, these types of vehicles cannot travel long distances without visiting a charging station, which is time and cost consuming.Therefore, the first concern of delivery companies is to find the optimal routing of EVs which minimizes the time and the cost of deliveries.This problem is a variant of the classical vehicle routing problem (VRP) and is commonly known as the electric vehicle routing problem (EVRP).The EVRP consists in using electric vehicles to complete the last mile delivery operations with a minimum cost route (Toth et Vigo., 2014).For large details on EVRP, please refer to (Qin et al., 2021).Most papers on EVRP addressed the problem in a deterministic environment.However, in the real-world case, unexpected events are usually present.For example, the travel time between two points depends on several factors such as the state of route, the weather, the traffic jump, or the possible accidents.These unexpected factors can increase the estimated travel time.Consequently, a drastic change in routing plans should be applied at the last moment.Hence, planners should anticipate these unexpected events to avoid disruptions in their initial planning.This work addresses how to find a routing plan of the electric vehicle under travel time uncertainty.A proactive strategy that considers the occurrence probability of a finite set of discrete scenarios that can happen during the route planning is proposed.The contributions of this paper are a novel mixed integer mathematical modeling for the U-EVRW and a modified clustering search based genetic algorithm with a Multi-Parent order Crossover operator (MCSGA) to solve it.The rest of the paper is structured as follows: The next section presents a literature review on the EVRP.Then, a mathematical model for U-EVRW is presented in Section 3. In Section 4, a detailed description of the proposed MCSGA is presented.Next, computational experiments and comparison results are exposed in Section 5. Finally, Section 6 concludes the paper.

Related works
Papers dedicated to logistics activities are increasingly a fertile area of operations research, with the vehicle routing problem (VRP) being one of the most studied and well-known transport problems that consider the increasing use of EVs (Afroditi et al., 2014).In general, studies on EVRP aim to provide decision support tools and to encourage the switch to green energy, including electric vehicles.Over the years, many interesting extensions of VRP have been proposed by including different real-world constraints.First, (Erdoğan & Miller-Hooks, 2012) proposed a mixed integer linear program for the Green Vehicle Routing Problem (GVRP) and solved this model by a modified Clarke and Wright savings heuristic and a Density-Based Clustering Algorithm.Then, Schneider et al. (2014) presented the study of the Electric Vehicle Routing Problem with Time Windows (E-VRPTW).They developed a hybrid heuristic that combines a variable neighborhood search algorithm with a tabu search heuristic to solve the problem.Later, the EVRP is studied by different authors.Nowadays EVRP is one of the most important extensions of VRP, due to its environmental benefits.This literature review will progressively highlight articles on VRP published since the year 2014 that reviewed aspects of green logistics and new solution approaches.(Allahviranloo et al., 2014) proposed three new formulations to address the selective vehicle problem under uncertain demand level.Three Parallel genetic algorithms are developed to solve these formulations at large scale level.The results proposed by the authors show that the Parallel genetic algorithm wherein communication between subpopulations, does not eliminate repeated chromosomes, outperforms other algorithms at level of computation time.
The authors in (Keskin & Çatay, 2016) opted for partial charging restriction instead of full charging; they proposed a binary mixed integer linear program and developed an adaptive large neighborhood search (ALNS) algorithm to solve it.In (Hiermann et al., 2019), the authors proposed a mixed integer programming model for the EVRP.They considered a fleet with different types of electric vehicles.For small instances, the authors proposed a branch-and-price algorithm as a solution method, while for large instances; they developed a hybrid ALNS-based heuristic with integrated local search.The work given in (Montoya et al., 2017) is the first which addressed the EVRP with a nonlinear charging function for the battery-charge level.According to the results the authors suggest that neglecting nonlinear charging may lead to infeasible or overly expensive solutions.(Cao et al., 2018) proposed a hybrid algorithm based on an artificial immune system for VRP with interval requests.The authors in (Zhang & Tang, 2009) proposed a hybridization of the ant colony optimization (ACO) with scatter search (SS) as a new metaheuristic solution to solve the VRP.In (Verma, 2018), the authors presented a new variant of the E-VRPTW by allowing the available stations to serve both as a traditional recharging station and as a battery swapping station.Hiermann et al. (2016) addressed the VRP with consideration of three types of vehicles, conventional, Plug-in hybrid, and electric vehicles.They developed a metaheuristic which combines a genetic algorithm with local and large neighborhood search to solve this problem.In (Macrina et al., 2019), the authors developed an iterative local search heuristic to solve the VRP with both electric and internal combustion engine vehicles.To respect environmental measures, limitations on the polluting emissions for the internal combustion engine vehicles are considered.(Cortés-Murcia et al., 2019) proposed a new variant of the electric vehicle routing problem in which an alternative vehicle (bikes, drones, segways, etc.) is used to customer visits while the electric vehicle is at recharging station.In (Koç et al., 2019) the authors address an extended version of the electric vehicle routing problem with nonlinear charging function by considering several companies that shared charging stations to determine the routes of their vehicle fleet.The objective of this paper is to find the best location and technology of the charging stations to minimize the cost of travel time of the vehicles of each company.(Pelletier et al., 2019) propose a robust mixed integer linear program for the electric vehicle routing problem under uncertainties in energy consumption due to several factors such as the weather, road conditions and driver behavior.They developed a two-phase heuristic method based on Large Neighborhood Search to solve larger instances of the problem.In (Zhao & Lu, 2019), the authors proposed a heuristic approach based on the Adaptive Large Neighborhood Search (ALNS) and integer programming to solve VRP taken into account customer time windows, limited loading capacities and limited battery capacities of the vehicles.Wang et al (Wang et al., 2019) extended the multi-depot VRPTW problem by considering shared transport resources to minimize the total carbon emission as well as the operational cost.They developed a hybrid metaheuristic which combines the Clarke and Wright Savings Heuristic Algorithm (CWSHA), the Sweep Algorithm (SwA) and the Multi-Objective Particle Swarm Optimization algorithm as a method of resolution.In (Mao et al., 2020), the authors developed a hybrid algorithm which combine An Improved Ant Colony Optimization Algorithm and a Local Search to solve the EVRP with multiple recharging options, which are partial recharging and battery swapping.(Messaoud, 2021) proposed an extension of the EVRP by considering the Stochastic Travel Times (EVRPSTT).To solve this problem, the author developed a Chance Constrained Programming (CCP) Model, as well as an Improved Large Neighborhood Search (ILS) algorithm and a Monte Carlo Sampling (MCS) procedure.In (Pan & Liu, 2023), the authors proposed a formulation based on a deep reinforcement learning framework to deal with the dynamic characteristics of the route in real-world urban logistics.A deep neural network with a dynamic attention mechanism is used to collect the data on the changes in customers' demands in real-time and a cutting-edge reinforcement learning algorithm is presented for better training the routing process's dynamics and uncertainty.

Mathematical model
This section describes an extension of the mathematical formulation presented in (Schneider et al., 2014) for the E-VRPWT.The proposed model incorporates a proactive strategy which considers the travel time uncertainty of the vehicles.In E-VRPWT a fleet of EVs with fixed loading capacities and limited battery power are used for the package delivery from a single depot to a set of customers.Customers are characterized by their demands, delivery time windows and service durations.To complete the delivery operation, the EVs may need to visit a station to recharge the consumed battery during the route.In (Schneider et al., 2014) model's, the authors considered a deterministic travel time between nodes.However, in real world case, due to many unforeseen factors that can affect the travel time of vehicles between nodes, that quantity is highly uncertain and may have a significant impact on the initial planning.In this work, we propose a proactive model that takes into account the travel time uncertainty of vehicles by the integration of the occurrence probability of unexpected scenarios.To clarify the proposed mathematical model, the following notations are considered:

Assumption
The assumptions represent the real-world applications: • Each customer will be served at most once.
• Each customer will be only visited by one vehicle.
• A station can be visited more than once.
• The total demands of customers visited by one vehicle cannot exceed its load capacity.
• Any vehicle can serve any customer.
• The demands, locations, and time windows of customers are known in advance.
• Routes begin and end at the depot.
• The vehicle speed is assumed the same for all vehicles and for simplification and to have distance=time the speed is assumed equal to 1. • The vehicle's traveling time is considered uncertain.

Sets
V: Set of customers vertices.i = (0, N+1) represents the same depot node, and every route starts at the depot 0 and ends at the depot N + 1. F: Set of recharging stations  ′ : Set of dummy vertices generated to enable several visits to each vertex in the set F.  ′ : Set of recharging stations including the starting depot 0.  ′ : Set of recharging stations including the ending depot N+1. : Set of customers including the starting depot 0.  ′ : Set of customer vertices and recharging stations,  ′ =  ∪ . ′ : Set of customers and recharging stations including the starting depot 0 ( ′ =  ′ ∪ 0). ′ : Set of customers and recharging stations including the ending depot N+1 ( ′ =  ′ ∪  + 1).
= Occurrence probability of each scenario k ∈ w dij =Distance between nodes i and j (we note that dij is the Euclidean distance between nodes i and j ) tij =Estimated travel time between nodes i and j  = Travel time between nodes i and j in scenario  ∈ w (uncertain travel time)

Decision Variables
x : is equal 1 if the route between i and j is traversed 0 otherwise, G : Service start time of node i, Gs : Service start time of node i in scenario k, U : Remaining cargo on arrival at vertex i, Us : Remaining cargo on arrival at vertex i in scenario k, bt : Remaining charge battery on arrival at vertex i, bts : Remaining charge battery on arrival at vertex i in scenario k Objective function subject to x ≤ 1 The main objective of the problem is to determine a proactive routing plan to anticipate unexpected scenarios which can alter the expected travel time between nodes.The objective function (1) minimizes both the cost of the total travel time and the cost of the fleet size used in the determined proactive routing plan.Constraints (2) and (3) ensure that each customer is visited exactly once and handle the connectivity of visits to recharging stations.Constraint (4) guarantees that if a customer or a station is visited by a vehicle this last should leave after the end of their services.Constraint (5) ensures that if a vehicle visits a customer j after a customer i, then the service start time of customer j will be greater than or equal to the service start time of customer i plus the service time spent with customer i plus the travel time between customer i and j.Constraint (6) does the same as Constraint (5) but in unexpected scenarios.Constraint (7) ensures that if a vehicle visits a customer j after a recharging station i, then the service start time of customer j will be greater than or equal to the service start time at recharging station i plus the recharge times for a complete recharge in station i plus the travelled time between i and j.Constraint (8) does the same as Constraint (7) but in unexpected scenarios.Constraint (9) guarantees that each customer is served in its time windows.Constraints (10-13) ensure that each vehicle does not exceed its maximum load capacity.Constraints (14-17) determine the battery level of each vehicle and guarantee that none of them run out of battery.Constraints (18-21) define the type of decision variables ( , G , Gs , U , Us , bt , bts ).

Proposed Modified Clustering Search Based Genetic Algorithm (MCSGA)
Due to the NP-hardness of the EVRP (Schneider et al., 2014), no exact method can solve it for large-sized instances in a reasonable time.For this reason, a Modified Clustering Search Based Genetic Algorithm (MCSGA) is proposed as a new alternative approach to solve it within a short running time.The proposed MCSGA is based on the Clustering Search algorithm (CSA) proposed by (Chaves & Lorena, 2010).In our case the CSA is combined with a genetic algorithm and a new crossover mechanism called the Multi-Parent Order Crossover operator (MPC) to improve exploration capacity of the algorithm.The MCSGA combines three main techniques: a genetic algorithm to generate potential feasible solutions for the problem, a clustering search method to group the solutions into different clusters and find promising regions in the feasible space of the problem, and a local search technique to intensify the search in the promising regions.A cluster (i) is characterized by a center (CTi), a capacity maximal (CAmax) and a number ri that represents the number of times that a local search is applied consecutively on the cluster i. (rmax defines the maximal number of times that a local search can be applied consecutively on the cluster (i) without no improvements in the solution.).(CTi) represents the position of the cluster (i) in the space search while (CAmax) is the number of solutions that can be grouped in any cluster (i).To begin the procedure, the MCSGA generates randomly (m) feasible solutions of the problem.Each generated solution represents the center (CTi) of the (m) different clusters.Then, a genetic algorithm is applied to generate potential new solutions (section 4.1).Next, a clustering process is executed to allocate the generated new solutions in the most adequate cluster (section 4.2).Finally, a local search technique is applied whenever a cluster is considered as a promising region.

Genetic algorithm
To generate good solutions for the MCSGA, a genetic algorithm (Holland, 1975) is developed and used as the main components in our CSGA algorithm.Following, the steps of the proposed GA are briefly described.

Solution representation and decoding process
The encoding of the solution is primordial to develop an efficient meta-heuristic.In the proposed algorithm, a solution is represented by an array of 2×(N+K+1) cells, where N is the number of customers and K is the number of times when the recharging stations are visited.The first and the last cell represent the depot "D" because each vehicle should start and end with the depot, the cells with odd indices may contain an integer number in the interval [1, N [.This integer number represents the customer visited by the vehicle or may contains the recharging station "Si" where i [1, M ] and M is the number of recharging stations.The cells with even indices may contain the depot "D" or the same integer number of the previous cell which means that the vehicle will visit another customer without the need for recharging the power battery.Fig. 1 shows an example of the solution representation:

Fig. 1. Solution representation
To ensure the feasibility of a solution, the following constraints should be respected: • Each customer should be visited in their windows time.
• A customer is visited by a vehicle if the battery capacity of the vehicle is at least enough to arrive and leave the customer i.
• The loading capacity of each vehicle should be respected.
The decoding of the solution presented in Fig. 1 is the following: three vehicles are used to complete the delivery.The first one starts its route at the depot then the customers 2, 3 and the recharging station S1 are visited consecutively and finally the vehicle returns to the depot "D".The second, begin its route at the depot, then the customers 4, 1 and 5 are visited and finally the vehicle return to the depot "D" and the last one, starts its route at the depot "D", then the customers 6, 7 are visited and finally the vehicle return to the depot "D".

Initial population
The initial population plays a primordial role in fast convergence of genetic algorithms towards the optimal solution because each new population is generated based on the component of the current population.Therefore, to start with potential solutions, the initial population of the used GA is randomly generated except one of them, which is based on the strategy to visit the D 2 2 3 3 S1 D 4 4 1 1 5 D 6 6 7 D maximum customers as possible without visiting a recharging station.

Fitness
In GA, Potential individuals are selected, and weak individuals are eliminated according to their fitness values.The EVRP is a minimization problem; thus, the higher fitness value must be the smaller objective function value.

Selection operator (Tournament process)
At each generation, a tournament process is applied to select the most interesting solutions and update the population.In this process, to give all individuals the chance to be selected, three, four, or five individuals are selected randomly from the current population and the best of them is selected to be part of the new population.The process is repeated until the population is completely updated.

Crossover process
The crossover operation has a primary impact on the efficiency and diversity of the algorithm.Two crossover operators have been applied in the proposed approach, a classic multiple crossover point operator and a multi-parent order crossover operator.
• Multiple Crossover point operator: This operator consists in mixing the principal components of two solutions (parents) based on multiple crossover points from the current population.The goal of this operator is the generation of two new solutions (offspring) (Fig. 2).

Fig. 2. Multiple Crossover point operator
• Multi-Parent order Crossover operator (MPC): The MPC consists randomly of a set of solutions of parents instead of two parents and mixed their components to have a single offspring (Sc).The procedure of MPC technique in our case is the followed: First, the number (cd) of candidate parents to mix is already fixed as an input.After that, we generate randomly (cd-1) crossover points in the interval [1, 2×(N+K+1)] where N is the number of customers and K is the number of times when the recharging stations are used.The (cd-1) crossover points divides each parent selected into (cd) parts.
Then, the components of the first part of the first parent are selected as the components of the first part of the offspring solution.Similarly, the components of the second part of the second parent are selected as the components of the second part of the offspring solution and so on until the components of the last part of the last parent which are selected as the components of the last part of the offspring solution (Fig. 3). The offspring solution is a legal solution i.e. the offspring respects all constraints of the problem.So, in this case we accept the new solution. The offspring solutions is illegal solution in this case a correction process is applied to ensure the feasibility of the new generated solutions (Fig. 2 and Fig. 3).
An illegal solution is caused by some of the following cases: • Maybe a customer is visited outside of their windows time, in this case a swapping strategy is applied to place the customer in a feasible window time.• Maybe a customer is visited by a vehicle, but the battery capacity of this latter does not allow him to reach the nearest station, in this case a swapping strategy is applied between the customer and the nearest charging station.

Mutation
To ensure the diversity of solutions, the MCSGA applies a swap mutation operator on the center CTi of each promising cluster.Two components from the center CTi are randomly selected and replaced one by another and vice versa to generate a new solution candidate (Fig. 4).

Clustering process (CP)
The CP consists in comparing each new solution generated by GA to the center CTi of each cluster i.Then, this new solution is assigned to the cluster with the most similar center.At each generation, the cluster with a new solution updates their center CTi to change their position within the search space.The similarity between solutions is measured by the distance between them.In this paper, the distance D(i,j) between two solutions i and j is equivalent to the number of cells with different values (Fig. 5).D(i,j) is a positive number that should increase according to the increasing distance between i and j.The purpose of CP is to determine the cluster that represents the most promising region of search.After each iteration of the MCSGA, the number of solutions grouped in each cluster is calculated and if a cluster reaches this maximal capacity (CAmax).Thus, this cluster can be considered as a potential region to explore.Then, a local search technique is applied on the center CT i of each potential cluster i to improve the solution quality.Finally, a mutation operator is applied on the center CTi of each potential cluster i to avoid possible local optimum.This mutation operator is applied only after the application of the local search technique a maximal number of times (rmax) consecutively and no improvement in solution quality is observed.Ending this CP, returns the GA that will generate a new solution.The MCSGA stops after a number maximal (Imax) of iterations.

Local search
After each iteration, a local search strategy is applied to improve the solution quality.This strategy is based on the following two procedures: • Swap two customers who are in the planned route of a given vehicle.The objective of this swap operation is to change the order in which the customers are visited.

Computational experiments
In this section, two computational experiments studies are presented.The first one evaluates the effectiveness of the proposed MCSGA approach and the second one analyzes the impact of the travel time uncertainty of the vehicles on the quality of solutions.

Parameter Setting Selection
Parameter settings are a primordial factor to improve the performance of metaheuristics.The number of clusters (CL), the size of population (SP), the crossover rate (CR), the number of iteration (Imax) and the number of candidate parents to mixed in multi-parent crossover operation (M) are the five primordial parameters in the proposed MCSGA.To find a suitable combination of the five mentioned parameters, different preliminary experiments based on the Irace package are realized (López-Ibáñez et al., 2016).Irace uses a set of possible values for each parameter and compares the performance of the algorithm on training instances.Table 1 summarizes the results obtained for each parameter tuned in the MCSGA.

Instances
In our experiments, the benchmark instances from (Schneider et al., 2014) are used.These instances are divided into three types: instances with random customer distribution (R), instances with clustered customer distribution (C), and instances with a mixture of both (RC).OV denotes the objective value, AOV denotes the average of the objective value after 30 runs and TR denotes the total run-time in seconds.The maximum duration for CPLEX was set to two hours.
To test the proposed metaheuristic and analyse the impact of travel time uncertainty of vehicles on the quality of solutions, is it necessary to know the optimal value or the upper bound of each instance used in the study.Therefore, just the instances on which we can obtain an optimal or an upper bound with CPLEX version 20.1.0 in a maximum running time of 2 hours are selected.1080 experiences for 36 instances are realized.Where, each of those instances are solved using the proposed MCSGA approach with 30 runs.

Performance of the proposed MCSGA approach on tested instances
In this section, the computational results obtained by the MCSGA algorithm were compared, firstly, with those obtained by CPLEX optimizer for the deterministic model.The deterministic model is the model presented in section 3 but without the constraints that consider the travel time uncertainty of the vehicles and, secondly, with both CPLEX optimizer and the classical CSA algorithm for the under-travel time uncertainty model formulation presented in section 3. To consider travel time uncertainty of vehicles, different scenarios and different probabilities of occurrence of these scenarios are used.In instance with five customers one scenario is used, in instance with ten customers, 2 scenarios are used and in instances with 15 or more customers, 3 scenarios are tested.The probabilities of occurrence considered for each scenario k are Pk=0, Pk=0.5 and Pk=1,(∀  ∈ ).This value of Pk has been selected randomly.The proposed MCSGA is implemented with python language on DELL computer with an Intel core i7, 8GB of ram memory.Table 2 provides the results found with CPLEX and MCSGA for deterministic models.The first two columns present the name and the type of instances respectively, column three contains the objective value and the running time obtained by CPLEX and the last column contains the average objective value from 30 runs and the running time obtained by MCSGA.According to the results presented in Table 2, we conclude that the MCSGA method can obtain the optimal solution for all instances tested after 30 runs.We observe that in all instances tested the optimal solution is achieved in each run among 30 runs tested except in "rc108C15", "rc204C15" and "rc103C15", This is because when the size of customers increases the optimal solution became hard to achieve.At the level of running time, the proposed MCSGA can obtain optimal solutions in a very short running time for all instances while CPLEX take approximately one hour in most of the case and two hours in some tested instances as "c103C15"and "rc108C15".

Table 3
Comparison between the results obtained with CPLEX and MCSGA for under uncertain travel time model OV denotes the objective value, AOV denotes the average of the objective value after 30 runs, TR denotes the total run-time in seconds and Pk denotes the probability of occurrence of scenario k (∀ k ∈w).The maximum duration for CPLEX was set to two hours.

Table 4
Comparison between the results obtained with Classical CSA and MCSGA for under uncertain travel time model OV denotes the objective value, AOV denotes the average of the objective value after 30 runs, TR denotes the total run-time in seconds and Pk denotes the probability of occurrence of scenario k (∀ k ∈w).
Table 3 and Table 4 present the results found with CPLEX, CSA (Chaves & Lorena, 2010) and MCSGA for the under-travel time uncertainty model.In table 3, the first column presents the name of instance, the second column presents the number of scenarios tested in each instance, the column three provide the objective value and the running time obtained by CPLEX for each probability tested (Pk=0, Pk=0.5 and Pk= 1) and the last column provide the objective value and the running time obtained by the MCSGA for each probability tested (Pk =0, Pk =0.5 and Pk =1).In Table 4, the first column presents the name of instance, the second column presents the number of scenarios tested in each instance, the column three provide the objective value and the running time obtained by the CSA for each probability tested (Pk =0, Pk =0.5 and Pk =1) and the last column presents the objective value and the running time obtained by the MCSGA for each probability tested (Pk =0, Pk =0.5 and Pk =1).
Table 3 and 4 show that both, the proposed MCSGA metaheuristic and CPLEX are able to find the optimal solution for the majority of the instances tested whatever the probability of the scenarios.This is not the case for CSA which is unable to find optimal solutions for the instances with 15 customers.The execution times achieved by the MCSGA method are much shorter than the times achieved by CPLEX, especially when 'P=1' and the number of customers is high.According to the tables 2, 3 and 4 we notice that the values of the objective function in under travel time uncertainty model are relatively high in comparison with those of deterministic model.However, with U-EVRWT model we obtain a more flexible route planning that avoids any disruptions in initial planning which is not possible in deterministic model.

Analyzing the impact of travel time uncertainty of vehicles on results
In this section, we analyze the impact of the travel time uncertainty of vehicles on three important components of the final solution which are the number of vehicles used in each routing plan, the objective function values and the running time.Fig. 6, Fig. 7 and Fig. 8 show the impact of travel time uncertainty on: the number of vehicles used in each routing plan, the objective function values and the running time consecutively.According to the figures, we notice that the average number of vehicles, the running time of the algorithm and the objective function values in U-EVRWT model is relatively high in comparison with the average number of vehicles, the running time of the algorithm and the objective function values in deterministic model, consequently, we conclude that, the general cost of vehicles routing plan in U-EVRWT model can be high in comparison with deterministic one, this is correct especially when the number of customers is high.However, we should note that this conclusion is valid just in absence of unexpected events.In cases where the unexpected events and travel time uncertainty are present, the deterministic model returns infeasible solutions while the U-EVRWT returns feasible and efficient solutions.

Impact of travel time uncertainty on the objective function values
c5 C10 C15

Conclusion
In this work an electric vehicle routing problem with time windows and under travel time uncertainty (U-EVRW) is addressed.
A mixed integer mathematical model was proposed to solve the problem and find the optimal proactive routing plan of the electric vehicles under the travel time uncertainty.A MCSGA was developed as a new effective metaheuristic to solve the U-EVRW.Comparison results have indicated that the proposed MCSGA outperform the CPLEX optimizer and CSA in both, solution quality and running time.In addition, the computational experiences have shown the impact of the travel time uncertainty on the average number of vehicles, the running time of the algorithm and the objective function values in U-EVRW model.In the future working on U-EVRW a reactive strategy to solve the problem can be developed.Another possible case of study is the development of an exact solution approach to solve the problem.

:
Set of customers and recharging stations including starting and ending depot ( , ′ =  ′ ∪ 0 ∪ { + 1}.W: Set of scenarios.Parameters LC = Loading capacity of EVs Q = Battery capacity of EVs CostV= Company usage cost of one vehicle (is assumed equal to 50$) CostT= Cost of travel time per minute (is assumed equal to 0.15 $) g= Recharging rate h= Charge consumption rate

Fig. 3 .
Fig. 3. Multi-Parent order Crossover operator (MPC) Correction process After generation of the offspring solution, the following cases can happen:

•
Swap two customers who are in the planning route of two different vehicles.With this operation the planning route of each vehicle can be altered and generate new possible solutions.The general MCSGA is described in Algorithm 1.

Table 1
Instances generated by Irace

Table 2
Comparison between the results obtained with CPLEX and MCSGA for the deterministic model