Keywords: Electric vehicle routing problems, Technician routing problems, Green vehicle routing problems, Fixed-route vehicle charging problems, Non linear charging functions, Metaheuristic, Matheuristic, Iterated local search

eVRP underlying graph G = (V, A). V : Set of vertices of G. 0: Depot I: Subset of V representing customers. F : Subset of V representing CSs. F : Set that contain the set F and β copies of each CS. A: Set of arcs in G. B: Set of breakpoints of the piecewise linear charging function. p i : Service time at customer i ∈ I. c ik : Charging time of the CS i ∈ F at the breakpoint k ∈ B. a ik : Charge level of the CS i ∈ F at the breakpoint k ∈ B.

e ij : Energy consumption between vertices i and j (i, j ∈ V ). t ij : Travel time between vertices i and j (i, j ∈ V ). T max : Tour duration constraint. Q: Tank capacity.

x ij : Binary variable, equal to 1 if an EV travels from vertex i to j and 0 otherwise. τ j : Depart time of the EV from vertex j ∈ V . y j : Charge level when the EV departs from vertex j ∈ V . q i : Charge level when an EV arrives at CS i ∈ F . o i : Charge level when an EV departs from CS i ∈ F

Acknowledgement

This research project has been possible with the help and support of many people. I would like to express my sincere gratitude to all people who were very helpful during my Ph.D. experience. This experience gave me the opportunity to know different countries and cultures, to meet a lot of friendly people, to do new friends, and to improve my academic and research skills.

I first thank to my wife, Diana, for her support, sacrifice, company, and love. She was always there, no matter the place, the weather or the language. I am fortunate to have her by my side. I would also like to thank my mother, Angela, my stepfather Orlando and my grandmother Blanca, for their company. They have supported me all the way in the development of this thesis and they always believed I could follow this dream.

I like to thank my thesis advisors, Christelle Guéret for her patience, support and advice in all this process; Jorge Mendoza I have learned from his many things that go far beyond optimization, vehicle routing problems, and operations research; and Juan Guillermo Villegas for his intellectual generosity and to be with me in great part of my academic formation (my master and Ph.D.).

I would like to thank the jury members Daniele Vigo, Dominique Feillet, Emmanuel Néron, and Fabien Lehuédé for accepting the invitation to be part of the jury.

I want to thank the University EAFIT for the support in this project, especially to Gabriel Arango, Director of Teaching; Alberto Rodriguez, Dean of the Faculty of Engineering; and Sergio Augusto Ramirez, Chief of Production Engineering department for their trust and support. I also want to thank Mario Velez for believing in me since I was doing my career; Jairo Maya for share his time and knowledge; and Gabriel Hincapié and Nora Cadavid for his support and helping. And I would like to thank the Universidad EAFIT scientific computing center (APOLO) for their support in the computational experiments I thank to people of Université d'Angers, especially to Michel Landron for his support and for introducing me to French culture; Simone Rees for her support and help; and my partners at LARIS, Achraf, Fally, Khadim, Ibrahim, Amin and Khaoula for their friendship.

I would like to thank everyone who supported me during my doctorate in Angers. Pedro and Natalia my "paisas brothers", who offered me their company, and supported me in the first two years. I also would thank Alexis, Evelin, Cristhope, Juan Pablo, Clemence, Ana, Silvan and Victoria for the unforgettable moments.

I want to thank the people who supported me in Colombia. To my in-laws and unconditional friends (Joa, Marc, James, Naty, David and Naty) for being with Diana when I was in France; Pauline for her French classes in Medellin; and Pablo Maya for his network flow course at the University of Antioquia.

I like to thank VeRoLog group for allowing me participating to the summer schools of 2014 and 2015. In those schools, I could improve my research skills and made a lot of friends.

Finally, I gratefully acknowledge the financial support provided by Programa de Movilidad Doctoral hacia Francia (Colfuturo -Emb. de Francia -ASCUN -Colciencias -Min. de Educación), and Programa Enlaza Mundos (Alcaldía de Medellín).

Chapter 1 Introduction

Electric vehicles (EVs), which include hybrids, plug-in hybrids, and battery electric vehicles, are one of the most promising technologies to reduce petroleum dependency and greenhouse gas emissions. Boosted by government-sponsored programs worldwide, the EVs segment has become a small but important part of the global automotive industry. Quoting a recent market study published by [START_REF] Berman | Selecting electric vehicles for fleets[END_REF], "sales of plug-in EVs will grow at a compound annual growth rate of nearly 40 percent over the remainder of the decade, while the overall auto market will expand by only two percent a year".

Freight and personal transportation is one of the segments in which the use of EVs is expected to generate higher economical and environmental impact. Studies suggest that transportation represents between one and two thirds of the total costs of logistics in a company [START_REF] Tseng | The role of transportation in logistics chain[END_REF]) and it is responsible for up to 15% of the sale price of goods [START_REF] Avella | Solving a fuel delivery problem by heuristic and exact approaches[END_REF]); a high portion of this cost is linked to vehicle acquisition and operating costs. While the acquisition costs of EVs (vehicles + charging infrastructure) remain significantly higher than those of traditional combustion vehicles (CVs), their operation costs are sensibly lower. Indeed, on average the operating cost of an EV is around 0.11e/km, which is much lower than the 0.22e/km of the average operation cost of a CV1 [START_REF] Bulk | Costs of the electric car[END_REF]). After analyzing a handful of studies comparing the whole-life cost (or vehicle's total cost of ownership) of EVs and CVs, Pelletier et al. (2016) conclude that in some contexts (e.g., urban distribution, low speeds, and frequent stops) EVs area financially viable alternative to CVs. With the advances on technology and the massive adoption of EVs their economical benefits are expected to rise in the years to come. Nonetheless, in the short-to-mid term the highest benefits derived from using EVs in freight transportation are the environmental returns. As a matter of fact, estimates by Energy Information Administration (2008) suggest that transportation is responsible for 33% of the total CO2 emissions in the United States (US), and [START_REF] Pedersen | Optimization models and solution methods for intermodal transportation[END_REF] estimates this figure to be around 28% in the European Union (EU). In addition to economical and environmental aspects, a massive implementation of EVs in both freight and personal transportation is a key strategic issue in terms of energy security for the EU.

The use of EVs in freight and personal transportation is starting to get momentum. Because of the economic benefits and the environmental regulations, several companies in different sectors have started to use EVs in their operations. The courier sector is probably the leader in EV adoption. For instance, since 2011, La Poste operates at least 250 EVs and signed orders for an additional 10,000 [START_REF] Kleindorfer | Fleet renewal with electric vehicles at La Poste[END_REF]. In 2011, UPS purchased 130 hybrid electric delivery trucks (HEVs) and added them to its fleet of 250 HEVs (UPS 2008). Likewise, the Portuguese postal company CTT operates with 10 small electric vans [START_REF] Post | CTT group invests 5M of euros in green fleet[END_REF].

In the service sector, several companies have started EVs adoption programs of EVs. For instance, the largest french electricity company ENEDIS has an ambitious fleet replacement program, they estimate to have 2,000 EVs between 2015 and 2016, which represents the 10% of their current fleet (ERDF 2014). The spanish electricity company Endesa expects to introduce 3,000 EVs in its fleet by 2020, which will represent 80% of its fleet (H&E 2016). In Colombia, a health care company, Sura-EPS, uses EVs for the home healthcare service, in which physicians and nurses visit their patients [START_REF] Loaiza | Convenio Sura-Celsia, por el bien del medio ambiente[END_REF].

Other interesting cases can be found in the freight transportation. In 2013 Coca-Cola deployed more than 30 EVs in selected cities across the U.S [START_REF] Priselac | A frigid fleet: Coca-Cola unveils first electric refrigerated trucks[END_REF]). Heineken have started to use the Europe's largest electric truck to distribute the beer within Amsterdam [START_REF] Heineken | Europe's largest electric truck will drive down emissions[END_REF]. In Mexico, the largest food company, Bimbo, signed orders to buy 100 EVs which will use the energy from a wind farm [START_REF] Ramirez | Bimbo inagura el segundo centro de ventas ecologico en mexico[END_REF].

Despite these encouraging implementations, the massive adoption of EVs for freight transportation is still hampered by technical constraints such as low driving ranges and long battery charging times. Indeed, on average an EV has 120 km of autonomy and its charging time ranges between 0.5 and 12 hours (depending on the charging technology). This autonomy is much lower than the 500 km of average autonomy of a CV, and the EV charging times are considerably longer than CV refueling times.

To overcome these constraints, considerable research has been, and still is, devoted via technological advances in battery capacity, electric efficiency, and battery charging speed. However, as pointed out by [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF], a considerable effort is still needed to develop fleet management tools able to cope with these restricting EV features. One of the fields in which the void is more critical is that of optimization techniques for solving the electric vehicle routing problems (eVRPs).

VRPs are concerned with the design of efficient routes to deliver goods and services from (to) central depots to (from) customer locations satisfying specific business constraints. From an economic perspective, the application of operations research (OR) techniques to solve VRPs has led to savings ranging from 5% to 20% in the global transportation cost in the US and the EU [START_REF] Toth | An overview of vehicle routing problems, in 'The vehicle routing problem[END_REF]. Similarly, the environmental benefits of efficient routing have been remarkable. In particular, a case study reported by [START_REF] Hall | On the road to mobility-2008 survey of vehicle routing software spotlights integration with portable phones[END_REF] points out that Paragon Routing identified in a United Kingdom brewery annual reductions of more than 3,700 tons of CO2 emissions resulting from more efficient routing and depot allocation.

In the last 50 years, a vast amount of research has been devoted to solve VRP variants and it is fair to say that excellent solution approaches exist for most of them (for a comprehensive overview the reader is referred to [START_REF] Toth | Vehicle routing: problems, methods, and applications[END_REF]). Unfortunately, most of those approaches are difficult, if not impossible, to use in the context in which routes are performed by EVs. Indeed, classical VRP variants consider that routes are performed by CVs. These vehicles have a long and easy-to-restore driving range (because petrol stations are available almost everywhere and tank refueling takes a negligible time). As a result, most classical routing algorithms focus on designing only one (long) route per vehicle and do not need to care for scheduling visits to refueling stations. In contrast, due to the short driving range and long battery recharging times of EVs and the limited availability of charging infrastructure (at least for 2016 standards), routing algorithms for eVRPs need to consider the detours of the EVs to the charging stations (CSs) and the decision of the amount of charge at each visited CS. In addition, the algorithms should take into account as much as possible EV features as: the nonlinear charging function, the compatibility constraint between the EVs and CS and the battery degradation cost, among others. To address these challenges, in recent years, the vehicle routing community has started to develop techniques specially tailored for eVRPs.

To the best of our knowledge, the pioneering work in the field of eVRPs is the recharging vehicle routing problem (RVRP) by [START_REF] Conrad | The recharging vehicle routing problem[END_REF]. The RVRP extends the distance-constrained VRP by assuming that routes are serviced by EVs. In this problem, the EVs recharge their batteries at some of the customers. In 2012, [START_REF] Erdogan | A green vehicle routing problem[END_REF] introduced the green vehicle routing problem (Green VRP). The Green VRP is the first problem where the EVs detour to the CSs to recharge their batteries. To solve the problem they proposed a constructive heuristic. That paper is probably the turning point in the study of eVRPs. Since then, different authors have devoted effort not only to propose more sophisticated methods for the Green VRP (Andelmin & Bartolini 2016, Koc & Karaoglan 2016, [START_REF] Montoya | A multi-space sampling heuristic for the green vehicle routing problem[END_REF], Schneider et al. 2015, 2014) but also to study extensions, mostly focusing on EV applications. For instance, Goeke et al. (2015), [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF], [START_REF] Keskin | Partial recharge strategies for the electric vehicle routing problem with time windows[END_REF] and [START_REF] Desaulniers | Exact algorithms for electric vehicle-routing problems with time windows[END_REF] extended the problem to consider customer time windows. [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF] studied a version of the problem where multiple charging technologies are available and partial battery charges are allowed. In their problem, each CSs has a subset of available technologies and the decision maker should plan the routes so a function of the travel, charging and battery degradation costs is minimized. [START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF] tackled a slightly different problem in which along with the routes, the decision maker must decide the size and composition of the EV fleet. To compose the fleet, the decision maker chooses between vehicles of different capacity, battery size, and acquisition cost. Goeke & Schneider (2015) work on a VRP with a fleet composed by EVs and CVs. Finally, [START_REF] Sassi | Iterated tabu search for the mix fleet vehicle routing problem with heterogenous electric vehicles[END_REF] studied the first rich eVRP. In the context of a joint project with La Poste, they tackled a realistic version of the problem where: CSs have different technologies, the fleet is made up of CVs and EVs of different types, there are compatibility constraints between the type of EV and the charging technology of the stations, charging costs depend on the time of the day, and partial charges are allowed.

Despite recent studies in the eVRP field, there is still the need to explore new solution approaches and to study new eVRPs that are closer to reality. Therefore, the purpose of this thesis is to contribute to eVRP field by proposing new methods for current eVRPs, and by studying new eVRPs. Following that objective, each chapter of this manuscript tackles a different eVRP. In consequence, each chapter is self contained in terms of the notation and literature review. A brief description of the dissertation structure follows.

In Chapter 2, we present a simple yet effective two-phase heuristic to tackle the Green VRP. In the first phase our heuristic builds a pool of routes via a set of randomized route-first cluster-second heuristics. In the second phase our approach assembles a Green VRP solution by solving a set partitioning formulation over the columns (routes) stored in the pool. At the core of our first phase lies a novel repair mechanism that optimally inserts visits to stations to restore the feasibility of routes violating the vehicle's autonomy constraint. This procedure solves a constrained shortest path problem on an auxiliary repair graph, that by construction includes only the trips between stations that respect the vehicle's autonomy constraint. To test our approach, we performed experiments on a set of instances from the literature.

The Green VRP and the most of the existing eVRPs in the literature assume that the EVs fully charge their batteries every time they reach a CS, and that the battery charge level is a linear function of the charging time. In practical situations, however, the amount of charge is a decision variable, and the battery charge level is a concave function of the charging time. In Chapter 3, we extend current eVRP models to consider partial charging and non-linear charging functions. We presented a computational study comparing our assumptions with those commonly made in the literature.

In Chapter 4, we propose an iterated local search (ILS) enhanced with heuristic concentration (HC) to tackle the eVRP with non-linear charging function (eVRP-NL). The ILS component uses a variable neighborhood descent (VND) procedure for the local search phase. Considering the importance of the charging decisions (where and how much to charge) in eVRPs, we present a particular local search operator for the charging decisions in the VND. This operator consists in solving a new fixed-route problem: the fixed-route vehicle-charging problem (FRVCP). To solve the FRVCP, we propose a heuristic and a mixedinteger linear programming formulation. The HC component assembles the final solution from the set of all routes found in the local optima reached by the ILS. To test our approach, we propose a new set of instances.

In order to tackle a real problem related with an operation with EVs, Chapter 5 introduces the technician routing problem with conventional and electric vehicles (TRP-CEV), which is inspired in the real operation of ENEDIS. This problem consists in routing a set of technicians to serve a set of geographically scattered requests within their time windows. The technicians use a fixed fleet composed of CVs and EVs. This problem considers the working schedule, lunch break and skills of the technicians. The TRP-CEV involves not only customer sequencing (i.e., routing) and charging decisions but also vehicle-to-technician assignments. The objective function in the TRP-CEV seeks to minimize the total cost, defined as the sum of the travel costs, battery charging fixed costs, the parking cost at CSs, and the fixed cost of using each technician. To tackle the TRP-CEV, we proposed a two-phase parallel matheuristic (PMa). In the first phase, the PMa builds a pool of feasible TRP-CEV routes solving, in parallel, a set of sub-problems. In the second phase the PMa assembles a TRP-CEV solution by solving an extended set covering formulation. Furthermore, we build a set of real instances based on data provided by ENEDIS. We analysed the solutions delivered by our PMa on those instances, and evaluated the quality and structure of the solution as a function of the percentage of EVs in the fleet.

Finally, Chapter 6 presents a general conclusion and future research directions. Appendices describe the detailed results of the algorithms proposed throughout the thesis. The work of this chapter was presented at the following conferences and seminars:

-Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2014). A route-first cluster-second heuristic for the Green Vehicle Routing Problem. In ROADEF 2014, Bordeaux, France.

-Montoya, A., Guéret, C., Mendoza, J. E., and

Introduction

The use of alternative fuel vehicles (AFVs) (e.g., electric vehicles, hydrogen vehicle, among others) in freight transportation leads to new optimization problems. One of these problems is the green vehicle routing problem (Green VRP)1 , introduced by Erdogan & Miller-Hooks (2012). The Green VRP is an extension of the well-known vehicle routing problem, arising when a fleet of AFVs based at a central depot services a set of geographically spread customers. The special feature of this VRP comes from the limited range of AFVs. To ensure the feasible completion of trips, the AFVs may visit alternative fuel stations (AFSs) en-route to refill the tank or recharge the battery.

Formally, the Green VRP is defined on an undirected and complete graph G = (V, E). The vertex set V = {0} ∪ I ∪ F = {0, 1, 2, ..., n + a} is made up of a depot (vertex 0), a set of customers I = {1, 2, ..., n}, and a set F = {n + 1, n + 2, ..., n + a} of a AFSs. It is assumed that the depot can also be used as a refueling station, and that all refueling stations can handle an unlimited number of vehicles. Each vertex i ∈ V has a service time τ i . If i ∈ I then τ i is the service time at the customer; and if i ∈ F ∪ {0} then τ i is the refueling time, which is assumed to be constant. The set E = {(i, j) : i, j ∈ V, i = j} corresponds to edges connecting vertices of V . Each edge (i, j) has two associated nonnegative attributes: a travel time t ij and a distance d ij . The travel speeds are assumed to be constant over the edges. In addition, there is no limit on the number of stops that can be made for refueling. When refueling occurs, it is assumed that the tank is filled to its maximum capacity. The customers are served using a fleet of homogeneous AFVs with tank capacity Q and consumption rate cr. The vehicle driving-range constraint is dictated by the fuel tank capacity and a tour duration constraint T max .

In the Green VRP the objective is to find a set of routes of minimum total distance such that each customer is visited exactly once; the level of the tank when the vehicle arrives at any vertex is nonnegative; each route satisfies the maximum-duration limit; and each route starts and ends at the depot. Figure 2.1 depicts a feasible solution to a Green VRP.

The Green VRP is an NP-hard problem. Indeed, [START_REF] Lenstra | Complexity of vehicle routing and scheduling problems[END_REF] showed that the classical VRP is NP-hard. Since the VRP is a special case of the Green VRP, we can conclude that the Green VRP is also NP-hard. Moreover, recent studies show that commercial solvers cannot solve to optimality instances of 20 customers in reasonable computational times [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF]. Therefore, to tackle industrial-scale Green VRP instances we need heuristic approaches. For the solution of the Green VRP we present a multi-space sampling heuristic (MSH), which is a simple yet effective heuristic introduced by [START_REF] Mendoza | A multi-space sampling heuristic for the vehicle routing problem with stochastic demands[END_REF]. The idea behind MSH is to sample different solution representation spaces and then to assemble a solution with (parts of) the sampled elements. Implementations of MSH have delivered competitive results to complex routing problems such as the VRP with stochastic demands [START_REF] Mendoza | A multi-space sampling heuristic for the vehicle routing problem with stochastic demands[END_REF]), the VRP with stochastic travel and service times [START_REF] Gómez | On modeling stochastic travel and service times in vehicle routing[END_REF], and the combined maintenance and routing problem [START_REF] Fontecha | Combining maintenance and routing models: the case of a water & sewage utility[END_REF]. The algorithm is built out of two main components: a set of sampling functions and an assembling procedure. The sampling functions are randomized route-first, cluster-second heuristics. Using these heuristics MSH draws a sample from the TSP solution representation space and extracts from it a sample of the route representation space. Later, MSH uses the sampled routes to assemble a final solution. The assembling procedure is a set partitioning model that runs over the set of routes sampled in the first phase. To implement our MSH, we adapted the randomized route-first, cluster-second heuristics proposed in [START_REF] Mendoza | A multi-space sampling heuristic for the vehicle routing problem with stochastic demands[END_REF] to the Green VRP. This adaptation is far from being trivial, since route feasibility in the Green VRP is difficult to assess. Indeed, a route that is fuel-infeasible (that is, infeasible for fuel autonomy) can be repaired in a number of ways by inserting one or more visits to AFSs. Therefore, extracting Green VRP routes from a giant TSP tour is much more complex than in problems previously tackled using MSH. Existing metaheuristics for the Green VRP rely on insertion-based heuristics and neighborhood schemes for repairing fuel-infeasible routes [START_REF] Erdogan | A green vehicle routing problem[END_REF], Schneider et al. 2014, 2015[START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF]. In general, these strategies consider only one insertion at a time and therefore are exposed to myopic choices. We propose an optimal procedure based on a reformulation of the repair problem as a constrained shortest path problem.

The main contributions of this chapter are threefold: i) we introduce an optimal repair procedure based on a constrained shortest path formulation that inserts refueling stations into Green VRP routes. One of the main advantages of our procedure is that it can be used as a building block in any other Green VRP heuristic; ii) we show how to use our procedure to build a simple and effective MSH heuristic for the Green VRP; and iii) we update best-known solutions (BKSs) to 8 out of 52 standard benchmark instances for the problem.

The remainder of this chapter is organized as follows. Section 2.2 reviews the related literature. Section 2.3 presents a detailed description of our version of the multi-space sampling heuristic. Section 2.4 presents a computational evaluation of the proposed method. Finally, Section 2.5 concludes the chapter. For easy reference, Section 2.6 summarizes the notation of this chapter.

Literature review

The Green VRP was introduced by [START_REF] Erdogan | A green vehicle routing problem[END_REF]. They propose a mixed-integerlinear programming (MILP) formulation and two heuristics. The first heuristic is a modified Clarke and Wright savings algorithm (MCWS) that repairs infeasible routes by inserting AFSs using a savings criterion, and removes redundant AFSs after merging the routes. The second heuristic is a density-based clustering algorithm (DBCA) that first builds clusters and then runs MCWS for each cluster. In their article, the authors propose two sets of test instances: the set of "small" instances has 40 test problems with 20 customers, and the set of "large" instances has 12 instances with 111 to 500 customers. Computational experiments showed that CPLEX 11.2 was unable to solve to optimality even small instances using the MILP model. They also indicated that there were no significant differences in performance between MCWS and DBCA.

The Green VRP is closely related to the classical distance-constrained VRP [START_REF] Laporte | Optimal routing under capacity and distance restrictions[END_REF], but the latter does not consider the possibility of extending the vehicle's distance limitation. Another problem that is closely related to the Green VRP is the multi-depot vehicle routing problem with inter-depot routes (MDVRPI) described by [START_REF] Crevier | The multi-depot vehicle routing problem with inter-depot routes[END_REF]; this problem considers intermediate depots at which vehicles can be replenished with goods. To solve this problem [START_REF] Crevier | The multi-depot vehicle routing problem with inter-depot routes[END_REF] presented a heuristic procedure that combines ideas from adaptive memory programming, tabu search, and integer programming. Tarantilis to avoid confusions with the already established and more-studied generalized vehicle routing problem (GVRP). et al. (2008) proposed a hybrid guided local search heuristic that outperforms the procedure proposed by [START_REF] Crevier | The multi-depot vehicle routing problem with inter-depot routes[END_REF]. [START_REF] Muter | A branch-and-price algorithm for the multidepot vehicle routing problem with interdepot routes[END_REF] proposed a branch and price algorithm for the MDVRPI; it solves to optimality some instances with up to 50 customers. [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF] introduced the electric vehicle routing problem with time windows and recharging stations (E-VRPTW), which is an extension of the VRP with a fleet of EVs. The E-VRPTW considers limited vehicle freight capacities, customer time windows, and the possibility of recharging at any of the available stations using an appropriate recharging scheme. [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF] presented an MILP formulation and a hybrid metaheuristic combining variable neighborhood search and tabu search (VNS/TS). Their VNS/TS explores infeasible solutions with respect to capacity, time windows, and battery-usage constraints. A dynamic penalizing scheme is used to guide the search toward feasible solutions. The VNS component explores 15 neighborhoods based on cyclic exchanges [START_REF] Thompson | Cyclic transfer algorithm for multivehicle routing and scheduling problems[END_REF] 2015) propose an adaptive variable neighborhood search (AVNS). Their AVNS uses a modified savings algorithm to generate an initial solution that is later improved with local search. The algorithm uses an adaptive shaking with twentyfour neighborhood structures, five route selection methods, three vertex-sequence selection methods, and an adaptive mechanism to choose the route and vertex-selection methods. The solution generated at the shaking step is subsequently improved by several greedy local searches. Furthermore, the AVNS has a dynamic penalization scheme to guide the search toward feasible solutions and a simulated annealing acceptance criterion. Since the Green VRP is a special case of the VRPIS, [START_REF] Schneider | An adaptive VNS algorithm for vehicle routing problems with intermediate stops[END_REF] tested their approach on the instances of [START_REF] Erdogan | A green vehicle routing problem[END_REF]. This method outperformed all previous methods both in terms of solution quality and computational time.

Recently, [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF] introduced the green vehicle routing problem with multiple technology and partial recharges (GVRP-MTPR). As the name suggests, in this problem charging stations may have different technologies (e.g., charging times) and EVs do not necessarily charge their batteries to the full capacity when they reach a charging point. In their article, the authors presented an MILP formulation, a local search method (48A), and a simulated annealing (SA). Their local search method uses 48 possible combinations of 6 different neighborhoods and selects the best overall solution. Their SA uses a relocate neighborhood to explore the solution space. Every time that the incumbent solution is updated, the SA uses a deterministic local search to try to further improve the solution. Since the GVRP-MTPR is an extension of the Green VRP, (2015). A plausible explanation for these results, is that their 48A algorithm and their SA are not specifically tailored to the Green VRP. [START_REF] Sassi | Vehicle routing problem with mixed fleet of conventional and heterogenous electric vehicles and time dependent charging costs[END_REF] tackled an electric vehicle routing problem with heterogeneous mixed fleet and time dependent charging costs. In their problem the authors consider a number of realistic features such as: different charging technologies, coupling constraints between vehicles and charging technologies, charging station availability time windows, and charging costs depending on the time of the day. Regarding charging station capacity, the authors impose the maximum admissible power charge constraint over the charging spots located at the depot; however, they assume that the other charging stations have an unlimited capacity. They proposed an MILP formulation, a charging routing heuristic and a local search heuristic. They performed computational experiments on real data instances.

Gonçalves et al. (2011) considered a VRP with pick-up and delivery (VRPPD) with a mixed fleet that consists of battery electric vehicles (BEVs) and vehicles with internal-combustion engines. The objective is to minimize the total costs (vehicle-related fixed and variable costs). They considered time and capacity constraints and assumed a time for recharging the BEVs, which they calculated from the total distance traveled and the range for one battery charge. However, they did not incorporate the location of the recharging stations into their model. Thus, they basically proposed a mixed-fleet VRPPD with an additional distance-dependent time variable. They performed computational experiments with an MILP formulation on instances with up to 17 customers. [START_REF] Conrad | The recharging vehicle routing problem[END_REF] introduced the recharging vehicle routing problem (RVRP) wherein vehicles with a limited range must service a set of customers but may recharge at certain customer locations before continuing their trip. They proposed an MILP formulation of the problem. They performed computational experiments on instances of 40 customers. [START_REF] Juan | Routing fleets with multiple driving ranges: Is it possible to use greener fleet configurations?[END_REF] discussed the vehicle routing problem with multiple driving ranges (VRPMDR), an extension of the classical routing problem where the total distance each vehicle can travel is limited and is not necessarily the same for every vehicle. The VRPMDR finds applications in routing electric and hybridelectric vehicles, which can cover limited distances depending on the running time of their batteries. They proposed an MILP formulation and a multi-round heuristic algorithm that iteratively constructs a solution for the problem.

Finally, [START_REF] Pelletier | Battery electric vehicles for goods distribution: a survey of vehicle technology, market penetration, incentives and practices[END_REF] presented an overview of the field of goods distribution with EVs, that includes a review of the main transportation science literature on EVs regarding fleet size and mix, vehicle routing problem, and optimal paths.

Multi-space sampling heuristic

Mendoza & Villegas (2013) originally proposed the multi-space sampling heuristic for the vehicle routing problem with stochastic demands. Despite its simple design, MSH obtains competitive results. MSH has two phases: sampling and assembling. In the sampling phase the algorithm uses a set of randomized TSP heuristics to draw a biased sample from the set K of TSP-like tours (i.e., giant tours visiting all customers). Following the route-first cluster-second principle [START_REF] Beasley | Route first-cluster second methods for vehicle routing[END_REF][START_REF] Prins | Order-first split-second methods for vehicle routing problems: A review[END_REF], MSH extracts every feasible route that can be obtained without altering the order of the customers of each sampled TSP tour. MSH uses these routes to build a set Ω ⊂ R, where R is the set of all feasible routes. In the assembling phase MSH follows the principle of petal heuristics [START_REF] Foster | An integer programming approach to the vehicle scheduling problem[END_REF]) and maps set Ω to a solution s (in the set of all feasible solutions to the problem S) by solving a set partitioning formulation. Note that in MSH the knowledge of the problem is embedded in two components: the procedure that extracts routes from the TSP tours during the sampling phase and the set-partitioning formulation used in the assembling phase. The former controls the feasibility and cost of each route while the latter controls the feasibility and cost of the whole solution. To adapt MSH to the Green VRP we designed a tailored route extraction procedure. To favor scalability we also modified the strategy used by the original MSH to build the route pool Ω. For the sake of completeness, in the remainder of this section we present all the components used in our MSH setting a special focus on those we specifically designed for the Green VRP.

General structure

Algorithm 1 describes the general structure of our MSH. The procedure starts by entering the sampling phase (lines [5][6][7][8][9][10][11][12][13][14][15][16][START_REF] Liao | The electric vehicle touring problem[END_REF][18][19]. At each iteration k ≤ K, the algorithm selects a sampling heuristic from a set H (line 6) and uses it to build a TSP tour tsp. Then, the algorithm uses a tour splitting procedure (known as split) to retrieve a solution s ∈ S. Differently to the original MSH, our version does not store in Ω all the routes evaluated by split during the partitioning process but only the routes belonging to some of the retrieved solutions. To decide if the routes of a solution s should join Ω we use the following condition (lines 9-16): f (s) ≤ f (s * ) • (1 + λ) where s * is the best solution found, λ a positive parameter, and f (•) denotes the objective function of a solution. The idea behind this choice is to favor computational scalability by reducing the size of Ω while assuring a good compromise between the diversity and the quality of the routes in the pool. In the assembly phase (line 20), the heuristic invokes a procedure called SetPartitioning to solve a set partitioning formulation over Ω using f (s * ) as an upper bound. The resulting solution R is reported by the heuristic (line 21).

Algorithm 1 Multi-space sampling heuristic: General structure

1: function MULTISPACESAMPLING(G, H, K, λ) 2: Ω ←-∅ 3: k ←-1 4:
while k ≤ K do 5:

for j = 1 to j = |H| do 6: h ←-H j 7: tsp ←-h(G) 8: s ←-split(G, tsp) 9:
if k = 1 and j = 1 then 10:

s * ← s 11: else if f (s) ≤ f (s * ) × (1 + λ) then 12: Ω ←-Ω ∪ s 13: if f (s) < f (s * )

Sampling heuristic

To sample K (line 7 in Algorithm 1), our approach uses randomized versions of three TSP constructive heuristics. Although the strategies used to generate the randomized versions of the three heuristics are directly borrowed from [START_REF] Mendoza | A multi-space sampling heuristic for the vehicle routing problem with stochastic demands[END_REF], for the sake of completeness we briefly describe them here.

Let tsp be an ordered set representing the TSP customer tour being built by a given sampling heuristic, W the set of customers visited by tsp, and N = I \ W an ordered set of nonrouted customers. For the sake of simplicity, we assume that sets W and N are updated every time a customer is added to tsp. Let l be a random integer in {1, ..., min{L, |N |}}, where parameter L denotes the randomization factor of each heuristic.

The first heuristic is randomized nearest neighbor (RNN). It initially sets tsp = (0) and u = 0. At each iteration, RNN identifies the customer v that is the lth nearest customer to u, appends v to tsp, and sets u = v. RNN stops when |N | = 0 and appends 0 to tsp to complete the tour. The second heuristic is randomized nearest insertion (RNI). It initializes tsp as a tour starting at the depot and performing a round trip to a randomly selected customer. At each iteration, RNI sorts N in nondecreasing order of d min (v), where d min (v) is defined to be min{d u,v |u ∈ W}. RNI then inserts v = N l (i.e., the lth element in the ordered set N ) into the tour tsp in the best possible position (i.e., the position generating the smallest increment in the cost of the tour). RNI stops when |N | = 0. The third heuristic is randomized best insertion (RBI). It initializes tsp as a tour starting at the depot and performing a round trip to a randomly selected customer. At each iteration RBI sorts N in nondecreasing order of ∆ min (v), where ∆ min (v) is defined to be min{d(u, v) + d(v, w) -d(u, w)|(u, w) ∈ tsp}, and inserts v = N l in tour tsp in the best possible position. RBI stops when |N | = 0.

Split

To extract a feasible solution s from tsp (line 8 in Algorithm 1), our approach uses an adaptation of the optimal tour splitting procedure for the VRP introduced by [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]. The splitting procedure builds a directed acyclic graph G * = (V * , A) composed of the ordered vertex set V * = (v 0 , v 1 , ..., v i , ..., v n ) and the arc set A. Vertex v 0 = 0 is an auxiliary vertex, while vertices v 1 , ..., v n ∈ tsp \ {0}. The vertices in V * are arranged in the same order in which they appear in tsp. Arc (v i , v i+nr ) ∈ A represents a feasible route r v i ,v i+nr with distance d rv i ,v i+nr starting and ending at the depot and visiting customers in the sequence v i+1 to v i+nr . Route r v i ,v i+nr may not satisfy the fuel constraint (i.e., the route's fuel consumption is greater than Q). If it does not, we try to repair it by inserting visits to AFSs. If the insertion of AFSs increases the duration of the route beyond T max then we do not include the arc associated with the route in G * . The insertion of visits to AFSs is accomplished by the optimal repair procedure explained in Section 2.3.4. To obtain a feasible solution s for the Green VRP, the procedure finds the set of arcs (i.e., routes) along the shortest path connecting 0 and v n in G * . It is worth noting that since G * is directed and acyclic, building the graph and finding the shortest path can be done simultaneously [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF].

Algorithm 2 shows the tour splitting procedure based on the split algorithm for the capacitated VRP [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]). After initializing the shortest path labels (lines 2-6), split enters the outer loop (lines 7-41). Each pass through the outer loop sets the tail of an arc and initializes an empty route (lines 8-11). Then we use the inner loop (lines 13-40) to scan all the arcs sharing the same tail node. At each inner-loop iteration, we explore a new arc by simply adding the next customer to the route (line 16). In the next steps (lines 17-23) we compute the weight and time of the arc (i.e., the total distance c r and total time tp r of the associated route). If an arc is feasible for the duration constraint (line 24) we check the feasibility of the associated route with respect to the fuel constraint (and store the result in a boolean variable f r ). When a route is infeasible for the fuel constraint but tp r ≤ T max , we try to repair the route using procedure Repair(r,G) (line 27). For an arc (i, j) to be added to the graph, its corresponding route must be feasible for fuel and time. If the shortest path label of the head node of arc (i, j) can be improved (i.e., label i-1 + c r ≤ label j ), then we update the shortest path label and the predecessor information (lines 35-38). The algorithm then moves to the next inner-loop iteration or exits the loop. After completing the outer loop we retrieve the solution using the tour tsp and the predecessor labels (for an algorithm to retrieve the solution we refer the reader to [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]).

To adapt the tour splitting procedure to the Green VRP, it is necessary to introduce two important functions: checkFuel (line 24), which evaluates the feasibility of a route with respect to the vehicle's fuel constraint, and Repair (line 26), which tries to repair the route, and returns a boolean variable (indicating whether or not the route could be repaired) and the distance of the repaired route. Note, however, that there are some cases where the repaired route is feasible for the fuel constraint but infeasible for the duration constraint. This occurs because the time needed to recover the feasibility of the route (i.e., the travel time to the AFSs and the service time there) increases its planned duration beyond T max ; in these cases we do not add the associated arc to the graph.

Figure 2.2 illustrates the tour splitting procedure. Figure 2.2a shows the TSP tour and the relevant information. Figure 2.2b depicts the auxiliary graph G * ; the arcs in bold correspond to the shortest path. Figure 2.2c is a table showing the evaluation of the arcs.

Repair procedure

Existing approaches for the Green VRP depend on insertion-based heuristics and neighborhood schemes for repairing fuel-infeasible routes. The MCWS algorithm proposed by Erdogan for i = 1 to n do 4:

label i ←-+∞ 5:
pred i ←-0 pred: predecessor labels 6: end for 7:

for i = 1 to n do continue ←true 13:

while j ≤ n and continue=true do 14:

v ←-tsp j Get customer in position j of tsp 15:

v ←-tsp j-1 Get customer in position j -1 of tsp 16:

r ←-r ∪ {v} 17: return s 44: end function starts inserting the AFS with the least insertion cost (i.e., distance added to the tour) in the fuel-infeasible back-and-forth routes. After the routes merging step, MCWS inserts into each fuel-infeasible route the AFS with the least insertion cost and then evaluates the possibility of eliminating redundant AFSs. Schneider et al. (2014) used a penalization scheme for the fuel-infeasible routes and proposed an operator that performs insertions and removals of AFSs. That operator uses a tabu list in order to control insertions already tested. [START_REF] Schneider | An adaptive VNS algorithm for vehicle routing problems with intermediate stops[END_REF] used a penalization scheme for the fuel-infeasible routes and presented three operators to improve the location of AFSs within the route. The first operator moves an AFS to a different position within a route. The second operator evaluates for each AFS visit of each route whether visiting a different AFS decreases the routing cost. The third operator aims at removing redundant AFS visits. Finally, Felipe et al. (2014) proposed a constructive heuristic, which considers the insertion of AFSs to fuel-infeasible routes. They presented a neigboorhood scheme to improve the location of an AFS in a route, without modifying the sequence of visits to the customers. Notice that all these strategies consider only one insertion at a time and therefore are exposed to myopic choices. In contrast, our repair procedure optimally chooses which stations to insert and where to insert them while considering all the possible combinations of insertions leading to an energy-feasible route. To accomplish its goal, our procedure relies on a reformulation of the repair problem as a constraint shortest path problem.

if j = i then 18: c r ←-d 0,v + d v,0 Update total distance 19: tp r ←-t 0,v + t v,0 + τ v + τ 0 Update total time 20: else 21: c r ←-c r -d v ,0 + d v ,v + d v
Let Π = {π 1 , ..., π i , ..., π j , ..., π nr+1 } be a route that violates the fuel constraint, and π 1 = 0. The feasibility of Π may be restored by inserting visits to AFSs. We use an optimal procedure that simultaneously decides which stations must be visited and their optimal insertion position within the route. The procedure can be seen as a constrained shortest path problem in a repair graph B = (Z, U ). Figure 2.3 depicts the structure of B2 . The vertex set Z = {α, ..., [i, k], ..., β} is made up of two dummy vertices (α and β) that act as copies of the depot, representing the source and sink vertices of B; and the vertices [i, k], which represent a visit to station k ∈ F after visiting vertex π i (i.e., the ith element in route Π), where 1 ≤ i ≤ n r + 1 and 1 ≤ k ≤ a. To define the edge set U , let us first introduce some key elements. Let P = {p 1 , ..., p j , ..., p |P | } be a path in G. For a given path P we define three metrics: its distance d(P ) = |P |-1 j=1 d p j ,p j+1 , its total planned time t(P ) = |P |-1 j=1 t p j ,p j+1 + |P | j=2 τ p j , and its fuel consumption q(P ) = cr × d(P ). The arc set U is composed of five types of arcs, that is, U = 5 i=1 U i , where:

-U 1 : the outgoing arcs of α. An arc (α, [i, k]) represents the path P = (0, F k ) if i = 1, and P = (0, ..., π i , F k ) if i > 1. Its cost and time are defined as c u:u∈U 1 = d(P ) and t u:u∈U 1 = t(P ) + τ 0 . -U 2 : the arcs connecting two stations without visiting any customer. Arc ([i, k], [i, l]) represents the path P = (F k , F l ). Its cost and time are defined as c u:u∈U 2 = d(P ) and t u:u∈U 2 = t(P ). -U 3 , the arcs connecting two stations and visiting some customers in between. Arc ([i, k], [j, l])

represents the path P = (F k , π i+1 , ..., π j , F l ). Its cost and time are defined as c u:u∈U 3 = d(P ) and t u:u∈U 3 = t(P ). -U 4 , the incoming arcs to β representing the return to the depot after visiting some customers since the last visit to a station. Arc ([i, k], β) represents the path P = (F k , π i+1 , ..., π nr+1 , 0). Its cost and time are defined as c u:u∈U 4 = d(P ) and t u:u∈U 4 = t(P ) -τ 0 . -U 5 , the incoming arcs to β representing a return to the depot directly from a station after visiting the last customer in route Π. The arc ([n r + 1, k], β) represents the path P = (F k , 0). Its cost and time are defined as c u:u∈U 5 = d(P ) and t u:u∈U 5 = t(P ) -τ 0 . We include in U only arcs with c u × cr ≤ Q and t u ≤ T max . With this graph construction procedure we ensure that all paths from α to β represent a route that visits customers in the same order they appear in Π and is feasible with respect to the fuel constraint. Note, however, that not every path connecting α and β in B represents a route that is feasible in terms of the duration constraint. Therefore, to find an optimal repair for route Π, we need to solve a constrained shortest path problem (CSP), where the maximum travel time T max is the constrained resource. It is interesting to observe that since all the arcs of B are feasible for fuel, the tank capacity is not a constraint for the CSP. To solve the CSP we use the pulse algorithm (Lozano & Medaglia 2013a). The algorithm is based on the idea of propagating pulses through a network from a source node to a sink node (α and β in our case). At the core of the algorithm lies the ability to effectively and aggressively prune pulses (i.e., prevent their propagation) without jeopardizing the optimal path. The pulse algorithm is one of the state-of-the-art methods for the solution of resource-constrained shortest path problems (Lozano & Medaglia 2013b). Figure 3.3 shows a detailed example of the repair procedure. Figure 2.4a illustrates an infeasible route with three customers and three AFSs for a Green VRP. Figure 2.4b shows the route after the optimal insertion of the AFSs using the repair procedure. Figure 2.4c shows the corresponding repair graph and the shortest path.

Set partitioning

In the assembly phase, MSH maps the set Ω to a solution in S by solving a set partitioning formulation min R⊆Ω r∈R d r : ∪ r∈R = V ; r i ∩ r j = 0 ∀r i , r j ∈ R . The objective is then to select the best subset of routes from Ω to build the set of routes R (i.e., the final solution) guaranteeing that each customer will be visited by exactly one route.

Computational experiments

To test our approach, we ran it on the 52-instance testbed proposed by [START_REF] Erdogan | A green vehicle routing problem[END_REF]. These instances consist of 5 sets. Four sets contain ten 20-customer instances ranging from 2 to 10 AFSs. The remaining set, resulting from a case study, consists of 12 instances with the number of customers ranging from 111 to 500 and between 21 and 28 AFSs. In this testbed there are instances with infeasible customers. Therefore 2014) filter the customers that either cannot be served directly within the maximum route duration or whose direct service requires visiting more than one refueling station. To allow comparison with previously published results, we followed the same convention and filter unfeasible customers. It is important to remark that in the problem definition an AFS vertex could have as successor another AFS, so it is possible to insert visits to more than one AFS between two customer visits. According to Erdogan & Miller-Hooks (2012), the geographical coordinates given for the customers have to be converted to Cartesian coordinates using the Haversine formula [START_REF] Bullard | Plane and Spherical Trigonometry: With Stereographic Projections[END_REF] with an average Earth radius of 4,182.45 miles. It is worth noting that although [START_REF] Erdogan | A green vehicle routing problem[END_REF] include an additional constraint on the number of vehicles in their Green VRP formulation, their experiments are conducted with an unlimitedfleet version of the problem. Therefore, our results are directly comparable to theirs.

We implemented our MSH in Java (jre V.1.7.0_51) and used Gurobi Optimizer (version 5.6.0) to solve the set partitioning problem. We set a time limit of 10 • n seconds on Gurobi to control the running time of the set partitioning problem. All the experiments were run on a computing cluster with 2.33 GHz Intel Xeon E5410 processors with 16 GB of RAM running under Linux Rocks 6.1.1. Each replication of the experiments was run on a single processor.

Results

After conducting a parameter tuning campaign, we set L RN N , L RN I , and L RBI to 2, and λ = 1. We found that these parameters lead to a pool of well-diversified routes. For the sake of brevity, we will not discuss these experiments.

Tables 2.1 and 2.2 summarize the results delivered by our MSH on the small and large Green VRP instances running with 3 configurations: K = 1, 000; K = 5, 000; and K = 10, 000 (defined as MSH(K)). We compare our results to the best result obtained by MCWS and DBCA by Erdogan & Miller-Hooks (2012), the VNS/TS by [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF], the AVNS by [START_REF] Schneider | An adaptive VNS algorithm for vehicle routing problems with intermediate stops[END_REF], and the 48A and SA by [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF]. The rows of Tables 2.1 and 2.2 indicate the number of times that each method found the BKS, the average and maximum gap between the average solution and the BKS (in %), the average best gap 3 (in %), the cumulative number of vehicles 4 , and the average and maximum computational times (in minutes). The results of VNS/TS, AVNS and our MSH are computed over 10 runs and the 48A and SA over a single run. [START_REF] Erdogan | A green vehicle routing problem[END_REF] reported the best solution of several runs with different parameters but did not give the exact number of runs. The detailed results are reported in Appendix A.1.

For small instances, MSH(5k) and MSH(10k) have competitive results with reference to MCWS/ DBCA, VNS/TS, AVNS, and 48A (Table 2.1). They have an average gap of only 0.04% and 0.01%, match the 40 BKSs, and use the lowest cumulative number of vehicles. In terms of CPU time, MSH(5k) and MSH(10k) are only outperformed by MSH(1k) and 48A. MSH(1k) is the fastest method, and outperforms MCWS/DBCA, VNS/TS, and 48A in terms of solution quality.

For large instances, MSH(5k) outperforms MCWSA/ DCBA, VNS/TS, 48A, and SA in terms of both solution quality (i.e., average best gap, and cumulative number of vehicles) and CPU time (Table 2.2). When compared to AVNS, MSH(10k) has a better average best gap with respect to BKSs (0.05% vs. 0.17%) and finds more BKSs (8/12 vs. 4/12), whereas AVNS is faster and it seems to scale better.

It is interesting to observe that for both the small and large instances our three MSH configurations use the lowest cumulative number of vehicles.

In order to give a graphical view of key metrics presented in Tables 2.1 and 2.2, Figure 2.5 presents a head-to-head comparison between the solution methods. It shows the trade-off between solution quality and CPU time 5 . Each method is represented as a point in the plot. The X coordinate of the point represents the average best gap with respect to the BKS, while the Y coordinate represents the average CPU time.

Notice that MCWS and DBCA were not included in the comparison since CPU times for these algorithms were not reported in by [START_REF] Erdogan | A green vehicle routing problem[END_REF].

For the small instances, Figure 2.5a shows that MSH(5k) and MSH(10k) dominate the AVNS and VNS/TS, while MSH(1k) dominates 48A and VNS/TS. Moreover, MSH(5k) has a remarkable performance on small instances, because its CPU time is close to that of 48A (which is the fastest approach in the literature) and its average best gap is 0.00%. On the other hand, on the large instances, MSH(5k) and MSH(10k) dominates SA, 48A and VNS/TS (Figure 2.5b). The comparison between MSH(10k) and AVNS shows that there is not clear dominator: while the former dominates in solution quality, the latter dominates in CPU time. Nonetheless, it is worth recalling that while AVNS uses elaborate components (i.e., an adaptive shaking with twenty-four neighborhood structures, five route selection methods, among others), our MSH uses simple building blocks that, with the notable exception of the repair procedure, are common in the literature. Therefore our MSH is probably easier to implement and extend to tackle other problems.

Finally, we analyzed the CPU time of the two MSH phases (sampling and assembly). Figure 2.6 presents the percentage of time spent on each phase for each instance size (in terms of number of customers). The show that, for all instances, assembly is the more time consuming phase, and that the percentage of time spent on the assembly phase increases with the instance size. However, the assembly phase contributes to improve the solutions on average by 4.10%, which is key to obtain competitive results with respect to existing approaches in the literature.

Conclusion

This chapter proposes a multi-space sampling heuristic for the Green VRP. This approach has three main components: a set of three randomized TSP heuristics, a tour partitioning procedure, and a set partitioning formulation. At the core of the tour partitioning procedure lies a novel repair mechanism that optimally inserts visits to refueling stations to restore the feasibility of routes violating the vehicle's fuel constraint. This procedure solves a constrained shortest path problem on an auxiliary repair graph, that by construction includes only the trips between refueling stations that respect the vehicle's fuel constraint.

The procedure is based on a reformulation of the repair problem as a constrained shortest path problem where i) the maximum travel time of a route is the constrained resource, and ii) each path from the source node to the sink node in the underlying graph models a possible way to repair the route by inserting visits to refuelling stations.

We tested our approach on a 52-instance public testbed for the Green VRP. Our approach found 8 new BKSs for the testbed and matched another 40. When compared to state-of-the-art metaheuristics, the multispace sampling heuristic obtains competitive results in terms of solution quality and computational time, and it is one of the simplest methods used to solve the Green VRP.

Nomenclature

τ i : Service time at a vertex. If i ∈ I then τ i is service time at customer; if i ∈ F ∪ {0} then τ i is refueling time at AFS. Q: Tank capacity.

Notation for multi-space sampling heuristic

Introduction

As mentioned in the Chapter 1, because of the short driving range, eVRP solutions frequently include routes with planned detours to charging stations (CSs) where the electric vehicles (EVs) recharge. In general, eVRP models make assumptions about the capacity of the CSs, the EV energy consumption, and the EV battery charging process. This chapter focuses on the latter. The objective is to evaluate the impact of different charging assumptions commonly used in the literature on the quality and feasibility of eVRP solutions.

The remainder of this chapter is organized as follows. Section 3.2 reviews the charging assumption commonly used in the literature. Section 3.3 presents a mixed-integer linear programming (MILP) formulation for a particular eVRP. With mild modifications, the model can be adapted to work with different charging assumptions. Section 3.4 presents a computational experiment to compare the battery charging assumptions. Finally, Section 3.5 concludes the chapter. For easy reference, Section 3.6 summarizes the notation of this chapter.

Charging asumption in the literature

To model the battery charging process, eVRP models make assumptions about the charging policy and the charging function approximation. The former defines how much of the battery capacity can be (or must be) restored when an EV visits a CS, and the latter models the relationship between battery charging time and charging level. With respect to the charging policies, the eVRP literature can be classified into two groups: studies assuming full and partial charging policies. As the name suggests, in full charging policies, the battery capacity is fully restored every time an EV reaches a CS. On the other hand, in partial charging policies, the amount of charge (and thus the time spent at each charging point) is a decision variable.

In general, the charging functions are nonlinear, because the terminal voltage and current change during the charging process. This process is divided into two phases. In the first phase, the charging current is held constant, and thus the state of charge (SOC) increases linearly with time until the battery's terminal voltage increases to a specific maximum value (see Figure 3.1). In the second phase, the current decreases exponentially and the terminal voltage is held constant to avoid battery damages. The SOC increases then concavely with the time [START_REF] Pelletier | Goods distribution with electric vehicles: Battery degradation and behaviour modeling[END_REF]. Although the shape of the charging functions are known, their exact modeling is very complex because it depends on many factors as: current, voltage, self-recovery and temperature, among others [START_REF] Wang | Estimation of state of charge of batteries for electric vehicles[END_REF]. The battery state of charge is then often described by differential equations. Since such equations are difficult to integrate to optimization models for transportation problems, researchers rely on approximations of the actual charging functions. We discuss these approximations bellow. For each of them, we present a graphical comparison with respect to real charging data provided by [START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF]. These authors conducted experiments to estimate the charging time for different charge levels with two types of EVs and three types of CSs. It is worth noting that each approximation can be used in a full charge policy (FC) or in a partial charge policy (PC). ). This approximation is optimistic, because it assumes that batteries charge up to Q faster than they do in reality. In the second approximation (L2) the charging rate is the slope of the line connecting the first and last observations (see Figure 3.3b) of the charging curve. This approximation tends to be pessimistic, because over a large portion of the curve, the charging rate is slower than in reality. Piecewise linear approximations (PL) : This approximation, proposed by Zündorf (2014) for a shortest path problem with EVs consists in approximating the charging function by a series of linear segments. This approximation is the closest to the real data (see Figure 3.4). To assess the quality of this approximation, we fit piecewise linear functions to the [START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF] data and obtained approximations with an average relative absolute error of 0.90%, 1.24%, and 1.90% for CSs of 11, 22, and 44 kW, respectively.

Battery level

Time

𝑄

Real data Approximation

Battery level

Time

Real data Approximation

Setting up the study: problem and formulations

To accomplish our study about the impact of charging assumptions commonly used in the literature on the quality and feasibility of eVRP solutions, we need two elements: an particular version of an eVRP and a solution method of the problem. Note that the current charging function approximations are particular cases of the piecewise linear. Therefore, if we define an eVRP with piecewise linear charging function and select a method to solve it, we have all the elements to accomplish the study. We thus introduce in this section the eVRP with non-linear charging function (eVRP-NL) (which has never been studied in the literature), and use its MILP formulation and a commercial solver as solution method. The remainder of this subsection describes formally the eVRP-NL and introduces its MILP formulation.

Problem description

Formally, the eVRP-NL can be defined on a directed and complete graph G = (V, A). The vertex set V = {0} ∪ I ∪ F is made up of a depot (vertex 0), a set of customers I, and a set of CSs F . Each customer i ∈ I has a service time p i . Each CS i ∈ F has a nonlinear charging function, which is modeled using a piecewise linear approximation. This approximation is defined by a set of breakpoints B, where each breakpoint k ∈ B is associated to a charging time c ik and a charge level a ik (see Figure 3.6a). The set A = {(i, j) : i, j ∈ V, i = j} corresponds to arcs connecting vertices in V . Each arc (i, j) has two associated nonnegative values: a travel time t ij and an energy consumption e ij . The customers are served using an unlimited and homogeneous fleet of EVs. All the EVs have a battery of capacity Q (expressed in kWh) and a maximum tour duration T max . It is assumed that the EVs leave the depot with a fully charged battery, and that all the CSs can handle an unlimited number of EVs simultaneously. Feasible solutions to this eVRP satisfy the following conditions: each customer is visited exactly once; each route satisfies the maximum-duration limit; each route starts and ends at the depot; and the battery level when an EV arrives to and departs from any vertex is between 0 and Q.

Since the distance is directly related to the energy consumption, most work on eVRPs with homogeneous fleet focuses on minimizing the total distance (Schneider et ). However, this objective function neglects the impact of charging operations in the cost of the solutions. This may lead to decisions such as: charging the batteries more than needed, or charging the batteries when their level is high. These decisions directly affect the battery long-term degradation cost (which can be as high as 3 times the energy cost according to Becker et al. In the charging functions, variables q i and o i specify the battery levels when an EV arrives at and departs from CS i ∈ F . The charging function maps these variables to charging times s i and e i , in order to estimate the time spent at the CS (∆ i ). In this example, Route 1 does not visit any CS, because its total energy consumption is less than the battery capacity. On the other hand Route 2 visits CS 8. In this route, the EV arrives at the CS with a battery level q 8 = 1.0, and it charges the battery to a level o 8 = 6.0. To estimate the time spent at the CS, we use the piecewise linear charging function: the charging time associated to q 8 and o 8 are s i = 0.8 and e i = 6.0, then the time spent at CS 8 is ∆ 8 = 6.0 -0.8 = 5.2. The duration of Route 2 is the sum of the travel time (13.0), the charging time (5.2), and the service time (1.0), that is, 19.2 which is lower than T max. The cost of this route is 18.2 (travel time + charging time). Finally, Route 3 visits CSs 10 and 9; and it spends ∆ 10 = 7.2 and ∆ 9 = 1.6 time units charging in these CSs, respectively.

MILP formulation

To formulate the eVRP-NL, we introduce set F containing the set F and β copies of each CS (i.e., |F | = |F | × (1 + β)). The value of 1 + β corresponds to the number of times that each CS can be visited. For this MILP formulation, we use the following decision variables: variable x ij is equal to 1 if an EV travels from vertex i to j, and 0 otherwise. Variables τ j and y j track the time and charge level when the EV departs from vertex j ∈ V . Variables q i and o i specify the charge levels when an EV arrives at and departs from CS i ∈ F , and s i and d i are the associated charging times (see Figure 3.6b). Variable ∆ i = d i -s i represents the time spent at CS i ∈ F . Variables z ik and w ik are equal to 1 if the charge level is between a i,k-1 and a ik , with k ∈ B \ {0}, when the EV arrives at and departs from CS i ∈ F respectively. Finally, variables α ik and λ ik are the coefficients of the breakpoint k ∈ B in the piecewise linear approximation, when the EV arrives at and departs from CS i ∈ F respectively. The MILP formulation follows:

𝑐 𝑖1 𝑐 𝑖2 𝑐 𝑖3 𝑐 𝑖0 𝑎 𝑖0 𝑎 𝑖1 𝑎 𝑖2 𝑎 𝑖3 = 𝑄

Battery level

Time

t ij x ij + i∈F ∆ i (3.1) subject to j∈V,i =j x ij = 1, ∀i ∈ I (3.2)
j∈V,i =j

x ij ≤ 1, ∀i ∈ F (3.3) j∈V,i =j x ji - j∈V,i =j x ij = 0, ∀i ∈ V (3.4) e ij x ij -(1 -x ij )Q ≤ y i -y j ≤ e ij x ij + (1 -x ij )Q, ∀i ∈ V, ∀j ∈ I (3.5) e ij x ij -(1 -x ij )Q ≤ y i -q j ≤ e ij x ij + (1 -x ij )Q, ∀i ∈ V, ∀j ∈ F (3.6) y i ≥ e i0 x i0 , ∀i ∈ V (3.7) y i = o i , ∀i ∈ F (3.8) y 0 = Q (3.9) q i ≤ o i , ∀i ∈ F (3.10) q i = k∈B α ik a ik , ∀i ∈ F (3.11) s i = k∈B α ik c ik , ∀i ∈ F (3.12) k∈B α ik = k∈B z ik , ∀i ∈ F (3.13) k∈B z ik = j∈V x ij , ∀i ∈ F (3.14) α ik ≤ z ik + z i,k+1 , ∀i ∈ F , ∀k ∈ B \ {b} (3.15) α ib ≤ z ib , ∀i ∈ F (3.16) o i = k∈B λ ik a ik , ∀i ∈ F (3.17) d i = k∈B λ ik c ik , ∀i ∈ F (3.18) k∈B λ ik = k∈B w ik , ∀i ∈ F (3.19) k∈B w ik = j∈V x ij , ∀i ∈ F (3.20) λ ik ≤ w ik + w i,k+1 , ∀i ∈ F, ∀k ∈ B \ {b} (3.21) λ ib ≤ w ib , ∀i ∈ F (3.22) ∆ i = d i -s i , ∀i ∈ F (3.23) τ i + (t ij + p j )x ij -T max (1 -x ij ) ≤ τ j , ∀i ∈ V, ∀j ∈ I (3.24) τ i + ∆ j + t ij x ij -(S max + T max )(1 -x ij ) ≤ τ j , ∀i ∈ V, ∀j ∈ F (3.25) τ j + t j0 ≤ T max , ∀j ∈ V (3.26) τ 0 ≤ T max (3.27) x ij = 0, ∀i, j ∈ F : m ij = 1 (3.28) τ i ≥ τ j , ∀i, j ∈ F : m ij = 1, j ≤ i (3.29) τ j ≤ T max i∈V x ij , ∀j ∈ F (3.30) i∈V x ih ≥ j∈V x jf , ∀h, f ∈ F : m hf = 1, h ≤ f (3.31) x ij ∈ {0, 1}, ∀i, j ∈ V (3.32) τ i ≥ 0, y i ≥ 0 ∀i ∈ V (3.33) z ik ∈ {0, 1}, w ik ∈ {0, 1}, α ik ≥ 0, λ ik ≥ 0, ∀i ∈ F , ∀k ∈ B (3.34) q i ≥ 0, o i ≥ 0, s i ≥ 0, d i ≥ 0, ∆ i ≥ 0, ∀i ∈ F (3.35)
The objective function (4.1) seeks to minimize the total time (travel times plus charging times). Constraints (3.2) ensure that each customer is visited once. Constraints (3.3) ensure that each CS is visited at most once. Constraints (3.4) impose the flow conservation. Constraints (3.5) and (3.6) track the battery charge level at each vertex. Constraints (3.7) ensure that, if the EV travels between a vertex and the depot, it has sufficient energy to reach its destination. Constraints (3.8) reset the battery tracking to o i upon departure from CS i ∈ F . Constraint (3.9) ensures that the battery charge level is Q at the depot. Constraints (3.10) couple the charge levels when an EV arrives at and departs from any CS. Constraints (3.11-3.16) define the charge level (and its corresponding charging time) when an EV arrives at CS i ∈ F (based on the piecewise linear approximation of the charging function). Similarly, constraints (3.17 This MILP formulation can be easily adapted to model the other approximations (FS, L1, and L2) and the different charging policies (FC and PC):

-Full charge and piecewise linear approximations (FC-PL): we replace constraints (3.17) and (3.18) by

o i = k∈B λ ik a ib , ∀i ∈ F (3.36) e i = k∈B λ ik c ib , ∀i ∈ F (3.37)
-Partial charge using the first segment (PC-FS): To run our MILP formulation with this assumption, we modify the input data to include only the first segment. -Partial charge and linear approximations (PC-L): To run our MILP formulation with PC-L1 and PC-L2, we modify the input data so that in the piecewise linear approximation there is a single segment with the corresponding charging rate.

Computational experiments 3.4.1 Experimental settings & Test instances

We conducted an experiment to compare the different charging assumptions commonly used in the literature, combining charging policies and charging function approximation. The evaluated assumption are: partial charging and piecewise linear approximation (PC-PL), full charging and piecewise linear approximation (FC-PL), partial charging and first segment (PC-FS), Partial charge and linear approximations (PC-L1 and PC-L2).

With the aim that the conclusions of the comparison are independent of the solution method and consequently of the quality solution, we just compare optimal solutions of the evaluated assumptions. To ensure the optimal solution we solve a set of small instances.

We ran the experiments on a set of 20 randomly generated instances with 10 customers and 2 or 3 CSs. We located the customers using either a random uniform distribution, a random clustered distribution, or a mixture of both. The CSs are located using two criteria: randomly or using a simple p-median heuristic. To generate the charging functions approximations we use the real charging data for a 16 kWh battery provided by [START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF]. Finally, the maximum route duration for every instance was fixed to 10 hours.

Experimental environment & Parameter setting

As mentioned in the last section, the MILP formulation uses β copies of the CSs to model multiple visits to the same CS. Although several authors followed this strategy ( 2016), they do not explain how the value of β is set. It is worth noting that β plays an important role in the definition of the solution space, and therefore it restricts the optimal solution of the model. For instance, an optimal solution found with β = 3 may not be optimal for β = 4. In practice, there is no restriction on the number of times that a CS can be visited, but large values of β result in models that are computationally intractable. To overcome this difficulty, we designed an iterative procedure to solve the MILP formulation for increasing values of β. Starting with β = 0, at each iteration, our procedure (i) tries to solve the MILP formulation to optimality with a time limit of 100 h, and (ii) sets β = β + 1. The procedure stops when the time limit is reached or an iteration ends with a solution s β satisfying f (s β ) = f (s * β-1 ), where f (•) denotes the objective function and * an optimal solution. We implemented our MILP formulation (and its variant) in Gurobi 5.6

Results

Table 3.1 presents the results. Considering that the PC-PL assumption is the closest assumption to the reality, the results obtained using the others assumptions are compared with reference to the results of the PC-PL. For each charging assumption, we give the objective function value (of ), the percentage gap between of and the PC-PL solution (G), the number of routes in the solution (r), and the value of β. Since in practice the charging time is controlled by the nonlinear charging function, the charging decisions of the PC-L solutions are evaluated a posteriori using the piecewise linear approximation. The last rows of Table 3.1 summarize the results. We present, for each assumption, the average and maximum percentage gap, the number of solutions employing more EVs than in the PC-PL solution, and the number of infeasible solutions.

The results show that solutions based on the full charging policy perform poorly in terms of both objective function (+20.11% on average) and number of routes (8 solutions use a larger fleet) with respect to those based on the partial charging policy. This is because the EVs spend more time than necessary at the CSs. The main motivation for a full charging policy is to avoid complex charging-quantity decisions. However, according to our results, the gain in simplicity does not offset the loss of solution quality.

In the PC-FS assumption EVs can charge their batteries up to only around 80% of the actual capacity. Artificially constraining the capacity may force EVs to detour to CSs more often than necessary when traveling to distant customers. Because the maximum route duration is limited, the time spent detouring and recharging the battery reduces the number of customers that can be visited. Consequently, more routes may be needed to service the same number of customers. Our results confirm this intuition: in 3 out of the 20 instances the PC-FS assumption increases the number of routes. Furthermore, in practice some distant customers may not be included in routes unless the EVs can fully use their battery capacity. In our experiments, 9 instances become infeasible under PC-FS. In conclusion, although PC-FS simplifies the problem (avoiding the nonlinear segment of the charging function) it may lead to solutions that are infeasible, or with larger fleets and (on average) 2.70% more expensive.

As mentioned before, PC-L1 assumes that batteries charge faster than they do in reality (Figure 3.3a). As a consequence, routes based on PC-L1 may in practice need more time to reach the planned charge levels. The extra time may make a route infeasible if there is little slack in the duration constraint. Indeed, the post-hoc evaluation shows that for 14 instances, the PC-L1 solutions are infeasible in practice. On the other hand, PC-L2 assumes that batteries charge slower than in reality (Figure 3.3b). Overestimating the charging times does not lead to feasibility issues, but the resulting routes may be overly conservative. For instance, in our experiments PC-L2 leads to solutions that are (on average) 1.45% more expensive, and it increases the number of routes in 2 instances.

of of of r r r β β β of of of G(%) G(%) G(%) r r r β β β of of of G(%) G(%) G(%) r r r β β β of of of r r r β β β of of of G(%) G(%) G(%) r r r β β β of of of r r r β β β of of of G(%) G(%) G(%) r r

Conclusion

In this chapter, we reviewed the different assumptions used in the literature to model the charging process in eVRPs. In order to know the impact of the charging assumption in the quality and feasibility of the solutions, we conducted a study comparing the optimal solutions for an eVRP under different battery charging assumptions. To conduct this study, we introduce the eVRP-NL, its MILP formulation (which can adapted to model different charging assumption), and a set of small instances with a real information about the charging functions.

Our results show that the full charging policy may lead to overly expensive solutions. Furthermore, they show that the linear charging and first segment approximation have a negative impact on the quality and feasibility of the solutions. In consequence, it is advisable to consider charging function approximations that capture the nonlinear behavior of the process when modelling eVRPs.

Nomenclature

Introduction

In Chapter 3, we showed the importance of considering the partial charging and charging function approximations that capture the nonlinear behavior of the charging process in the eVRPs. However, it is necessary to propose an approach to solve industrial-sized instances of the electric vehicle routing problem with non linear charging function (eVRP-NL). We thus propose a hybrid metaheuristic combining iterated local search (ILS) and heuristic concentration (HC). The ILS component uses a variable neighborhood descent (VND) procedure for the local search phase. The VND uses three local search operators. The first two operators are classical relocate and 2-Opt moves. On the other hand, the third one is a specialized operator for improving the charging decisions. This operator relies on solving the fixed-route vehiclecharging problem (FRVCP). A new problem defining the charging decisions (i.e., where and how much to charge) of an energy-infeasible fixed-route (i.e., a route infeasible with respect to the energy autonomy of the vehicle). To solve the FRVCP, we propose a heuristic and a mixed-integer linear programming (MILP) formulation. Finally, the HC component assembles the final solution from the set of all routes found in the local optima reached by the ILS.

The main contributions of this chapter are threefold. First, we introduce and study a new variant of the fixed-route vehicle-charging problem (FRVCP) that appears as a subproblem of the eVRP-NL. The FRVCP consists in finding the optimal charging decisions (i.e., where and how much to charge) for a route servicing a fixed sequence of customers. Second, we propose a hybrid metaheuristic to solve the eVRP-NL and assess its performance in a new set of benchmark instances. Third, we analyze our solutions and provide some insight into the characteristics of good eVRP-NL solutions.

The remainder of this chapter is organized as follows. Section 4.2 introduces our hybrid metaheuristic. Section 4.3 discusses the FRVCP and presents two approaches to solve it. Section 4.4 presents a computational evaluation of the proposed method. Finally, Section 4.5 concludes the chapter. For easy reference, Section 4.6 summarizes the notation of this chapter.

Hybrid metaheuristic

To solve the eVRP-NL we developed a hybrid metaheuristic combining ILS [START_REF] Lourenço | Iterated local search: Framework and applications[END_REF]) an HC (Rosing & ReVelle 1997). Figure 4.1 presents the general structure of the proposed approach (hereafter referred to as ILS+HC).

To find an initial solution we follow a sequence-first split-second approach which uses a constructive heuristic to build a TSP tour visiting all the customers and a splitting procedure to retrieve an eVRP-NL solution. Then, at each iteration of the ILS we improve the current solution using a variable neighborhood descent (VND) [START_REF] Mladenović | Variable neighborhood search[END_REF] with three local search operators: relocate, 2-Opt, and global charging improvement (GCI). At the end of each ILS iteration, we update the best solution and add the routes of the local optimum to a pool of routes Ω. To diversify the search, we concatenate the routes of the local optimum to build a new TSP tour, and then perturb the new TSP tour. We start a new ILS iteration by splitting the perturbed TSP tour. After K iterations the ILS component stops, and we carry out the HC. In this phase, we solve a set partitioning problem over the set of routes Ω to obtain an eVRP-NL solution. In the remainder of this section, we describe the main components of our method.

Initial solution

We generate the initial TSP tour using the simple and well-known nearest neighbor heuristic (NN). For a description of NN see Rosenkrantz et

Split

To extract a feasible solution from a TSP tour, our approach uses an adaptation of the splitting procedure introduced by [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]. The splitting procedure builds a directed acyclic graph G * = (V * , A * ) composed of the ordered vertex set V * = (v 0 , v 1 , ..., v i , ..., v n ) and the arc set A * . Vertex v 0 = 0 is an auxiliary vertex, and each vertex v i represents the customer in the ith position of the TSP tour. Arc (v i , v i+nr ) ∈ A * represents a feasible route r v i ,v i+nr with an energy consumption e rv i ,v i+nr , starting and ending at the depot and visiting customers in the sequence v i+1 to v i+nr .

Note that since the TSP tour only includes customers, route r v i ,v i+nr may be energy-infeasible (i.e., the total energy needed to cover the route is greater than Q). In that case, we solve a FRVCP to obtain an energy-feasible route by inserting visits to CSs. If visiting CSs increases the duration of the route beyond T max , we do not include the arc associated with the route in A * . Finally, to obtain a feasible eVRP-NL solution, the splitting procedure finds the set of arcs (i.e., routes) along the shortest path connecting 0 and v n in G * .

Variable neighborhood descent

To improve the solution generated by the splitting procedure we use a VND based on three local search operators. The first two operators, namely, relocate and 2-Opt1 , focus on the customer sequencing decisions. In other words, these two operators only alter the sequence of customers and do not have the ability to insert, remove, or change the position of CSs. To update the charging times after a relocate or 2-Opt move we use the rule proposed by Felipe et al. ( 2014): when visiting a CS, charge the strict minimum amount of energy needed to continue the route until reaching the next CS (or the depot if there is no other CS downstream). If reaching the next CS (or the depot) is impossible, even with a fully charged battery, the move is deemed infeasible. Similarly, if after updating the charging times the resulting route is infeasible in terms of the maximum-duration limit, the move is simply discarded. It is worth nothing that the Felipe et al. (2014) rule is optimal when all CSs are homogeneous; nonetheless, that is not the case in our eVRP-NL.

As its name suggest, the third operator, global charging improvement or GCI, focuses on the charging decisions. GCI is applied to every route visiting at least one CS. The operator works as follows. First, GCI removes from the route all visits to CSs. If the resulting route is energy-feasible, the operator stops. On the other hand, if the route is energy-infeasible, GCI solves a FRVCP trying to optimize the charging decisions (where and how much to charge) for the concerned route. Depending on the configuration, our ILS+HC solves the underlying FRVCP either heuristically or optimally. Full details on the FRVCP and the solution techniques embedded in our approach to solve it are given in Section 4.3.

Perturb

To diversify the search our approach concatenates the routes of the current best solution to build a TSP tour. Then, we perturb the resulting TSP tour with a randomized double bridge operator [START_REF] Lourenço | Iterated local search: Framework and applications[END_REF]) and then apply the split procedure to obtain a new eVRP-NL solution.

Heuristic concentration

Finally, the heuristic concentration component solves a set partitioning formulation over the pool of routes Ω. The objective is then to select the best subset of routes from Ω to build the final solution guaranteeing that each customer is visited by exactly one route.

The fixed-route vehicle-charging problem

As mentioned above our ILS+HC relies on solving a variant of the FRVCP at two specific points of its execution: when evaluating routes in the split procedure, and when searching the GCI neighborhood. The FRVCP is a variant of the well-known fixed-route vehicle-refueling problem (FRVRP). The FRVRP seeks the minimum-cost refueling policy (which fuel stations to visit and the refueling quantity at each visited station) for a given origin-destination route [START_REF] Suzuki | A variable-reduction technique for the fixed-route vehicle-refueling problem[END_REF]). Most of the research carried on the FRVRP and its variants applies only to internal combustion vehicles (which have negligible refueling times). Nonetheless, a few FRVCP variants have received attention in the literature. Most of these variants assume full charging policies (Montoya et 2014) assume a partial charging policy. Their problem differs from ours in three fundamental ways: i) they do not take into account the charging times (because their objective is to minimize the energy and degradation costs), ii) they do not deal with maximum route duration constraints, and iii) in their problem the CSs are already included in the fixed route and no detours are to be planned. In the remainder of this section we introduce our FRVCP and the techniques embedded in our approach to solve it. For the sake of simplicity, in the remainder of the manuscript to refer to our FRVCP variant simply as the FRVCP.

Let Π = {π(0), π(1), ..., π(i), ..., π(j), ..., π(n r )} be an energy-infeasible route, where π(0) and π(n r ) represent the depot. The route has a total time t, which is the sum of the travel times plus the service times. The feasibility of Π may be restored by inserting visits to CSs. As mentioned in Section 3.3, each CS j ∈ F has a piecewise linear charging function defined by a set of breakpoints B. The piecewise linear function is composed by a set of segments. Each segment is defined between the breakpoints k -1 and k ∈ B, it has a slope ρ jk (representing a charging rate), and it is bounded between the battery levels a jk-1 and a jk (see Figure 4.2a). Figure 4.2b shows the fixed-route Π and the possible visits to the CSs between two vertices of Π. In Figure 4.2b the values e π(i-1)π(i) and t π(i-1)π(i) represent the energy consumption and the travel time between vertices π(i -1) and π(i) ∈ Π. Similarly, e π(i-1)j and t π(i-1)j represent the energy consumption and the travel time between vertex π(i -1) ∈ Π and CS j ∈ F , and e jπ(i) and t jπ(i) represent the energy consumption and the travel time between the CS j ∈ F and vertex π(i) ∈ Π.

In the FRVCP the objective is to find the charging decisions (where and how much to charge) that minimize the sum of the charging times and detour times while satisfying the following conditions: the level of the battery when the EV arrives at any vertex is nonnegative; the charge in the battery does not exceed its capacity; and the route satisfies the maximum-duration limit. Since the FRVRP is NP-hard [START_REF] Suzuki | A variable-reduction technique for the fixed-route vehicle-refueling problem[END_REF]) and the FRVCP generalizes the FRVRP, we can conclude that the FRVCP is also NP-hard.

Mixed-integer linear programming formulation

We formulate the FRVCP using the following decision variables: variable ε π(i)j is equal to 1 if the EV charges at CS j ∈ F before visiting vertex π(i) ∈ Π. Variable φ π(i) tracks the battery level. If ε π(i)j = 0, φ π(i) is the battery level when the EV arrives at vertex π(i). On the other hand, if ε π(i)j = 1, φ π(i) is the battery level when the EV arrives at CS j ∈ F right before visiting vertex π(i). Variable θ π(i)jk is equal to 1 if the EV charges on the segment defined by breakpoints k -1 and k ∈ B at CS j ∈ F before visiting vertex π(i) ∈ Π. Finally, variables δ π(i)jk and µ π(i)jk are (respectively) the amount of energy charged and the battery level when the charging finishes on the segment between breakpoints k -1 and k ∈ B at CS j ∈ F before the visit to vertex π(i) ∈ Π. The MILP formulation of the FRVCP follows:

min π(i)∈Π\{π(0)} j∈F k∈B\{0} δ π(i)jk ρ jk + π(i)∈Π\{π(0)} j∈F ε π(i)j (t π(i-1)j + t jπ(i) -t π(i-1)π(i) ) (4.1) Subject to φ π(1) = Q - j∈F ε π(1)j e π(0)j -e π(0)π(1)   1 - j∈F ε π(1)j   ( 4
φ π(i) = φ π(i-1) + j∈F k∈B\{0} δ π(i-1)jk - j∈F ε π(i-1)j e jπ(i-1) - j∈F ε π(i)j e π(i-1)j -e π(i-1)π(i)   1 - j∈F ε π(i)j   ∀π(i) ∈ Π \ {π(0), π(1), π(n r )} (4.
µ π(i)j1 = φ π(i) + δ π(i)j1 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (4.6) µ π(nr)j1 = φ π(nr-1) + l∈F k∈B\{0}
δ π(nr-1)lk -l∈F e lπ(nr-1) ε π(nr-1)l -ε π(nr)j e π(nr-1)j + δ π(nr)j1 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (4.7)

µ π(i)jk = µ π(i)j,k-1 + δ π(i)jk ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0, 1} (4.8) µ π(i)jk ≥ a jk-1 θ π(i)jk ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0, 1} (4.9) µ π(i)jk ≤ a jk θ π(i)jk + (1 -θ π(i)jk )Q ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (4.10) j∈F ε π(i)j ≤ 1, ∀π(i) ∈ Π \ {π(0)} (4.11) θ π(i)jk ≤ ε π(i)j ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (4.12) δ π(i)jk ≤ θ π(i)jk Q ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (4.13) t + π(i)∈Π\{π(0)} j∈F k∈B\{0} δ π(i)jk ρ jk + π(i)∈Π j∈F ε π(i)j (t π(i-1)j + t jπ(i) -t π(i-1)π(i) ) ≤ T max (4.14) φ π(i) ≥ 0, ∀π i ∈ Π \ {π(0)} (4.15) ε π(i)j ∈ {0, 1}, ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (4.16) θ π(i)jk ∈ {0, 1} ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (4.17) δ π(i)jk ≥ 0 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (4.18) µ π(i)jk ≥ 0 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (4.19)
The objective function (4.1) seeks to minimize the total route time (including charging and detour times). Constraints (4.2-4.5) define the battery level when the EV arrives at vertex π(i) ∈ Π if ε π(i)j = 0; or to CS j ∈ F before visiting vertex π(i) ∈ Π, if ε π(i)j = 1. Constraints (4.6-4.8) define the battery level when the EV finishes charging at CS j ∈ F in the segment between breakpoints k -1 and k ∈ B before visiting vertex π(i) ∈ Π. Constraints (4.9-4.10) ensure that if the EV charges on a given segment, the battery level lays between the values of its corresponding break points (a j,k-1 and a jk ). Constraints (4.11) state that only one CS is visited between any two vertices of the fixed route. Constraints (4.12) ensure that the EV only uses segments of visited CSs. Likewise, constraints (4.13) ensure that the EV charges only at selected segments of visited CSs. Constraint (4.14) represents the duration constraint of the route. Finally, constraints (4.15-4.19) define the domain of the decision variables.

Solving the FRVCP

To solve the FRVCP, we embedded into our method two different solution approaches. The first consists in solving the MILP introduced in the previous subsection using a commercial solver. The second consist in solving the problem using a greedy heuristic adapted from the literature. The remainder of this subsection describes these approaches.

Approach 1: commercial solver Because of its verbosity, at first glance the model introduced in Section 4.3.1 may seem too complex to be efficiently solved using out-of-the-box software. Nonetheless, since in practice the number customers per route and the number of available CSs tend to be low, the resulting MILP formulations are within the scope of commercial solvers. For instance, on an problem with 8 CSs, the MILP formulation for a fixed route serving 10 customers has: 539 continuous variables, 352 integer variables, and 1,263 constraints 2 . We optimally solved that model using Gurobi Optimizer (version 5.6.0) in 0.06s. Based on this observation we decided to embed into our ILS+HC a component that relies on a commercial solver to tackle the FRVCP.

To further reduce the size of the MILP formulation and consequently improve solver's performance, we propose four preprocessing strategies that eliminate infeasible CSs insertions. Our strategies rely on the two following premises: (i) the energy consumption and the travel time between vertices satisfy the triangular inequality; and (ii) since the piecewise linear charging function is concave (i.e., ρ j,k-1 ≥ ρ jk ), the first segment has the fastest charging rate.

We propose two types of strategies. The first three strategies filter CS insertions that are infeasible independently of how the customers are packed and sequenced in the routes. These strategies are applied only once before running the ILS+HC. On the other hand, the fourth strategy filters CS insertions that are infeasible for a particular fixed route.

Strategy 1 : This strategy estimates the minimum time τ needed to visit CS j ∈ F between two vertices i and h ∈ I ∪ {0}. This time is defined as the sum of a lower bound on the travel time (u) and a lower bound on the charging time (v) of any route serving customers i and h. Note that in any route serving customer i, t 0i is a lower bound on the route's duration from the start up to customer i. Similarly, note that in any route servicing customer h, t h0 is a lower bound on the time needed to complete a route from vertex h. Based on these two observations we can compute the minimum duration of a route visiting CS j between vertices i and h as u = t 0i + t ij + t jh + t h0 + p i + p h . To compute the minimum charging time v, we need to compute the minimum amount of energy (ec) that an EV coming from vertex i and traveling to vertex h must charge at CS j. This amount is the charge needed to recover the energy consumed to make the detour to j, that is, ec = e ij + e jh -e ih .

Because the battery level when the EV arrives at i in any eVRP-NL solution is unknown a priori, we consider that the battery is charged at j using the fastest charging rate (ρ 0j ). Then v = ec ρ 0j . It is clear that if τ = u + v > T max any route visiting CS j between customers i and h is infeasible. We therefore forbid this insertion in our MILP. Strategy 2 : This strategy computes a lower bound on the remaining energy that an EV must have at arrival to customer i to be able to visit CS j right after. Note that in terms of energy remaining, the best way to reach vertex i is to visit it right after fully charging at CS c(i) = arg min l∈F ∪{0} e li . If Q -e c(i)i < e ij any route visiting j after i is energy-infeasible and we can safely forbid this insertion in our MILP. Strategy 3 : Note that o = Q -e ji is a lower bound on the energy remaining when an EV arrives at customer i right after charging at CS j. Note also that o must be enough to at least close the route (reach the depot) or reach the closest CS in terms of energy consumption. If o < e ic(i) , where c(i) = arg min l∈F ∪{0} e il , any route visiting j before i is energy-infeasible and we can safely forbid this insertion.

Strategy 4 : This strategy estimates a lower bound on the new duration of a given fixed route if CS j ∈ F is inserted between vertices π(i) and π(i + 1) ∈ Π, i = n r . This bound, t , is defined as the sum of the new travel time (u) and a lower bound on the charging time (v). It is easy to see that u = t + t π(i)j + t jπ(i+1) -t π(i)π(i+1) . Similarly to Strategy 1, to estimate v, we consider that the battery is charged at j using the fastest charging rate. Therefore v = ec ρ 0j , where ec = e π(i)j + e jπ(i+1) -e π(i)π(i+1) is the charge needed to recover the energy consumed in the detour to j. If t = (u + v) > T max , inserting j between vertices π(i) and π(i + 1) leads to an infeasible route, so we can safely forbid this insertion.

Applying our preprocessing strategies to the MILP formulation for the 10-customer route of the example above, the model reduces to: 71 continuous variables, 40 integer variables, and 165 constraints. We solved the model in Gurobi in 0.02s.

Approach 2: greedy heuristic

Existing metaheuristics for eVRPs use different approaches to make charging decisions. One popular approach is the recharge relocation operator proposed by [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF] for the green vehicle routing problem with multiple technologies and partial recharges (GVRP-MTPR). This approach considers the insertion of only one CS per route. Starting from an energy-feasible fixed-route the procedure first deletes the current CS. Then, it tries to improve the charging decisions by inserting each CS into each arc of the fixed-route. To decide how much energy to charge at the inserted CS, the algorithm applies a simple rule: charge the minimum amount of energy needed to reach the depot (i.e., to complete the route). We propose here a heuristic to solve the FRVCP based on this approach.

Our heuristic works in two phases: location of CSs and charge setting. In the first phase, the heuristic iteratively inserts CSs into the arcs of the fixed-route Π in order to ensure the feasibility in terms of energy. In the second phase, the heuristic improves the charging decisions by adjusting the energy charged at each visited CS. Algorithm 3 describes the structure of our heuristic. It uses four important procedures trackBattery(•), sumNegative(•), totalTime(•), and copyAndInsert(•). Procedure trackBattery(•) computes the battery level Y i at each vertex i ∈ Π, assuming that the EV fully charges its battery at each visited CS. Note that Y i may take negative values. Procedure sumNegative(•) computes the sum of the battery levels with negative values (i.e., s = i∈Π min{0, Y i }). Procedure totalTime(•) computes the total time t of the route (supposing a full charging policy). Finally, procedure copyAndInsert(•) takes as input a fixed route, a CS, and a position in the route; and returns a copy of the fixed route with the CS inserted at the given position.

The heuristic starts by the location phase (line 2-27). After computing s for the current fixed-route Π, the heuristic enters the outer loop (line 7-27). In each pass through the inner loop (line 8-25), the heuristic: i) evaluates the insertion of a CS into each arc of Π assuming that the EV fully charges its battery, and ii) selects the insertion that maximizes s (lines 14-18). If s = 0 (i.e., the route is energy-feasible), the heuristic selects the insertion that minimizes t (line 19-23). Then, the heuristic performs the selected insertion (line 26). If the route Π is still energy-infeasible (i.e., s < 0), the heuristic starts again at line 8 and tries to insert additional CSs until feasibility in terms of energy is reached.

In the charge setting phase (line 28), the heuristic invokes procedure ruleMinEnergy(Π) to set the energy charged at each CS following the Felipe et al. (2014) rule. Finally, the heuristic evaluates if the route satisfies the maximum-duration constraint [START_REF] Villegas | Grasp/vnd and multi-start evolutionary local search for the single truck and trailer routing problem with satellite depots[END_REF](30)(31)(32)(33). The heuristic returns a boolean variable indicating whether or not the fixed-route is feasible (f ), and the route Π with the newly inserted CSs.

Computational experiments

In this section, we present three computational studies. The first study compares the quality of the solutions obtained by two versions of our metaheuristic. The second study evaluates the CPU time of our if s = 0 and t < t then return f , Π 35: end function metaheuristic, and assesses the impact of the preprocessing strategies. Finally, the third study analyses the charging decisions of the best solutions found.

Test instances for the eVRP-NL

To test our approach, we generated a new 120-instance testbed built using real data of EV configuration and battery charging functions. In order to ensure feasibility, we opted to generate our instances instead of adapting a existing dataset from the literature. To build the instances we first generated 30 sets of customer locations with {10, 20, 40, 80, 160, 320} customers. For each instance size, we generated 5 sets of customers location. We located the customers in a geographic space of 120 x 120 km using either a random uniform distribution, a random clustered distribution, or a mixture of both. For each of the 30 sets of locations we chose the customer location strategy using a uniform probability distribution. Our main motivation to choose a 120 x 120 geographic area was to build instances representing a semi-urban operation. These operations are the best suited applications for eVRPs. Indeed, in city operations routes tend to be sufficiently short to be covered without mid-route charging. On the opposite side, in rural operations routes tend to be long enough to required multiple mid-route charges but access to charging infrastructure is very limited (at least for 2016 standards).

From each of the 30 sets of locations we built 4 instances varying the level of charging infrastructure availability and the strategy used to locate the CSs. We considered two levels of charging infrastructure availability: low and high. To favor feasibility, for each combination of number of customers and infrastructure availability level we handpicked the number CSs as a proportion of the number of customers. We located the CSs either randomly or using a simple p-median heuristic. Our p-median heuristic starts from a set of randomly generated CS locations and iteratively moves those locations trying to minimize the total distance between the CSs and the customers. In our instances, we included three types of CSs: slow, moderate, and fast. For each CS we randomly selected the type using a uniform probability distribution.

The EVs in our instances are Peugeot Ion. This EV has a consumption rate of 0.125 kWh/km, and a battery of 16 kWh. Note that an EV with this characteristics is well suited to service applications such as homecare routing. In reality the exact energy consumption on an arc (e ij ) varies with parameters such as the cumulative elevation gain, the external temperature, the speed, and the use of peripherals (e.g., the radio). Nonetheless, for the sake of simplicity we followed the classical approach in the literature and assumed that the energy consumption on an arc is simply the EV's consumption rate multiplied by the arc's distance. To generate the charging functions we fit piece-wise linear functions to the real charging data for a 16 kWh battery provided by [START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF]. Figure 4.3 depicts our piece-wise linear approximations. Finally, the maximum route duration for every instance was fixed to 10 hours. Our 120 instances are publicly available at www.vrp-rep.org (Mendoza et al. 2014) 3 .

Parameter settings & experimental environment

As discussed in Section 4.3.2 our ILS+HC can be configured to solve the FRVCP using two different approaches, namely, the commercial solver and the greedy heuristic. We tested our algorithm with two different configurations. The first configuration (ILS(S)+HC) uses the commercial optimizer to solve the FRVCP in the GCI neighborhood, while the second configuration (ILS(H)+HC) uses the greedy heuristic. In both configurations, the algorithm solves the FRVCP in the split procedure using the greedy heuristic. This choice was guided by computational performance.

To fine tune the number of iterations of the ILS component (i.e., K), we conducted a short computational study. Our results showed that K = 80 provides the best trade-off between solution quality and computational performance.

We implemented our ILS in Java (jre V.1.8.0) and used Gurobi Optimizer (version 5.6.0) to solve the FRVCP and the set partitioning problem in the HC component. We set a time limit of 800 seconds on

Solution accuracy: optimal vs. heuristic charging decisions

Since the eVRP-NL is a new problem, there are no results or algorithms to benchmark against. To get an idea of the quality of the solutions delivered by our ILS+HC, we ran its two versions, ILS(S)+HC and ILS(H)+HC, on the twenty 10-customer instances and compared their results with optimal solutions found using the MILP presented in Chapter 3. Table 4.1 summarizes the results comparing the two algorithms under the light of four metrics: the number of optimal solutions found, the average and maximum gap with respect to the optimal solution over the 10 runs, and the average best gap4 . Appendix A.2 presents detailed results for the 20 instances.

The results suggest that our algorithms are able to deliver high-quality solutions for the eVRP-NL. As shown in the table, ILS(S)+HC matched the 20 optimal solutions while ILS(H)+HC matched 15. The spread in the number of optimal solutions and the values reported for Avg. Best Gap tip the balance towards ILS(S)+HC over ILS(H)+HC in terms accuracy. This result is not surprising since one could expect that making optimal, instead of heuristic, charging decisions would translate into higher quality solutions. A close look at the Avg. and Max. gaps reveals a less foreseeable result: ILS(S)+HC exhibits a significantly To compare the performance of ILS(S)+HC and ILS(H)+HC on more industrial-sized instances we ran both algorithms on the remaining 100 benchmarks. Table 4.2 summarizes the results. In the comparison we employed the same metrics introduced above, replacing the number of optimal solutions by the number of best known solutions (BKSs) found. As a reference, we included in the table the results delivered by the commercial solver running the MILP for 10h. It is worth mentioning that Gurobi reported integer solutions for only 25 out of the 100 instances. None of these 25 solutions has a certificate of optimality. Detailed results for each instance can be found in Appendix A. 2 The results in Table 4.2 confirm the conclusions drawn on our first experiment: optimally solving the FRVCP in the ICG neighborhood leads to a more accurate and stable solution method. According to our data, ILS(S)+HC not only found the 100 BKSs but also reported significantly lower Avg. and Max. gaps (1.51% and 4.44% vs. 7.51% and 28.68%).

The cost of making optimal charging decisions

In Section 4.4.3 we showed the benefits of making optimal charging decisions when solving the eVRP-NL. This benefits, however, comes at the price of losing computational performance. To estimate this loss, we measured the execution times of three versions of our algorithm on the 120 instances: ILS(S)+HC, ILS(H)+HC, and ILS(S')+HC. The latter is a version of ILS(S)+HC running without our preprocessing strategies (see Section4.3.2). Table 4.3 reports for each instance size, the average execution time (in seconds) of each algorithm. In addition, the last column of the table reports the average speedup between ILS(S )+HC and ILS(S)+HC (measured as the ratio of their execution times).

Not surprisingly, the results show that ILS(H)+HC is the fastest approach. It is interesting to note that in relative terms, the gap in computational performance between ILS(H)+HC and ILS(S)+HC decreases with the size of the instance. This behavior can be explained by the positive correlation between the instance size and the time needed to solve the set partitioning model in the HC phase. Indeed, the latter heavily depends on the size of the pool Ω which, in turn, is more correlated to the size of the instance than to the strategy used to generate the routes. In other words, on large instances both algorithms spend around the same (considerable) amount of time solving the HC phase. As a consequence, the performance gap in the ILS phase becomes less remarkable.

A second interesting conclusion from this study is that the preprocessing strategies have a remarkable positive impact on the computational performance of the algorithm. As the table shows, on average, the strategies are responsible for a 1.65 speedup on the execution time. The data also suggests that the speedup is independent of the size of the instance. A plausible explanation for this behavior is that for a given instance both ILS(S)+HC and ILS(S')+HC evaluate around the same number of routes during the ILS phase.

Characteristics of good eVRP-NL solutions

We analyze in this section the characteristics of the BKSs found in our experiments. We aim to provide the reader with some insight that may be useful when designing new solution methods for the eVRP-NL. In total, our BKSs are made up of 1,426 routes. Our first analysis concerns the fraction of those routes that exploit mid-route charging. Figure 4.4 presents the percentage of routes with and without visits to CSs grouped by instance size. The data shows that on average 71.47% of the routes in the BKSs visit at least one CS. This percentage is roughly the same for each instance size. This figure provides two insights. First, mid-route charging is a key element of good eVRP-NL solutions (probably because it gives algorithms flexibility to better pack the customers into the routes). Second, since charging decisions concern most of the routes making up a good solution, they play a critical role in its quality. The latter observation helps explaining the spread on accuracy between ILS(S)+HC and ILS(H)+HC found in the computational studies reported in Section 4.4.3.

The second analysis concerns the number of mid-route charges per route. Figure 4.5a and 4.5b present histograms of the number mid-route charges per route and the maximum number of mid-route charges per route on a solution. Figure 4.5a shows that among the routes performing mid-route charging, 58.58% do it once, 40.00% twice, and 1.43% three times. Although a large portion of the routes perform a single midroute charge, 85.83% of the solutions contain at least one route performing more than one mid-route charge (Figure 4.5b). This figures suggest that models and methods for the eVRP-NL can benefit from relaxing the at-the-most-one-visit-to-a-CS-per-route constraint that is sometimes used in the eVRP literature.

The third analysis concerns the energy recovered through mid-route charges. histogram of the average battery level (in % of the total battery capacity) after a mid-route charge. The numbers show that over 90% of the mid-route charges are partial charges (i.e., they do not fully charge the battery). This figure tips the importance of embedding components capable of making partial charging decisions into eVRP-NL solution methods 5 . A second interesting observation from the data displayed in Figure 4.6 comes from the percentage (around 12%) of mid-route charges that restore the battery above 80% of its capacity. As mentioned in the study reported in Chapter 3, one common assumption in the eVRP literature is that the battery can only be charge on the linear segment of the charging curve (which ends at roughly 80% of the battery capacity). Our data suggest that good eVRP-NL solutions often include routes with mid-route charges that take the battery level up to the non-linear part of the charging function.

Conclusion

This chapter introduces an extension of the electric vehicle routing problem which considers morerealistic assumptions about the battery charging process: the electric vehicle routing problem with nonlinear charging function (eVRP-NL). To solve the problem we propose an iterated local search (ILS) enhanced with heuristic concentration (HC). At the heart of the proposed method lays a neighborhood scheme consisting in solving a new variant of the fixed-route vehicle-charging problem (FRVCP). This problem consists in optimizing the charging decisions (where and how much to charge) of a route serving a fixed sequence of customers. Depending on the configuration, our ILS+HC solves the underlying FRVCP either using a greedy heuristic or a commercial solver. To improve the performance of the solver, we proposed four preprocessing strategies that eliminate infeasible detours to CSs. To assess the performance of our method, we built a set of 120 instances based on real EV and battery charging data. Tested on 20 small 10-customer instances, our method matched the optimal solutions for every instance. Experiments conducted on larger instances proved the value of equipping eVRP-NL algorithms with components capable of making optimal charging decisions. Finally, we analyzed the solutions delivered by our method aiming to provide fellow researchers with some insight into the characteristics of good eVRP-NL solutions. Our analysis concluded that good eVRP-NL solutions tend to use multiple mid-route charges, exploit partial recharges, and employ the non-linear segment of the battery charging function.

Introduction

Nowadays, service and utility companies have started to use electric vehicles (EVs) on their operations. Considering that EVs are a new technology, and that the companies currently have a fleet of conventional vehicles (CVs), the incorporation of EVs into their fleets is made gradually. For instance, the french electricity distribution company ENEDIS has started to incorporate 2,000 EVs in its fleet of 20,000 vehicles (ERDF 2014). Similarly, the spanish electricity company Endesa has started to introduce EVs in its fleet; they expected to have 3,000 EVs by 2020 (H&E 2016). The incorporation of EVs into existing fleets of CVs gives birth to what is commonly known in the literature as vehicle routing problems with electric and conventional vehicles (Goeke & Schneider 2015).

In this chapter, we introduce what we call the technician routing problem with conventional and electric vehicles (TRP-CEV). The TRP-CEV is largely inspired by a problem faced by ENEDIS in their daily operations. The problem consists in routing a set of technicians to serve a set of geographically scattered requests within their time windows. The problem also considers the working schedules, the lunch breaks, and skills of the technicians. The technicians use a fixed fleet composed of CVs and EVs of different types. Each type of EV has its own autonomy charging protocol. The latter determines in which charging stations (CSs) the EV can charge its battery. We assume that the charging functions are nonlinear as proposed in Chapter 3. In the TRP-CEV decisions include the vehicle-to-technician assignment, the sequencing of the routes followed by the technicians, and the battery charging program for the EVs (where and how much to charge). The objective function seeks to minimize the total cost, defined as the sum of the travel costs, battery charging fixed costs, the parking cost at CSs, and the fixed cost of each technician. Following the definition of [START_REF] Lahyani | Rich vehicle routing problems: From a taxonomy to a definition[END_REF] the TRP-CEV lays in the category of rich vehicle routing problem (RVRP) 1 .

It is widely accepted that large-scale instances of complex RVRPs are out of the reach of exact methods. Moreover, in industrial contexts in which the problem must be solved on a daily (or even weekly) basis, companies tend to prefer the more-balanced trade off between solution quality and speed offered by metaheuristics. Since one of our primary objectives is to provide ENEDIS with an optimization algorithm that can eventually be deployed to production, we leaned towards designing a parallel matheuristic (PMa) for their problem. The proposed approach works in two phases. In the first phase, the algorithm decomposes the TRP-CEV into a set of "easier-to-solve" vehicle routing problems with time windows and lunch breaks (VRPs-TWLB) and solves these problems in parallel using a greedy randomized adaptive search procedure (GRASP). In the second phase, the algorithm assembles a solution to the TRP-CEV by solving a set-covering formulation over a set of high-quality routes found for the independent VRP-TWLB.

The main contribution of this chapter are fourfold. First, we introduce the TRP-CEV. Second, we propose a parallel matheuristic to solve the TRP-CEV which can be extended to other problems with competitive results. Third, we develop a set of real instances based on the ENEDIS operation. Fourth, we study the impact of different fleet compositions on relevant variables for ENEDIS.

The remainder of this chapter is organized as follows. Section 5.2 reviews the related literature. Section 5.3 describes the problem. Section 5.4 introduces our parallel matheuristic. Section 5.5 introduces the fixedroute vehicle-charging problem with time windows, which arises as a subproblem of TRP-CEV. Section 5.6 presents the results of a computational experiment with ENEDIS real instances. Finally, Section 5.7 concludes the chapter. For easy reference, Section 5.8 summarizes the notation of this chapter.

Literature review

The TRP-CEV combines features from different problems in the literature: the technician routing problem (TRP), the vehicle routing problem with lunch breaks (VRP-LB) and the heterogeneous vehicle routing problem (HVRP). In the remainder of this chapter, we present a brief review of literature on these problems.

One of the first problems in the field of TRPs is the field technician scheduling problem (FTSP) by [START_REF] Xu | Effective heuristic procedures for a field technician scheduling problem[END_REF]. The FTSP consists in routing a set of technicians to serve a set of geographically scattered requests. Each request has a time window, and a priority (represented by a weight). Each technician has a working schedule, and a skill proficiency level for each request. The objective functions for the FTSP is to maximize the number of served requests weighted their priority and minimize the total working time. [START_REF] Dohn | The manpower allocation problem with time windows and job-teaming constraints: A branch-and-price approach[END_REF] proposed an interesting variant of the TRP, where some requests require collaboration between technicians, and therefore the technicians must initiate the execution simultaneously. [START_REF] Pillac | A parallel matheuristic for the technician routing and scheduling problem[END_REF] studied a version of the TRP where the technician-task compatibility includes spare parts and tools. The technicians have the opportunity to replenish their tools and spare parts at a central depot at any time. [START_REF] Kovacs | Adaptive large neighborhood search for service technician routing and scheduling problems[END_REF] consider that technicians have different proficiency levels in several skills, and the option that technicians can be grouped into teams to serve specifics request. Recently, [START_REF] Zamorano | Branch-and-price approaches for the multiperiod technician routing and scheduling problem[END_REF] extend the TRP considering multiple periods. In their problem, technicians are assigned every day to a team, and the composition of the team remains fixed for the duration day. For a detailed review of the TRP and other related problems the reader is referred to (Castillo-Salazar et al. 2016). As can be seen the TRPs have many variations. In the case of the TRP-CEV, it only consider the technician-task compatibility that depends of the required skills of the requests and the skills of the technicians.

VRPs with lunch breaks arise mainly in service industries. For instance, in home healthcare operations where nurses and doctors visit a set of customers at home [START_REF] Liu | Mathematical model and exact algorithm for the home care worker scheduling and routing problem with lunch break requirements[END_REF][START_REF] Cheng | A home health care routing and scheduling problem[END_REF]). In the reparation or maintenance service, where the technicians serve a set of requests [START_REF] Kovacs | Adaptive large neighborhood search for service technician routing and scheduling problems[END_REF][START_REF] Bostel | Multiperiod planning and routing on a rolling horizon for field force optimization logistics[END_REF]). To the best of our knowledge, all studies tacking VRP-LB consider a time window to start the lunch break and a specified duration. Futhermore, most of the studies assume that the lunch break is taken at a customer location [START_REF] Liu | Mathematical model and exact algorithm for the home care worker scheduling and routing problem with lunch break requirements[END_REF][START_REF] Cheng | A home health care routing and scheduling problem[END_REF][START_REF] Kovacs | Adaptive large neighborhood search for service technician routing and scheduling problems[END_REF]). The only exception is the study by [START_REF] Bostel | Multiperiod planning and routing on a rolling horizon for field force optimization logistics[END_REF], where the authors assume that the lunch break is taken at specified locations (restaurants). In contrast with the previous studies, for the TRP-CEV the lunch break is defined between two time points in the working time, where the technician must stops any activity to take the lunch break.

In the last 30 years, the HVRPs have been one the most studied problem in the VRP field (Koç et Despite in the literature there are several studies about the problems that compose the TRP-CEV, as far as we know no other authors have studied a problem that combines all those components. Moreover, the TRP-CEV is the only problem that introduces EVs in the technician operation context, taking into account the existing real applications of EVs in operation with technician (ERDF 2014, H&E 2016).

Problem description

Formally, the TRP-CEV is defined on a directed and complete graph G = (N, A). The vertex set N = {0} ∪ C ∪ F is made up of a depot (vertex 0), a set of requests C, and a set of CSs F . Each request i ∈ C has a service time p i and a time window [ec i , lc i ], and requires a technician with specific skills from the set of skills S. Each CS i ∈ F has a parking cost pc i . The request are served by a set of T technicians. Each technician u ∈ T has the following characteristics: a fixed cost f u ; a subset of skills s u ⊆ S; a working schedule between st u and et u , that includes a lunch break between sl u and el u ; and a consumption factor f c u associated to his or her driving profile (e.g., sport, regular and eco). The technicians drive a fixed fleet of heterogeneous vehicles composed by different types of CVs and EVs. The set of vehicles types K is made up of the subset of types of CVs (K c ) and the subset of types of EVs (K e ). There are m k vehicles of type k ∈ K available in the fleet. Each type of vehicle k ∈ K has an unitary travel cost tc k (expressed in e/km). Additionally, EVs of type k ∈ K e have a fixed cost g k for recharging the battery (expressed in e and representing the degradation cost), and a battery capacity Q k (expressed in kWh). Due to technological incompatibilities, each type of EV k ∈ K e can charge only in a subset F k of the CS. The charging function of EVs is modelled using a discrete function. Parameter a ikl represents the charging time of a vehicle of type k ∈ K e charging from zero to l percent of the battery at CS i ∈ F (with l ∈ B = 1, 2, ..., 100). Finally, the set A = {(i, j) : i, j ∈ N, i = j} corresponds to arcs connecting vertices in N . Each arc (i, j) has three associated nonnegative values: a travel time t ij , a distance d ij , and an energy consumption e ijk for each type of EV k ∈ K e .

In the TRP-CEV the objective is to find a set of routes of minimum total cost, defined as the sum of the total travel cost, fixed charging cost, parking cost, and the fixed cost of each technician, satisfying the following conditions: each request is served exactly once within its time window by a technician with the required skills; the level of the battery when the EVs arrive at any vertex is nonnegative; the EVs only charge at CSs with compatible technology; each technician works only during his or her working schedule; each technician takes a lunch break; and each route starts and ends at the depot.

Mixed-integer linear programming formulation

To allow several visits to the CSs, we introduce the set F that contains the set F and copies of each CS, and the set N = {0} ∪ C ∪ F . Furthermore, since we need to track the information of each vehicle, we introduce the set of vehicles V = V c ∪ V e that contains all the CVs (V c ) and EVs (V e ) available (i.e.,

|V c | = k∈Kc m k and |V e | = k∈Ke m k ).
Furthermore, we introduce the following binary parameters: n il is equal to 1 if request i ∈ C needs a technician with skill l ∈ S, and 0 otherwise. Likewise, m ul is equal to 1 if technician u ∈ T has skill l ∈ S, and 0 otherwise. Finally, h ik is equal to 1 if the EV k ∈ V e can charge at CS i ∈ F . In addition, the real parameter S max is the maximum charging time.

To model the TRP-CEV using a mixed-integer linear programming (MILP) formulation, we use the following decision variables: x ijku is equal to 1 if technician u ∈ T travels from vertex i to j ∈ N in vehicle k ∈ V , and 0 otherwise. Variable τ iu tracks the arrival time of the technician u ∈ T at vertex i ∈ N . Variable y ik tracks the battery level of EV k ∈ V e upon departure from vertex i ∈ N . Variables q ik and o ik represent the battery level when EV k ∈ V e arrives at and departs from CS i ∈ F . Variables α ikl and β ikl are equal to 1 if the battery level is l ∈ B when vehicle k ∈ V e arrives at and departs from the CS i ∈ F respectively. Variable ∆ ik represents the time spent by vehicle k ∈ V e at CS i ∈ F . Variable w iju is equal to 1 if technician u ∈ T takes the lunch break on his way from i to j ∈ N , and 0 otherwise. Variable z iju is equal to 1 if the technician u ∈ T takes the lunch break after travelling from vertex i to j ∈ N and 0 otherwise. Variable z iju is equal to 1 if the technician u ∈ T takes the lunch break before travelling from vertex i to j ∈ N , and 0 otherwise. The MILP follows:

Min i∈N j∈N k∈V u∈T x ijku • d ij • tc k + i∈F k∈Ve ∆ ik • pc i + i∈N j∈F k∈Ve u∈T x ijku • g k + j∈N k∈V u∈T x 0jku • f u (5.1)
Subject to:

i∈N k∈V u∈T x ijku = 1, ∀j ∈ C (5.2) i∈N x ijku ≤ 1, ∀j ∈ F , ∀k ∈ V e , ∀u ∈ T (5.3) k∈V u∈T x ijku ≤ 1, ∀i, j ∈ N (5.4) j∈N u∈T x ijku ≤ 1, ∀i ∈ N , ∀k ∈ V (5.5) j∈N k∈V x ijku ≤ 1, ∀i ∈ N , ∀u ∈ T (5.6) j∈N x 0jku ≤ 1, ∀k ∈ V, ∀u ∈ T (5.7) j∈N x jiku - j∈N x ijku = 0, ∀i ∈ N , ∀k ∈ V, ∀u ∈ T (5.8) i∈N k∈V u∈T x ijku • m ul ≥ n jl , ∀j ∈ C, ∀l ∈ S (5.9) i∈N x ijku ≤ h jk , ∀j ∈ F , ∀k ∈ V e , ∀u ∈ T (5.10) e ijk • f c u • x ijku -(1 -x ijku ) • Q k ≤ y ik -y jk ≤ e ijk • f c u • x ijku + (1 -x ijku ) • Q k , ∀i ∈ N , ∀j ∈ C, ∀k ∈ V e , ∀u ∈ T (5.11) e ijk • f c u • x ijku -(1 -x ijku ) • Q k ≤ y ik -q jk ≤ e ijk • f c u • x ijku + (1 -x ijku ) • Q k , ∀i ∈ N , ∀j ∈ F , ∀k ∈ V e , ∀u ∈ T (5.12) y ik ≥ e i0k • f c u • x i0ku , ∀i ∈ N , ∀k ∈ V e , ∀u ∈ T (5.13) y 0k = Q k , ∀k ∈ V e (5.14) y ik = o ik , ∀i ∈ F , ∀k ∈ V e (5.15) q ik ≤ o ik ,
∀i ∈ F , ∀k ∈ V e (5.16)

q jk ≤ Q k • i∈N u∈T x ijku , ∀j ∈ F , ∀k ∈ V e (5.17)
o jk ≤ Q k • i∈N u∈T x ijku , ∀j ∈ F , ∀k ∈ V e (5.18)
q ik Q k • 100 ≥ l∈B l • α ikl , ∀i ∈ F , ∀k ∈ V e (5.19)
q ik Q k • 100 -1 ≤ l∈B l • α ikl , ∀i ∈ F , ∀k ∈ V e (5.20) l∈B α jkl ≤ i∈N u∈T x ijku , ∀j ∈ F , ∀k ∈ V e (5.21) o ik Q k • 100 ≤ l∈B l • β ikl , ∀i ∈ F , ∀k ∈ V e (5.22) o ik Q k • 100 + 1 ≥ l∈B l • β ikl , ∀i ∈ F , ∀k ∈ V e (5.23) l∈B β jkl ≤ i∈N u∈T
x ijku , ∀j ∈ F , ∀k ∈ V e (5.24)

∆ ik = l∈B a ikl • β ikl - l∈B a ikl • α ikl , ∀i ∈ F , ∀k ∈ V e (5.25)
i∈N j∈N

w iju = i∈N k∈V
x 0jku , ∀u ∈ T (5.26)

w iju ≤ k∈V x ijku , ∀i ∈ N , ∀j ∈ N , ∀u ∈ T (5.27) τ iu + w iju • (el u -sl u ) + (t ij + p i ) • x ijku - (et u + (el u -sl u )) • (1 -x ijku ) ≤ τ ju , ∀i ∈ C, ∀j ∈ N , ∀k ∈ V e , ∀u ∈ T (5.28) τ iu + w iju • (el u -sl u ) + ∆ ik + t ij • x ijku - (S max + et u + (el u -sl u )) • (1 -x ijku ) ≤ τ ju , ∀i ∈ F , ∀j ∈ N , ∀k ∈ V e , ∀u ∈ T (5.29) τ iu + p i ≤ sl u • j∈N w iju + et u • (1 - j∈N w iju ), ∀i ∈ C, ∀k ∈ V, ∀u ∈ T (5.30) τ iu + ∆ ik ≤ sl u • j∈N w iju + et u • (1 - j∈N w iju ), ∀i ∈ F , ∀k ∈ V e , ∀u ∈ T (5.31) τ iu ≥ el u • j∈N w jiu + st u • (1 - j∈N w jiu ), ∀i ∈ N , ∀k ∈ V, ∀u ∈ T (5.32) sl u • z iju -τ iu -p i + et u • (1 -z iju ) ≥ t ij • z iju , ∀i, j ∈ N , ∀u ∈ T (5.33) τ ju -el u • z iju ≥ t ij • z iju , ∀i, j ∈ N , ∀u ∈ T (5.34) z iju + z iju = w iju , ∀i, j ∈ N , ∀u ∈ T (5.35) τ iu -t 0i ≥ st u , ∀i ∈ N , ∀u ∈ T, (5.36)
τ iu + p i + t i0 ≤ et u , ∀i ∈ C, ∀u ∈ T (5.37) τ iu + ∆ ik + t i0 ≤ et u , ∀i ∈ F , ∀k ∈ V e , ∀u ∈ T (5.38) ec i ≤ τ iu ≤ lc i , ∀i ∈ C, ∀u ∈ T (5.39) x ijku ∈ {0, 1} ∀i, j ∈ N , ∀k ∈ V e , ∀u ∈ T (5.40) x ijku ∈ {0, 1} ∀i ∈ N , j ∈ N \ F , ∀k ∈ V c , ∀u ∈ T (5.41) α ikl ∈ {0, 1}, β ikl ∈ {0, 1}, ∀i ∈ F , ∀k ∈ V e , ∀l ∈ B (5.42) τ iu ≥ 0, ∀i ∈ N , ∀u ∈ T (5.43) y ik ≥ 0, ∀i ∈ N , ∀k ∈ V e (5.44) q ik ≥ 0, o ik ≥ 0, ∆ ik ≥ 0, ∀i ∈ F , ∀k ∈ V e (5.45) w iju ∈ {0, 1}, z iju ∈ {0, 1}, z iju ∈ {0, 1}, ∀i, j ∈ N , ∀u ∈ T (5.46)
The objective function (5.1) seeks to minimize the total cost (travel cost, parking cost, fixed charging cost, fixed cost of the technicians). Constraints (5.2) ensure that each request is visited once. Constraints (5.3) ensure that each copy if a CSs is visited at most once. Constraints (5.4) ensure that each arc is used by at most one vehicle and technician. Constraints (5.5-5.7) ensure that at most only one technician is allocated to one vehicle (and vice versa). Constraints (5.8) impose the flow conservation. Constraints (5.9) ensure that each request is visited by a technician with the appropriate skills. Constraints (5.10) ensure that EV i ∈ V e only charges at CSs with compatible technology. Constraints (5.11) and (5.12) track the battery level at each vertex for the EVs. Constraints (5.13) ensure that, if the EV travels between a vertex and the depot, it has enough remaining energy to reach its destination. Constraints (5.14) impose that the EVs depart from the depot with fully charged battery. Constraints (5.15) reset the battery tracking to o i upon departure from CS i ∈ F . Constraints (5.16) couple the battery level when an EV arrives to and departs from any CS. Constraints (5.17) and (5.18) ensure that the charge level is bounded by the battery capacity when an EV arrives at and departs from CS j ∈ F . Constraints (5. 19-5.24) define the battery level in percentage when vehicle k ∈ V e arrives at or departs from the CS i ∈ F . Constraints (5.25) define the time spent charging at a CS. Constraints (5.26) ensure that the technician u ∈ T has to take the lunch break. Constraints (5.27) ensure that the lunch break is taken between two visited vertices. Constraints (5.28) and (5.29) track the arrive time at each vertex. Constraints (5.30)-(5.32) ensure that technician u ∈ T takes the lunch break between sl u and el u . Constraints (5.33-5.35) ensure that the technician takes the lunch after or before travelling to vertex j ∈ N . Constraints (5.36-5.38) impose the working schedule of the technician. Constraints (5.39) ensure that every request is served within its time window. Finally, constraints (5.40-5.46) define the domain of the decision variables. Note that variables x ijku are not defined for j ∈ F and k ∈ V c

Parallel matheuristic

Considering the complexity of th TRP-CEV in terms of the number of constraints, we decide to use an approach based on decomposition. As the name suggest, a decomposition approach decomposes the problem in a set of sub-problems that consider less information or less constraints than the original problem. Each sub-problem is solved individually, and the solutions found for each sub-problem are used to find a solution of the original problem. This kind of approaches has been successful in the solution of different rich VRPs [START_REF] Doerner | Survey: Matheuristics for Rich Vehicle Routing Problems[END_REF].

One alternative to address the set of sub-problems is to solve them in parallel using metaheuristics. Parallel metaheuristics are commonly used to solve VRPs [START_REF] Crainic | Parallel solution methods for vehicle routing problems[END_REF]). They focus either on decomposing the solution approach (e.g., using different local search strategies in parallel), or on decomposing the problem domain (e.g., into geographic zone, or type of vehicle). We follow the latter approach. Successful examples of this approach include [START_REF] Taillard | Parallel iterative search methods for vehicle routing problems[END_REF]) that proposes a parallel tabu search for the capacitated VRP. In his method, the problem is divided spatially into regions and a problem is solved individually for each region with tabu search. [START_REF] Doerner | A parallel version of the d-ant algorithm for the vehicle routing problem[END_REF] follow the same decomposition principle but use and savings based ants system [START_REF] Reimann | D-ants: Savings based ants divide and conquer the vehicle routing problem[END_REF]) for each sub-problem. [START_REF] Taillard | A heuristic column generation method for the heterogeneous fleet VRP[END_REF] propose a heuristic column generation for the HVRP. This approach decomposes the HVRP in a set of homogeneous VRP (one for each type of vehicle), and solve each VRP using an adaptive memory programming procedure. Then, set partitioning formulation assembles a HVRP solution.

Our method shares the spirit of [START_REF] Taillard | A heuristic column generation method for the heterogeneous fleet VRP[END_REF]. Algorithm 4 describes the general structure of our PMa. In order to avoid symmetry problems, the algorithm starts by calling the procedure groupTechnicians(T ) -line 2-. This procedure groups the technicians that have the same characteristics (i.e., skills, fixed cost, working schedule and lunch break), and generates the set E of technician types. Then, the algorithm invokes the procedure buildCouples(E,K)-line 3-. It builds the set P of all possible couples computing the cartesian product P = E × K = {(u, k)|u ∈ E, k ∈ K}. Then, the algorithm starts in a parallel phase -lines 5-9-. In this phase, for each couple p ∈ P it solves a VRP-TWLB p or an eVRP-TWLB p on each thread (the problem to be solved depends of the type of vehicle used a CV or an EV). To solve each (e)VRP-TWLB p we ignore the fleet size and the limit on the number of technicians. Before to solve each (e)VRP-TWLB p , we call the procedure feasibleRequests(G,p) that identifies the feasible requests C p for the couple p ∈ P . The details of this procedure are explained in Section 5.4.1. Then, the algorithm invokes the metaheuristic GRASP(p,C p ,G) to solve the corresponding (e)VRP-TWLB p -line 7-. This metaheristic returns a set Ω p of feasible routes for the underlying problem. The set Ω p joins the global set Ω. After completing the parallel phase, the algorithm calls the procedure SetCovering(G,Ω,K,E)line 10-, which solves an extended set covering formulation over Ω to find a feasible TRP-CEV solution. This formulation includes constraints on the number of technicians and vehicles. Finally, the resulting solution σ is reported by PMa -line 11-. E ←-groupTechnicians(T )

3: P ←-buildCouples(E,K) 4: Ω ←-∅ 5:
parallel for each p ∈ P 6:

C p ←-feasibleRequests(G,p) 7: Ω p ←-GRASP(p,C p ,G) 8: Ω ←-Ω ∪ Ω p 9: end for 10: s ←-SetCovering(G,Ω,K,E) 11:
return σ 12: end function

Identifying feasible requests

Considering that each couple p ∈ P has specific features associated to the type of technician, and the type of vehicle, before solving the (e)VRP-TWLB p we identify the feasible requests C p for the couple p ∈ P using a set of criteria. To explain those criteria, we use the notation introduced in Section 5.3, and we denote u * and k * be the type of technician and type of vehicle of couple p ∈ P .

Skill criterion: If the type of technician u * has at least the required skills of the request i ∈ C (i.e., m u * l ≥ n il , ∀l ∈ S), then the request i ∈ C is part of the set of feasible requests (i.e, C p = C p ∪{i}).

Time window -working schedule criterion: If the latest time lc i of time window request i ∈ C minus the travel time from depot to request i is greater than the start time of the working schedule of the type of technician u * (i.e, lc i -t 0i ≥ st u * ); and the earliest time lc i of time window request i ∈ C plus the its service time plus the travel time from request i to depot is lower than the end time of the working schedule of the type of technician u * (i.e,

ec i + p i + t i0 ≤ et u * ), then C p = C p ∪ {i}. Time window -lunch break criterion: If the time window [ec i , lc i ] of request i ∈ C is not completely contained in the lunch break of the type of technician u * (i.e., [ec i , lc i ] ⊂ [sl u , el u ]), then C p = C p ∪ {i}.
Autonomy criterion: If the battery capacity Q k * of the type of EV k * is greater than the minimum energy needed to reach request i ∈ C from any CSs or the depot (i.e., Q k * ≥ min j∈F ∪{0} {e jik * }), then C p = C p ∪ {i}.

GRASP

To solve the (e)VRP-TWLB p associated to the couple p ∈ P , we developed a GRASP [START_REF] Feo | Greedy randomized adaptive search procedures[END_REF]. This metaheuristic is a memory-less multi-start method in which local search is applied to several initial solutions constructed with a greedy randomized heuristic. GRASP can be hybridized in different ways, for instance by replacing the local search with an other metaheuristic such as tabu search, simulated annealing, variable neighborhood search (VND), among others [START_REF] Resende | Metaheuristic hybridization with greedy randomized adaptive search procedures[END_REF]. In this case, to build the initial solution we use a randomized Solomon heuristic (RSH). Then at each iteration of the GRASP, we improve the current solution using a VND [START_REF] Mladenović | Variable neighborhood search[END_REF]. At each iteration we add a set of routes to the pool of routes Ω p . To favour the diversity of the routes in this pool, we add to Ω p the routes of the initial solution found by RSH and the routes of the local optimum after the VND. The execution of GRASP stops after I iterations and we add the routes of Ω p to the global pool Ω (i.e., Ω = Ω p ∪ Ω). A brief description of the components of the GRASP follows:

Initial solution

We generate initial solutions for the (e)VRP-TWBL p using a randomized Solomon heuristic [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF]. In this heuristic, the seed request is randomly selected, and instead of selecting the best insertion, an insertion is randomly selected among the κ best insertions.

To evaluate the feasibility of a request insertion, we first check if it satisfies the time window, working schedule, and lunch break constraints. Using the concepts introduced by Savelsbergh (1992), this verification is done in constant time. Then, if request is to be inserted on a route covered by an EV, we evaluate if the insertion satisfies the energy constraint. If it does not, we try to turn the insertion feasible by inserting visits to CSs into the route. If we cannot achieve our goal, we discard that request insertion.

The problem of inserting visits to CSs into an energy-infeasible route (i.e., a route that is infeasible in terms of the energy autonomy of the vehicle) can be modelled using a fixed-route vehicle-charging problem (FRVCP) as in Chapter 4. In this case, the FRVCP has to consider time windows and the working schedule. We refer to this extension of the FRVCP simply as the FRVCP with time windows (FRVCP-TW). To favour the computational efficiency during the construction of initial solutions, we use a simple heuristic to solve the FRVCP-TW. This heuristic is explained in Section 5.5.

Variable neighborhood descent

To improve the solution generated by the RSH we use a VND based on three local search operators 2 . The first two operators, namely, relocated and swap, focus on the request sequencing decisions. On the other hand, the third operator is an ejection chain, aimed at reducing the number of routes [START_REF] Rousseau | Using constraint-based operators to solve the vehicle routing problem with time windows[END_REF]. To evaluate the feasibility of the moves performed by the operators, we use the same steps employed during the generation of the initial solutions.

Global charging improvement

If the GRASP is solving an eVRP-TWLB, after the VND we try to improve the charging decisions using a procedure based on the approach introduced in Section 4.3.2. This procedure called global charging improvement, or GCI, is applied to every route visiting at least one CS. The procedure works as follows. First, it removes from the route all visits to CSs. If the resulting route is energy-feasible, the procedure stops. If the route is energy-infeasible, the procedure solves a FRVCP-TW on the concerned route. In this case, to solve the FRVCP-TW we use a commercial optimizer to solve an MILP formulation of the problem. Full details on the MILP formulation of the FRVCP-TW are given in Section 5.5.

Set covering

To find a feasible TRP-CEV solution, we solve a set-covering (SC) formulation over the pool of routes Ω. To formulate the SC we use the following parameters: c r is the cost of the route r ∈ Ω, b ri is equal to 1 if route r ∈ Ω serves request i ∈ C and 0 otherwise, pv rk is equal to 1 if route r ∈ Ω is performed by a vehicle of type k ∈ K and 0 otherwise, h ru is equal to 1 if route r ∈ Ω is performed by a technician of type u ∈ E and 0 otherwise, and n u is the number of technicians of type u ∈ E. We denote by χ r a binary decision variable that takes the value of 1 if route r ∈ Ω is selected and 0 otherwise. We formulate the SC as follows:

min r∈Ω c r • χ r (5.47) Subject to r∈Ω b rj • χ r ≥ 1 ∀j ∈ C (5.48) r∈Ω pv rk • χ r ≤ m k ∀k ∈ K (5.49) r∈Ω h ru • χ r ≤ n u ∀u ∈ E (5.50) χ r ∈ {0, 1} ∀r ∈ Ω (5.51)
The objective function (5.47) seeks to minimize the total cost. Constraints (5.48) ensure that each request is served at least once. Constraints (5.49) impose the maximum number of vehicles of each type. Constraints (5.50) impose the maximum number of technicians of each type. Finally, constraints (5.51) define the domain of the decision variables.

Considering that the requests must be served exactly once, a set-partitioning problem may seem more appropriate. However, as pointed out by [START_REF] Pillac | A parallel matheuristic for the technician routing and scheduling problem[END_REF] for a similar technician problem, finding a good combination of routes that visit all requests exactly once could be complex in some cases. Therefore, we opted to use a SC formulation instead of a set-partitioning formulation. In the case that the solution found by the SC formulation serves a request more than once, we remove the most costly duplicate visits using a greedy heuristic.

The fixed-route vehicle charging problem with time windows

As previously mentioned, our GRASP for the eVRP-TWLB relies on solving the FRVCP-TW when we have to take the charging decision of a given route. This problem is an extension of the FRVCP introduced in Section 4.3, that includes the time window constraints, and the working schedule of the technician.

Let Π = {π(0), π(1), ..., π(i), ..., π(j), ..., π(n r )} be a route that violates the energy constraint, where π(0) and π(n r ) represent the depot. The feasibility of Π may be restored by inserting visits to CSs. The route Π is performed by a technician u ∈ T using an EV of type k ∈ K e . The technician has a working schedule [st u , et u ] and a lunch break [sl u , el u ]. For the route Π, we consider as fixed the decisions between which vertices the technician takes the lunch break. Let π(s) and π(e) be the requests served right before and after the lunch break. To model the lunch break into the route Π, we create a dummy vertex π(l), which has a time window [sl u , sl u ], and a service time p l = el u -sl u . Then, we define Π = {π(0), ..., π(s), π(l), π(e), ..., π(n r )} as the route that includes the dummy vertex π(l). Using the notation introduced in section 5.3, if the technician takes the lunch break before travelling from π(s) to π(e), then d π(s)π(l) = 0, t π(s)π(l) = 0, e π(s)π(l)k = 0, and d π(l)π(e) = d π(s)π(e) , t π(l)π(e) = t π(s)π(e) , e π(s)π(l)k = e π(s)π(e)k . On the contrary, if the technician takes the lunch break after travelling from π(s) to π(e), then d π(s)π(l) = d π(s)π(e) , t π(s)π(l) = t π(s)π(e) , e π(s)π(l)k = e π(s)π(e)k , and d π(l)π(e) = 0, t π(l)π(e) = 0, e π(s)π(l)k = 0.

In the FRVCP-TW the objective is to find the charging decisions that minimize the sum of the travel cost, the fixed charging cost, and the parking cost, satisfying the following conditions: the level of the battery when the vehicle arrives at any vertex is nonnegative; the requests must be visited within their time windows; and the route satisfies the working schedule and the lunch break of the technician.

We model the FRVCP-TW using an MILP formulation with the following decision variables: variable ε π(i)j is equal to 1 if the vehicle charges at CS j ∈ F before visiting vertex π(i) ∈ Π. Variable φ π(i) tracks the battery level. If ε π(i)j = 0, φ π(i) is the battery level when the EV arrives at vertex π(i). On the other hand, if π(i)j = 1, φ π(i) is the battery level when the vehicle arrives at CS j ∈ F right before visiting vertex π(i). Variable γ π(i)j represents the energy charged at CS j ∈ F before visiting vertex π(i) ∈ Π. Variables q π(i)j and o π(i)j represent the battery level when the EV arrives at and departs from CS j ∈ F before visiting vertex π(i) ∈ Π. Variables α π(i)jl and β π(i)jl are equal to 1 if the battery level is l ∈ B when the EV arrives at and departs from the CS j ∈ F before visiting vertex π(i) ∈ Π, respectively. Variable δ π(i)j represents the time spent at CS j ∈ F before visiting vertex π(i) ∈ Π. Variable τ π(i) represent the arrival time at the vertex π(i) ∈ Π. The MILP formulation of the FRVCP-TW follows:

min π(i)∈Π\{π(0)} j∈F tc k • (d π(i-1)j + d jπ(i) -d π(i-1)π(i) ) • ε π(i)j + g k • ε π(i)j + pc j • δ π(i)j
(5.52) Subject to

φ π(1) = Q k - j∈F e π(0)jk • ε π(1)j -e π(0)π(1)k • (1 - j∈F ε π(1)j ) (5.53) φ π(i) = φ π(i-1) + j∈F γ π(i-1)j - j∈F e jπ(i-1)k • ε π(i-1)j - j∈F e π(i-1)jk • ε π(i)j -e π(i-1)π(i)k • (1 - j∈F ε π(i)j ) ∀π(i) ∈ Π \ {π(0), π (1)
φ π(i) -Q k • (1 -ε π(i)j ) ≤ q π(i)j ≤ φ π(i) + Q k • (1 -ε π(i)j ) ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.57) o π(i)j = q π(i)j + γ π(i)j ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.58) l∈B l • α π(i)jl ≤ q π(i)j Q k • 100 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.59) l∈B l • α π(i)jl ≥ q π(i)j Q k • 100 -1 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.60) l∈B l • β π(i)jl ≥ o π(i)j Q k • 100 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.61) l∈B l • β π(i)jl ≤ o π(i)j Q k • 100 + 1 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.62) l∈B α π(i)jl ≤ ε π(i)j ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.63) l∈B β π(i)jl ≤ ε π(i)j ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.64) δ π(i)j = l∈B a jkl β π(i)jl - l∈B a jkl α π(i)jl ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.65) τ π(i) ≥ τ π(i-1) + p π(i-1) + j∈F ε π(i)j • t π(i-1)j + t jπ(i) + t π(i-1)π(i) • 1 - j∈F ε π(i)j + j∈F δ π(i)j ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (5.66) ec π(i) ≤ τ π(i) ≤ lc π(i) ∀π(i) ∈ Π \ {π(0), π(n r )} (5.67) st u ≤ τ π(i) ≤ et u ∀π(i) ∈ Π (5.68) j∈F ε π(i)j ≤ 1 ∀π(i) ∈ Π (5.69) γ π(i)j ≤ Q k • ε π(i)j ∀π(i) ∈ Π, ∀j ∈ F (5.70) γ π(i)j ≥ 0, δ π(i)j ≥ 0 ∀π(i) ∈ Π, ∀j ∈ F (5.71) ε π(i)j ∈ {0, 1} ∀π(i) ∈ Π, ∀j ∈ F (5.72) φ π(i) ≥ 0 ∀π(i) ∈ Π (5.73) α π(i)jl ∈ {0, 1}, β π(i)jl ∈ {0, 1}, ∀π(i) ∈ Π, ∀j ∈ F, ∀l ∈ B (5.74) q π(i)j ≥ 0, o π(i)j ≥ 0, ∀π(i) ∈ Π, ∀j ∈ F (5.75)
The objective function (5.52) seeks to minimize the total cost given by the sum of the travel cost, the fixed charging cost, and the parking cost. Constraints (5.53-5.56) define the battery level when the EV arrives at vertex π(i) ∈ Π if ε π(i)j = 0, or to CS j ∈ F before visiting vertex π(i) ∈ Π, if ε π(i)j = 1. Constraints (5.57 and 5.58) define the battery level when the EV arrives at and departs from CS j ∈ F . Constraints (5.59-5.64) define the battery level in percentage when the EV arrives at or depart from CS j ∈ F . Constraints (5.65) define the time spent at each visited CS. Constraints (5.66) track the arrival time at each vertex. Constraints (5.67) ensure that every vertex is visited within its time window. Constraints (5.68) impose the working schedule of the technician. Constraints (5.69) state that only one CS is visited between any two vertices of the fixed route. Constraints (5.70) ensure that the EV charges only at a visited CSs. Finally, constraints (5.71-5.75) define the domain of the decision variables.

To solve the FRVCP-TW, we use two different solution approaches. The first solves the MILP formulation introduced using a commercial optimizer (we use this approach in the GCI component). The second is a greedy heuristic (we use this approach in the initial solution and the VND component).

Greedy heuristic

We extended here the heuristic proposed in Section 4.3.2 to include the time window constraints. Algorithm 5 describes the structure of our heuristic. It uses five important procedures trackTimeVariables(•), trackBattery(•), sumNegative(•), totalCost(•), and copyAndInsert(•). Procedure track-TimeVariables(•) computes the earliest departure time D π(i) and the latest feasible arrival times L π(i) at each vertex π(i) ∈ Π. Procedure trackBattery(•) computes the battery level Y π(i) at each vertex π(i) ∈ Π. Note that Y π(i) may take negative values. Procedure sumNegative(•) computes the sum of the battery levels with negative values (i.e., sn = π(i)∈Π min{0, Y π(i) }). Procedure totalCost(•) computes the total cost c of the route. Finally, procedure copyAndInsert(•) takes as input a fixed route, a CS, and a position in the route; and returns a copy of the fixed route with the CS inserted at the given position.

The heuristic starts computing D π(i) , L π(i) for each vertex, and the sum sn of the battery levels with negative values for the current fixed-route Π (line 2-7). Then the heuristic enters the outer loop (line 8-38). In each pass through the inner loop (line 9-31), the heuristic: i) computes a slack time st at CS F j between the vertices π(i) and π(i + 1) ∈ Π, and evaluates if the slack time is positive (line [11][12], ii) computes Y and sn, considering that the amount of charge at the inserted CS F j is min{ω(st), ec, Q k }, where ω(st) is the the amount of energy that the EV can charge in st time units, and ec is the energy needed to complete the fixed route (reach the depot) from CS F j , iii) selects the insertion that maximizes sn (lines [18][19][20][21][22]. If sn = 0 (i.e., the route is energy-feasible), the heuristic selects the insertion (CS and position) that minimizes the cost of the route c (line 23-28). If there is at least one feasible insertion (i.e, insertion=true), then the heuristic performs the selected insertion (line 32-34). When no feasible insertion exists, the route can not be repaired and the algorithm stops (line 35-36). If after inserting the CS, the route Π is still energy-infeasible (i.e., sn < 0), the heuristic starts again at line 9 and tries to insert additional CSs until feasibility in terms of energy is reached. for i = 0 to n r -1 do

11: st ←-L π(i+1) -D π(i) -t π(i)F j -t F j ,π (i+1) 12:
if st ≥ 0 then 13:

Π ←-copyAndInsert(Π, F j , i)

14:

Y ←-trackBattery(Π )

15:

sn ←-sumNegative(Y )

16:

c ←-totalCost(Π ) return true, Π 40: end function 5.6 Computational experiments

Test instances for the TRP-CEV

We generate a set of large instances using the data of ENEDIS operation. Maintenance and repair operation of ENEDIS is divided into geographic zones. Each zone has a depot where the technicians start and end their routes. We have 14 test instances classified according to the type geographic zone: four urban, five semi-urban and five rural instances. The number of requests ranges between 54 and 167. The location, time windows, service time, and required skills of the requests come from ENEDIS data. To complete the instances, we added the actual CSs located in the geographic zone of the requests. We obtained the information about of the location and type of technology (i.e., protocol and charging mode) from the open data platform of the French government [START_REF] Quest | Carte des bornes de recharge pour véhicules électriques[END_REF]). To compute the distance and travel time matrices, we use road network data from Open Street Maps (OSM). Furthermore, to build the energy consumption matrices, we use the energy consumption equations model of De Cauwer et al. (2015). As input of those equations, we use fine-grained information on the roads (e.g., distance, speed, elevation) from OSM and NASA, weather information from Metéo France, and real EV parameters (e.g, mass, rolling resistance).

The number of available technicians ranges between 9 and 12. The number of considered skills is two. The skills, working schedule, and lunch break of each technician correspond to real data of the operation. For confidentiality reasons, the cost of this problem was modified but we respect the magnitude and the structure of the ENEDIS operation costs. We consider one type of CVs (Renault Kangoo) and two types of EVs (Renault Kangoo ZE and Peugeot Ion). We took all data related to EVs features from technical sources and specialized websites (Renault 2014, Peugeot 2015, Automobile-propre 2015, [START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF]. In order to better understand the results of the following sections, it is important to know that the travel cost of a Renault Kangoo is much higher than that of a Renault Kangoo ZE and a Peugeot Ion (i.e., 1 vs. 0.08 and 0.06 e/km). For each instance, we generate 5 scenarios changing the percentage of EVs in the fleet (0%, 20%, 40%, 60%, 80% and 100%).

In order to compare the results of our PMa against optimal solutions, we generate a set of small instances with 5 or 10 requests, and 2 or 3 technicians and vehicles. These small instances can be solved using the MILP formulation of the TRP-CEV. To generate each small instance, we randomly selected a large instance, and then randomly drawn from it, requests, technicians, and types of vehicles (ensuring that at least one vehicle is an EV). Our large and small instances are publicly available at www.vrp-rep.org [START_REF] Mendoza | VRP-REP: a vehicle routing community repository[END_REF]).

Parameter settings & experimental environment

The parameters of our PMa is the number of iterations of the GRASP component (i.e., I) and the κ best insertion of the randomized Solomon heuristic. To fine tune I and κ, we conducted a short computational study. For the sake of brevity we do not discuss the study in this chapter. Our results showed that I = 50 and κ = 10 provides the best trade-off between solution quality and computational performance. Although, I = 50 may seem as a small number of iterations, notice that the GRASP component is used to feed the pool of routes Ω for each couple < k, u >. A large pool of routes impacts the performance of the solution of the SC formulation.

We implemented our PMa in Java (jre V.1.8.0) and used Gurobi Optimizer (version 6.0) to solve the FRVCP-TW and the SC formulation. We set a time limit of 800 seconds on Gurobi to control the running time of the SC. All the experiments were run on a computing cluster with 2.33 GHz Intel Xeon E5410 processors with 16 GB of RAM running under Linux Rocks 6.1.1.

Performance of PMa

Since the TRP-CEV is a new problem, there are no results or algorithms to benchmark against. To get an idea of the quality of the solutions delivered by our PMa, we conducted different experiments. We ran our PMa on the small instances and compared their results with optimal solutions found using the MILP formulation. Furthermore, we ran our PMa on a set of 14 private industrial-scale instances of ENEDIS and compared our results with solutions found by the routing tool used by the company. Finally, we ran our PMa on a set of instances of a problem related with the TRP-CEV, and compared our results with the results obtained by the approach proposed in the literature.

Results on small TRP-CEV instances Table 5.1 compares the results of our PMa with the optimal solution on the small instances of the TRP-CEV. In this table, we compare the results under the light of four metrics: the number of optimal solutions found, the average gap with respect to the optimal solution over the 10 runs, the average best gap3 , and the average computing time. The results of this table suggest that PMa is able to deliver high-quality TRP-CEV solutions for small instances. PMa matched the 10 optimal solutions and the values reported for Avg. Gap shows that our approach is stable for those instances. Appendix A.3 presents detailed results for the 10 instances. Results on private industrial-scale instances Table 5.2 compares the results of our PMa with the results of the professional routing tool currently deployed at ENEDIS on the set of large instances. It is worth noting that their tool cannot deal with EVs. Therefore, the instances used in this section only consider CVs. Table 2 uses the same metrics of Table 1, but the number of optimal solutions is replaced by the number of best known solutions (BKSs) found by the routing software or by our PMa. Since we only have one solution of the routing tool for each instance, we only compute the average gap for the routing software4 . The results in Table 5.2 show that our PMa consistently outperforms the current routing tool used by the company. It not only found the 14 BKSs but also reported a significantly lower average gap (2.71% vs. 6.56% ). 2016) assume that customers have a positive demand, the EVs have limited hauling capacity, all CSs have the same technology, the EVs charge fully their battery at CSs, the charging function is linear, the consumption rate is constant, and the EVs have a fixed cost. The objective function of E-FSMFTW seeks to minimize the sum of the cost of the fixed cost of used vehicles plus the total travelled distance.

Results for the E-FSMFTW instances

We adapted our PMa to solve the E-FSMFTW instances. It is worth noting that we do not aim to establish new state-of-the-art results for the E-FSMFTW. We rather want to show that our PMa can also be competitive on a related problem with only mild modifications.

Hiermann et al. ( 2016) proposed a set of 108 small instances with the number of customers ranging from 5 to 15, and a set of 168 large instances with 100 customers and 21 CSs. To solve those instances, [START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF] propose an adaptive large neighbourhood search (ALNS) and a branch-and-price (BnP). Their BnP was able to obtain the optimal solution within the time limit of two hours for all the small instances.

Tables 5.3 and 5.4 summarize the results delivered by our PMa and their ALNS on the small and large E-FSMFTW instances. For the small instances, we compare the results obtained by PMa and ALNS with the optimal solutions obtained by BnP. For the large instances, we compare the results obtained by PMa and ALNS with the BKSs taken from Hiermann et al. ( 2016) 5 and updated with some new BKSs found by our PMa. The results of the PMa and ALNS are computed over 10 runs. It is important to mention, that Hiermann et al. (2016) only report the average cost of the solutions that they obtained for each large instances. Appendix A.4 reports the detailed results of this experiment. On the small instances, PMa delivered competitive results with respect to ALNS. Although, PMa just matched 81 optimal solutions, its average gap is smaller than the achieved by ALNS (0.32% vs 0.55%). In terms of CPU time, PMa is faster than the ALNS.

The results for large instances show that PMa obtains competitive results, especially considering that it was not originally conceived to tackle the E-FSMFTW (Table 5.4). Our PMa matched 69/168 BKSs; 61 are new BKSs6 . However, the PMa has a higher average gap (1.18% vs 2.20%). After carefully analyzing the results, we believe this behavior may be explained by the characteristic of the Hiermann et al. (2016) instances. Indeed, in their instances, since some of the EVs have an autonomy of about 50 km (which is much lower than in reality). In consequence, it is possible that in order to reach some isolate customers it may be interesting (from a cost point of view) to visit more than one CSs between two customers. Now, this possibility is ruled out in our PMa. Note that in the definition of our FRVCP-TW visiting more than one CS between two customers (hereafter called multiple consecutive visits) is forbidden. This was a deliberated design decision motivated by three interrelated reasons. First, forbidding multiple consecutive visits vastly simplifies and speeds up solving the FRVCP-TW (either with the heuristic or the solver). Second, we conceived our algorithm to work on a problem where there is fixed cost for charging that largely offsets the driving cost of an EV. Therefore, multiple consecutive visits are naturally discourage. Third, our algorithm is meant to be used in a real-world context, where EVs have larger autonomies and, therefore, are less likely to require multiple consecutive visits to reach a customer. We believe that the trade off between missing some high-quality routes during the search and having a faster and simpler algorithm is positive.

Managerial insight

ENEDIS is starting the migration from a fleet of CVs to a mixed fleet of CVs + EVs. This process raises a number of managerial questions regarding the fleet composition and the best strategies to exploit EVs. To help managers answering those questions, we conducted a study analyzing the impact of the fleet composition and the access to charging infrastructure in the feasibility and quality of TRP-CEV solutions.

To conduct this study, we used the 84(= 14 × 6) instances described in Section 5.6.1. For each instance we ran our PMa 10 times and report average results. We analyze the data under the light of 4 metrics that are important to ENEDIS managers: the total cost (and how it is split between fixed and variable costs), the total CVs emissions7 , the number of visited CSs, and the number of routes.

Fleet composition analysis

We analyse in this section the behavior of the metrics when the proportion of EVs in the fleet changes. Figure 5.1 shows of each fleet composition: the average fixed, variable and total cost. The symbol N.F. indicates that in none of the 10 runs, our PMa could find a feasible solutions for the instance under the given fleet composition. The latter is the case for every rural instance when the fleet is entirely composed of EVs. Notice that a pathological case is observed on instance rural_21. A close look at the instance data reveals that the requests are particularly scattered over a large geographical region. In consequence, EVs need to charge their batteries quite often. Since charging the battery is a time-consuming process, the resulting routes are long (in terms of travel time) and difficult (if not impossible) to conciliate with the technicians' availability schedules.

The data on the fixed costs reported in Figure 5.1 shows an interesting and somehow unexpected behavior. It is worth reminding that in our problem the fixed costs are directly (and exclusively) related to the number of routes. Intuition says that on instances in which the fleet is composed mainly of EVs more routes would be needed to service requests. The reason is that routes covered with EVs tend to serve fewer requests than routes covered with CVs (because the EV spends a portion of their available time on charging operations). The data, however, shows that (at least on our instances) the fixed cost remain constant across the different fleet compositions.

A more intuitive result comes from the data on the variable costs. As Figure 5.1 shows, the variable costs decrease when the proportion of EVs in the fleet composition increases. This result is largely explained by the considerably higher cost per kilometre allocated to CVs with respect to EVs. It is worth noting that variable costs are higher on rural instances than in semi-urban and urban instances. The reason is simply that in rural operations vehicles travel larger distances. For instances, the average travelled distance per served request is 9.52 km on rural instances compared to 3.74 km and 1.94 km for the semi-urban and urban instances, respectively.

Access to charging infrastructure

Up to today, ENEDIS does not own charging infrastructure outside their depots. Therefore, the midroute recharges take place in public CSs. For a number of reasons ranging from difficult access to high charging costs, ENEDIS would prefer to avoid mid-route charging at public stations. To help them evaluate the impact of disallowing mid-route charging outside the depot on the quality and feasibility of TRP-CEV

Conclusion

In this chapter, we introduce what we call the technician routing problem with conventional and electric vehicles (TRP-CEV). The TRP-CEV is largely inspired by a problem faced by ENEDIS in their daily operations. This problem consists in routing a set of technicians to serve a set of requests, using a fleet composed by EVs and CVs of different types. This problem considers several real-life constraints such as: time windows on the requests; technician skills, working schedules and lunch breaks; skills incompatibility 5.8.4 Notation for the FRVCP-TW Π: Route that violates the energy constraint. Π: Route that includes the dummy vertex associated to the lunch break. ε π(i)j : Binary variable, equal to 1 if the vehicle charges at CS j ∈ F before visiting vertex π(i) ∈ Π. If π(i)j = 1, φ π(i) is the battery level when the vehicle arrives at CS j ∈ F right before visiting vertex π(i). φ π(i) : Battery level if ε π(i)j = 0, φ π(i) is the battery level when the EV arrives at vertex π(i). α π(i)jl : Binary variable, equal to 1 if the battery level is l ∈ B when the vehicle arrives at j ∈ F before visiting vertex π(i) ∈ Π. β π(i)jl : Binary variable, equal to 1 if the battery level is l ∈ B when the vehicle arrives at and departs from CS j ∈ F before visiting vertex π(i) ∈ Π. γ π(i)j : Energy charged at CS j ∈ F before visiting vertex π(i) ∈ Π. q π(i)j : Battery level when the EV arrives CS j ∈ F before visiting vertex π(i) ∈ Π. o π(i)j : battery level when the EV departs from CS j ∈ F before visiting vertex π(i) ∈ Π δ π(i)j : Spent time at CS j ∈ F before visiting vertex π(i) ∈ Π. D π(i) : Earliest departure time at vertex π(i) ∈ Π. L π(i) : Latest feasible arrival time at vertex π(i) ∈ Π. Y π(i) : Battery level at vertex π(i) ∈ Π. c: Total cost of the route. sn: Sum of the battery levels with negative values.

Chapter 6

General conclusions & perspectives

This final chapter summarizes the contributions of the thesis and offers future research directions. This thesis is devoted to study the electric vehicle routing problems (eVRPs). These problems extend classical VRPs to consider the main features of electric vehicles (EVs). Despite some researchers have started to study these kind of problems, the gap with classical models and, most importantly, real-world applications is still significant. Therefore, the purpose of this thesis was to contribute to closing that gap.

In Chapter 2, we proposed a simple yet effective two-phase heuristic to tackle the Green VRP. In the first phase our heuristic builds a pool of routes via a set of randomized route-first cluster-second heuristics with an optimal station insertion procedure. In the second phase our approach assembles a Green VRP solution by solving a set partitioning formulation over the columns (routes) stored in the pool. To test our approach, we performed experiments on a set of 52 instances from the literature. The results showed that our heuristic is competitive with state-of-the-art methods. Our heuristic unveiled 8 new best-known solutions, matched another 40, and delivered solutions with an average gap of 0.14% for the 4 remaining instances.

The Green VRP and most of the existing eVRP studied in the literature rely on one (or both) of the following assumptions: (i) the vehicles fully charge their batteries every time they reach a charging station (CS), and (ii) the battery charge level is a linear function of the charging time. In practical situations, however, the amount of charge is a decision variable, and the battery charge level is a concave function of the charging time. In Chapter 3, we extended current eVRP models to consider partial charging and nonlinear charging functions. We presented a computational study comparing our assumptions with those commonly made in the literature. Our results suggested that neglecting partial and nonlinear charging may lead to infeasible or overly expensive solutions To tackle large-scale instances of the eVRP with nonlinear charging function (eVRP-NL), in Chapter 4 we proposed an iterated local search (ILS) enhanced with heuristic concentration (HC). At each iteration, our ILS applies a variable neighborhood descent (VND) procedure with three different local search operators and stores the routes of the local optimum in a pool of routes. When the ILS finishes, the HC component assembles the final solution by recombining routes from the pool. One original feature of our approach is the use in the VND component of a local search operator that aims at improving the charging decisions of each route (where and how much to charge). To achieve this goal, the operator solves a new optimization problem: the fixed-route vehicle-charging problem (FRVCP). Our implementation solves the FRVCP using a heuristic and a commercial solver. To test our approach, we proposed a 120-instance testbed for the e-VRP-PNL . We discussed experiments assessing the performance of the proposed approach and showing the relevance of making optimal charging decisions.

Finally, in Chapter 5, we introduced the technician routing problem with conventional and electric vehicles (TRP-CEV), which is inspired in the real operation of ENEDIS . This problem consists in routing a set of technicians to serve a set of geographically scattered requests within their time windows. This problem considers the working schedules, the lunch breaks, and skills of the technicians. The technicians use a fixed fleet composed of conventional and electric vehicles. The EVs have a limited autonomy, and there are incompatibility constraints between EVs and CSs. In the TRP-CEV decisions include the vehicle-CHAPTER 6. GENERAL CONCLUSIONS & PERSPECTIVES to-technician assignment, the sequencing of the routes followed by the technicians , and the battery charging program for the EVs (where and how much to charge) considering that the charging functions are nonlinear. The objective function of TRP-CEV seeks to minimize the total cost, defined as the sum of the travel costs, battery charging fixed costs, the parking cost at CSs, and the fixed cost of using each technician. To tackle the TRP-CEV, we proposed a two-phase parallel matheuristic (PMa). In the first phase, the PMa builds a pool of feasible TRP-CEV routes solving in parallel a set of sub-problems. In the second phase the PMa assembles a TRP-CEV solution by solving an extended set covering formulation. To test our approach, we ran it on different benchmarks. The experiments showed that our PMa has competitive results. Furthermore, we applied (tested) the approach on a set of real instances based on the ENEDIS operation. We analysed the obtained solutions by our PMa on those instances, in order to evaluate the quality and structure of the solution as a function of the percentage of EVs in the fleet.

As a general conclusion, in this thesis we have highlighted the importance to make optimal charging decision in the eVRPs. We have found that the optimal charging decisions are essential to have competitive results in terms of solution quality. Therefore, for all the eVRPs addressed in this thesis, we proposed different approaches for making optimal charging decisions. On the other hand, in a large part of this thesis, we included realistic assumptions about the battery charging process. Those assumptions are the partial charging and the non-linear charging function. The incorporation of those assumptions into the eVRPs and the proposed approaches to tackle these assumptions allow to model more realistic eVRP problems. Other notable conclusion for the context of the eVRP is that the more complicated eVRPs (i.e., where there are charging decisions) are in rural zones. In different experiments conducted during this thesis, we found that only in rural areas there are day-long routes, where the travelled distance is larger than the EV autonomy.

To further our research we plan to design an exact method for the TRP-CEV to have a better base of comparison for our results. Recently, researchers have proposed exact methods for problems directly related to the TRP-CEV. For instance, a branch-and-price method for the technician dispatching problem [START_REF] Cortés | Branch-and-price and constraint programming for solving a real-life technician dispatching problem[END_REF]), a branch-price-and-cut method for the eVRPs with time windows [START_REF] Desaulniers | Exact algorithms for electric vehicle-routing problems with time windows[END_REF]) and a branch-and-price method for a VRP with lunch break [START_REF] Liu | Mathematical model and exact algorithm for the home care worker scheduling and routing problem with lunch break requirements[END_REF]). Taking advantage of recent studies, the idea behind this perspective is to integrate and extend those last advances to propose an exact method for the TRP-CEV.

Despite this thesis contributes to bringing the eVRP closer to reality in terms of the battery charging assumptions, there is a pending research about other assumptions of the eVRPs. We suggest that further research should be undertaken in the following assumptions: the capacity of the CSs, and the uncertain CS availability and energy consumption. In terms of the methods, an interesting perspective for the eVRPs context is to explore alternative approach for the fixed-route charging-problems that ensures the optimal solution with computing time faster. Despite we tried -without success-to find a polynomial time algorithm, we still think that it is possible to find a better approach than the proposed in this thesis. Considering that in the road network the consumption of energy, the travelled distance, and the travel time are not proportional, another interesting perspective is to model the eVRP as a multi-graph. The idea behind this perspective is that existing more than one edge between two nodes, where each edge represents attributes as the energy consumption, the travel time, or the travelled distance.

Résumé

Etant donné leur faible impact environnemental, l'utilisation des véhicules électriques dans les activités de service a beaucoup augmenté depuis quelques années. Cependant, leur déploiement est freiné par des contraintes techniques telles qu'une autonomie limitée et de longs temps de charge des batteries. La prise en compte de ces contraintes a mené à l'apparition de nouveaux problèmes de tournées de véhicules pour lesquels, en plus d'organiser les tournées, il faut décider où et de combien charger les batteries. Dans cette thèse nous nous intéressons à ces problèmes au travers de quatre études. La première concerne le développement d'une metaheuristique en deux phases simple mais performante pour résoudre un problème particulier appelé "Green VRP". Dans la seconde, nous nous concentrons sur la modélisation d'un aspect essentiel dans ces problèmes : le processus de chargement des batteries. Nous étudions différentes stratégies pour modéliser ce processus et montrons l'importance de considérer la nature non linéaire des fonctions de chargement. Dans la troisième étude nous proposons une recherche locale itérative pour résoudre des problèmes avec des fonctions de chargement non linéaires. Nous introduisons un voisinage dédié aux décisions de chargement basé sur un nouveau problème de chargement sur une tournée fixée. Dans la dernière étude, nous traitons un problème réel de tournées de techniciens avec des véhicules classiques et électriques. Ce problème est résolu par une metaheuristique qui décompose le problème en plusieurs sous-problèmes plus simples résolus en parallèle, puis qui assemble des parties des solutions trouvées pour construire la solution finale.

Abstract

Electric vehicles (EVs) are one of the most promising technologies to reduce the greenhouse gas emissions. For this reason, the use of EVs in service operations has dramatically increased in recent years. Despite their environmental benefits, EVs still face technical constraints such as short autonomy and long charging times. Taking into account these constraints when planning EV operations leads to a new breed of vehicle routing problems (VRPs), known as electric VRPs (eVRPs). In addition, to the standard routing decisions, eVRPs are concerned with charging decisions: where and how much to charge. In this Ph. D thesis, we address eVRPs through 4 different studies. In the first study, we tackle the green vehicle routing problem. To solve the problem, we propose a simple, yet effective, two-phase matheuristic. In the second study, we focus a key modelling aspects in eVRPs: the battery charging process. We study different strategies to model this process and show the importance of considering the nonlinear nature of the battery charging functions. In the third study, we propose an iterated local search to tackle eVRP with non-linear charging functions. We introduce a particular local search operator for the charging decisions based on a new fixedroute charging problem. The fourth and last study considers a real technician routing problem with conventional and electric vehicles (TRP-CEV). To tackle this problem, we propose a parallel matheuristic that decomposes the problem into a set of easier-to-solve subproblems that are solved in parallel processors. Then the approach uses parts of the solutions found to the subproblems to assemble final solution to the TRP-CEV.

Mots clés

Problèmes de tournées de véhicules électriques ; Problèmes de tournées de techniciens ; Green vehicle routing problems ; Problème de chargement sur une tournée fixée ; Fonctions de chargement non linéaires ; Méta-heuristique ; Matheuristique, Recherche locale itérative

Chapter 2 A

 2 multi-space sampling heuristic for the green vehicle routing problem The work presented in this chapter has been published in the Transportation Research Part C: Emerging Technologies Journal. The full reference follows: -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2015). A multi-space sampling heuristic for the green vehicle routing problem. Transportation Research Part C: Emerging Technologies, 2016, 70, 113 -128.

Figure 2

 2 Figure 2.1 -Example of a feasible Green VRP solution

  that transfer between routes sequences of customers of arbitrary length. Furthermore, Schneider et al. (2014) improved the MILP formulation proposed by Erdogan & Miller-Hooks (2012), and they evaluated their VNS/TS and MILP approaches on the 52-instance testbed proposed by Erdogan & Miller-Hooks (2012). Computational experiments showed that CPLEX 12.2 was unable to solve to optimality instances with 20 customers using their MILP model, and VNS/TS outperformed the constructive heuristics proposed by Erdogan & Miller-Hooks (2012). Schneider et al. (2015) introduced the vehicle routing problem with intermediate stops (VRPIS) that generalizes the Green VRP and the MDVRPI. To solve the VRPIS Schneider et al. (

  Felipe et al. (2014) tested their approaches on the instances of Erdogan & Miller-Hooks (2012). Their methods outperformed the constructive heuristics proposed by Erdogan & Miller-Hooks (2012), but are not competitive with the solution approaches of Schneider et al. (2014) and Schneider et al.

Figure 2 . 2 - 2

 222 Figure 2.2 -Splitting a TSP tour into a Green VRP solution

Figure 2 . 3 -

 23 Figure 2.3 -Outline of the structure of the repair graph B = (Z, U )

Figure 2 . 4 -

 24 Figure 2.4 -Optimal repair example for the three-customer sequence r = {0, A, B, C, 0}, three AFSs.

  , Erdogan & Miller-Hooks (2012), Schneider et al. (2014), Schneider et al. (2015), and Felipe et al. (

Figure 2 . 5 -

 25 Figure 2.5 -Trade-off between solution quality and CPU time

Figure 2 .

 2 Figure 2.6 -Percentage share of CPU time by MSH phases

2. 6 . 1

 61 Notation for problem definition G: Green VRP underlying graph G = (V, E). V : Set of vertices of G. 0: Depot I: Subset of V representing customers. F : Subset of V representing alternative fuel stations. E: Set of edges in G. d ij : Distance between vertices i and j (i, j ∈ V ). t ij : Travel time between vertices i and j (i, j ∈ V ). n: Number of customers. a: Number of alternative fuel stations. T max : Tour duration constraint. cr: Fuel consumption rate.

3 A

 3 [i, k]: Vertex of set Z, which represents visit to station k ∈ F after visit to vertex r i (i.e., ith element in route Π). c: Arc cost of repair graph. t: Arc time of repair graph.Chapter comparative study of charging assumptions in electric vehicle routing problemsThe work presented in this chapter has been submitted in the Transportation Research Part B: Methodological Journal. The full reference follows: -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2016). A hybrid metaheuristic for the electric vehicle routing problem with partial charging and nonlinear charging function. Transportation Research Part B: Methodological. Under second round review. The work of this chapter was presented at the following conferences and seminars: -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2015). The electric vehicle routing problem with non-linear charging functions. In ODYSSEUS 2015, Ajaccio, France. -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2016). A comparative study of charging assumptions in electric vehicle routing problems. In ROADEF 2016, Compiègne, France. -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2016). A comparative study of charging assumptions in electric vehicle routing problems. In EURO 2016, Poznan, Poland. -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2016). Contributions to the electric vehicle routing problems. PGMO Seminars, L'École polytechnique, Paris, France, March 18, 2016 -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2016). Electric vehicle routing problems with non-linear charging functions. Cirrelt seminar, Université de Montréal Campus, Montreal, Canada, October 18, 2016

Figure 3 .

 3 Figure 3.1 -Typical charging curve, where i and u represent the current and terminal voltage respectively. (Source Hõimoja et al. (2012))

FirstFigure 3 . 2 -

 32 Figure 3.2 -First segment approximation (FS)

Figure 3 . 3 -

 33 Figure 3.3 -Linear approximations of charging functions.

Figure 3 . 4 -

 34 Figure 3.4 -Approximation vs real data

2 Charging 2 ChargingFigure 3 . 5 -

 2235 Figure 3.5 -Example of a feasible eVRP-NL solution

  Figure 3.6 -Piecewise linear approximation for the charging function.

  -3.22) define the charge level (and its corresponding charging time) when an EV departs from CS i ∈ F . Constraints (3.23) define the time spent at any CS. Constraints (3.24) and (3.25) track the departure time at each vertex, where S max = max i∈F {c ib }. Constraints (3.26) and (3.27) ensure that the EVs return to the depot no later than T max . Constraints (3.28) and (3.31) help to avoid the symmetry generated by the copies of the CSs. The parameter m ij is equal to 1 if i and j ∈ F represent the same CS. Finally, constraints (3.32-3.35) define the domain of the decision variables.

  Conrad & Figliozzi 2011, Erdogan & Miller-Hooks 2012, Schneider et al. 2014, Sassi et al. 2014, Goeke & Schneider 2015, Hiermann et al.

3. 6 . 1 4 A

 614 Notation for the MILP formulation Chapter hybrid metaheuristic for the electric vehicle routing problem with partial charging and nonlinear charging function The work presented in this chapter has been submitted in the Transportation Research Part B: Methodological Journal. The full reference follows: -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2016). A hybrid metaheuristic for the electric vehicle routing problem with partial charging and nonlinear charging function. Transportation Research Part B: Methodological. Under second round review. The work of this chapter was presented at the following conferences and seminars: -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2015). The electric vehicle routing problem with non-linear charging functions. In ODYSSEUS 2015, Ajaccio, France. -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2016). Contributions to the electric vehicle routing problems. PGMO Seminars, L'École polytechnique, Paris, France, March 18, 2016 -Montoya, A., Guéret, C., Mendoza, J. E., and Villegas, J. G. (2016). Electric vehicle routing problems with non-linear charging functions. Cirrelt seminar, Université de Montréal Campus, Montreal, Canada, October 18, 2016

  al. (1974)

Figure 4 .

 4 Figure 4.1 -General structure of ILS+HC

Figure 4 . 2 -

 42 Figure 4.2 -Piecewise linear charging function and fixed-route for the FRVCP

3 )

 3 φ π(nr) = φ π(nr-1) + j∈F k∈B\{0} δ π(nr-1)jk + j∈F k∈B\{0} δ π(nr)jk -j∈F ε π(nr-1)j e jπ(nr-1) -j∈F ε π(nr)j (e π(nr-1)j + e jπ(nr) ) -e π(nr-1)π(nr) nr-1)jk -j∈F e jπ(nr-1) ε π(nr-1)j -j∈F e π(nr-1)j ε π(nr)j ≥ 0 (4.5)

Algorithm 3

 3 Greedy heuristic 1: function GREEDYHEURISTIC(Π 0 ,F ) for j = 1 to |F | do 9:for i = 0 to n r -1 do 10:Π ←-copyAndInsert(Π, F j , i)

  copyAndInsert(Π, F u , i)

Figure 4 . 4 -

 44 Figure 4.4 -Percentage of the routes with/without visits to CSs by instance size.

Figure 4

 4 Histogram of number of visits to CS per route. visits to CSs per route (b) Histogram of the maximum number of visited CSs in the routes of each solution.

Figure 4 . 5 -

 45 Figure 4.5 -Analysis of the number of visits to CSs

  al. 2016, Baldacci et al. 2008). However, only recently researchers have started to study the HVRP with EVs. Goeke & Schneider (2015) introduced the electric vehicle routing problem with time windows and mixed fleet (E-VRPTWMF), which considers a fleet of CVs and EVs. Sassi et al. (2015) extended the E-VRPTWMF considering different features motivated by a real-life industrial application. In their problem, the CSs have different technologies and time windows, the charging cost depends on the time of the day, and EVs can charge partially their batteries. Finally, Hiermann et al. (2016) introduced the electric fleet size and mix vehicle routing problem with time windows and recharging stations (E-FSMFTW), where the fleet is composed of different types of EVs. The objective of the E-FSMFTW is to select the size and composition of the EV fleet. Although, other previous studies have assumed mixed fleet with EVs and CVs, the TRP-CEV is the only to assume the nonlinear charging behavior of the charging process.

Algorithm 4

 4 Parallel matheuristic: General structure 1: function PARALLELMATHEURISTIC(G,T ,K) 2:

  , π(n r )} (5.54) φ π(nr) = φ π(nr-1) + j∈F γ π(nr-1)j + j∈F γ π(nr)j -j∈F e jπ(nr-1)k • ε π(nr-1)j -j∈F e π(nr-1)jk + e jπ(nr)k • ε π(nr)je π(nr-1)π(nr)k • 1 -j∈F ε π(nr)j (5.55) φ π(nr-1) + j∈F γ π(nr-1)j -j∈F e jπ(nr-1)k • ε π(nr-1)j -j∈F e π(nr-1)jk • ε π(nr)j ≥ 0 (5.56)

  Hiermann et al. (2016) recently introduced the closely-related E-FSMFTW. This problem consists in routing an unlimited heterogeneous fleet of EVs to serve a set of geographically scattered requests within their time windows. For this problem Hiermann et al. (

  . N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F. N.F.

Figure 5 .

 5 Figure 5.1 -Average fixed (i.e., the fixed cost of each technician), variable (i.e., the sum of the total travel cost, fixed charging cost and parking cost), and total cost for each instance for each fleet composition.

Figure 5 .

 5 Figure 5.2a shows the Kg of CO 2 emitted (on average) by the CVs employed in the solutions to each instance for each fleet composition. As expected, the emissions decrease when the proportion of EVs in the fleet increases. Notice that rural instances represent the most significant potential for emission reductions because vehicles travel longer distances. Finally, Figure 5.2b shows the maximum number of visited CSs in a solution over the 10 runs of each instance for each fleet composition. Note that only the solutions of the rural instances include routes visiting CSs. This behavior is expected because requests are geographically dispersed in those instances. This result confirms eVRPs arise mainly in the rural and semi-urban operations.

Figure 5 . 3 -

 53 Figure 5.3 -Average gap between the cost of the solutions with and without the option of visiting the CSs

Figure 5 . 4 -

 54 Figure 5.4 -Factors explaining the increment on the cost when the visits to the CSs are forbidden

  then

	14:	s * ←-s
	15:	end if
	16:	end if
	17:	k ←-k + 1
	18:	end for
	19:	end while
	20:	R ←-SetPartitioning(G, Ω, s * )
	21:	return R
	22: end function

Table 2 .

 2 1 -Summary results and comparison of our MSH with other methods on the small instances of[START_REF] Erdogan | A green vehicle routing problem[END_REF].

	Metric	MCWS	VNS	AVNS	48A	MSH(1k) MSH(5k) MSH(10k)
		/DBCA	/TS
	Number of BKS	2/40	38/40	40/40	29/40	35/40	40/40	40/40
	Avg. Gap (%)	NR	NR	0.15%	NR	0.38%	0.04%	0.01%
	Max. Gap (%)	NR	NR	1.73%	NR	3.61%	0.47%	0.13%
	Avg. Best Gap (%)	8.72%	0.63%	0.00%	0.46%	0.09%	0.00%	0.00%
	Cum. Number of Veh.	245	223	NR	225	222	222	222
	Avg. Time (min)	NR	0.65	0.17	0.02	0.01	0.04	0.07
	Max. Avg. time (min)	NR	0.88	0.38	0.04	0.02	0.06	0.12
	Computer	Pentium 4 Core I5 3.2GHz 2.67GHz 2.67GHz 2.8GHz Core I5 Core I5		XEON E5410 2.33 GHz
	Runs	NR	10	10	1	10	10	10
	NR: not reported

Table 2 .

 2 2 -Summary results and comparison of our MSH with other methods on the large instances of Erdogan & Miller-Hooks (2012).

	Metric		MCWS	VNS	AVNS	48A	SA	MSH	MSH	MSH
				/DBCA	/TS				(1k)	(5k)	(10k)
	Number of BKS	0/12		0/12	4/12	0/12	0/12	0/12	0/12	8/12
	Avg. Gap (%)	NR		NR	0.92%	NR	NR	2.64% 1.48% 1.02%
	Max. Gap (%)	NR		NR	1.84%	NR	NR	5.77% 4.21% 3.62%
	Avg. Best Gap (%)	15.97%	1.38%	0.17%	4.50%	4.97% 1.42% 0.40% 0.05%
	Cum. Number of Veh.	508		461	NR	466	459	454	445	444
	Avg. Time (min)	NR		159.58	6.20	157.03 156.05	27.92	31.47	35.04
	Max. Avg. time (min)	NR		525.52	19.51	514.68 456.26	80.11	84.95	89.95
	Computer		Pentium 4 Core I5 3.2GHz 2.67GHz 2.67GHz 2.8GHz 2.8GHz Core I5 Core I5 Core I5	XEON E5410 2.33GHz
	Runs		NR		10	10	1		1	10	10	10
	NR: not reported
		100.000%							100.00%
	gap (%)	1.000% 10.000%	48A			VNS/TS		gap (%)	10.00%	SA	48A
	Average best	0.001% 0.010% 0.100%	MSH(1k)					Average best	0.10% 1.00%	AVNS	MSH(1k) MSH(10k) MSH(5k)	VNS/TS
			MSH(5k) MSH(10k)	AVNS
		0.000%							0.01%
		0.01	0.1		1			1.00	10.00	100.00	1000.00
			Average time (min)				Average time (min)
			(a) Small instances (%)

(b) Large instances

  al. 2014, Desaulniers et al. 2014, Hiermann et al. 2016, Keskin & Catay 2016

Table 3 .

 3 1 -Comparison of our charging assumptions with charging assumptions from the literature

	Instance	PC-PL	FC-PL	PC-FS	PC-L1		PC-L2
				Solution	Evaluation	Solution	Evaluation

  al. 2015, Hiermann et al. 2016, Liao et al. 2016). To our knowledge, only Sweda et al. (

Table 4 .

 4 1 -Comparison of the two versions of the metaheuristic on small instances with proven optima

		18.00
		16.00
	(kWh)	12.00 14.00
	level	10.00
	Battery	6.00 8.00
		4.00
		2.00
		0.00
		0.00	0.50	1.00	1.50		2.00	2.50
				Time (hours)
			Type of CSs		Slow	Moderate	Fast
			Charging power (kWh/h)	11	22	44
			Piecewise linear approximation
	Figure 4.3 -Piecewise linear approximation for different types of CS charging an EV with a battery of 16
	kWh.
		Metric			ILS(H)+HC ILS(S)+HC
		Number of optimal solutions	15/20	20/20
		Avg. Gap (%)		1.97	0.34
		Max. Gap (%)		17.07	1.87
		Avg. Best Gap (%)		1.20	0.00

Gurobi to control the running time of the HC phase. All the experiments were run on a computing cluster with 2.33 GHz Inter Xeon E5410 processors with 16 GB of RAM running under Linux Rocks 6.1.1. The results delivered by the two ILS+HC configurations are compared over 10 runs. Each replication of the experiments was run on a single processor.

Table 4 .

 4 2 -Comparison of the two versions of the metaheuristic on large instances

	Metric	Gurobi ILS(H)+HC ILS(S)+HC
	Number of solutions 25/100	100/100	100/100
	Number of BKSs	7/100	7/100	100/100
	Avg. Gap (%)	41.36	7.51	1.51
	Max. Gap (%)	NA	28.68	4.44
	Avg. Best Gap (%)	NA	5.28	0.00
	more stable behavior than ILS(H)+HC.

Table 4 .

 4 3 -Average computing time (in seconds) of different variants of the metaheuristic

	Instance size ILS(H)+HC ILS(S )+HC ILS(S)+HC		Speedup
	10	0.64		8.54	5.62		1.52
	20	1.75		17.47	10.56		1.65
	40	8.48		64.16	35.35		1.82
	80	39.35	148.76	80.11		1.86
	160	289.08	976.84	568.02		1.72
	320	2,568.94	5,759.67	4,397.64		1.31
	Average	484.71	1,162.57	849.55		1.65
	Max	4,766.36 10,335.56	7,636.50		2.53
	Min	0.49		3.71	2.36		1.16
	ILS(S)+HC and ILS(S )+HC are the hybrid metaheuristic with and without preprocessing strategies, respectively
	100%
	90%
	80%
	70%
	60%
	50%
	40%
	30%
	20%
	10%
	0%
	10	20	40		80	160	320
			Number of customers
		Routes with CSs	Routes without CSs

Table 5 .

 5 1 -Comparison of the PMa on small instances with proven optima

	Metric	PMa
	Number of optimal solutions 10/10
	Avg. Gap (%)	0.00
	Avg. Best Gap (%)	0.00
	Avg. Time (s)	0.55

Table 5 .

 5 2 -Comparison of the PMa with the routing software used by ENEDIS

	Metric	Routing software PMa
	Number of BKSs	0/14 14/14
	Avg. Gap (%)	6.56 2.71
	Avg. Best Gap (%)	-0.00
	Avg. Time (min)	-3.39

Table 5 .

 5 3 -Comparison of our PMa with the ALNS by Hiermann et al. (2016) on small instances

	Metric	ALNS	PMa
	Number of optimal solutions	108/108	81/108
	Avg. Gap (%)	0.55	0.32
	Avg. Best Gap (%)	0.00	0.31
	Avg. Time (min)	0.32	0.06
	Computer	Core2 Quad XEON E5410 with 2.4 GHz with 2.33 GHz

Table 5 .

 5 4 -Comparison of our PMa with the ALNS by Hiermann et al. (2016) on large instances

	Metric	ALNS	PMa
	Number of BKS	NR	69
	Avg. Gap (%)	1.18%	2.20%
	Avg. Best Gap (%)	NR	1.52%
	Avg. Time (min)	22.66	25.17
	Computer	Intel Core2 Quad CPU Q6600 with 2.4 GHz with 2.33 GHz XEON E5410

Table A .

 A 6 -Results of PMa on large instances of TRP-CEV .31 20.02 4158.41 4260.54 17.93 4026.93 4599.33 19.62 4028.21 4745.61 21.18 4024.50 4585.72 18.70 4024.97 4345.65 20.14 semi_urbain_19 3395.19 3400.45 19.51 3278.80 3285.17 19.24 3212.82 3213.21 19.34 3212.06 3212.36 19.98 3212.09 3212.39 20.85 3212.08 3212.28 20.28 semi_urbain_20 6898.40 7598.07 18.17 7503.12 7522.85 20.03 6582.82 7288.81 19.26 7275.06 7285.98 18.11 7233.03 7235.50 18.55 7231.85 7234.76 19.63 semi_urbain_21 6899.58 7537.84 22.18 7490.38 7522.74 19.38 6598.19 7264.17 20.31 7265.63 7274.67 19.14 6436.39 7078.40 19.05 7232.37 7234.52 18.63 semi_urbain_22 6096.59 6132.66 19.95 5874.59 5890.09 20.20 5683.24 6260.91 18.19 5631.99 6125.93 21.51 5634.75 6114.44 18.59 5630.73 6032.78 20.85 urbain_18 3355.42 3362.49 19.93 3331.72 3359.44 23.86 3211.44 3227.53 20.73 3210.76 3211.37 19.11 3210.76 3211.37 21.05 3210.76 3211.37 18.00 urbain_19 4104.48 4111.86 19.12 4090.84 4111.68 19.31 4055.55 4614.59 20.06 4031.57 4457.53 10.44 4816.51 4824.00 17.35 4816.74 4818.88 20.99 urbain_21 3342.50 3358.59 13.70 3256.36 3260.79 16.41 3211.39 3211.95 16.27 3210.53 3211.52 15.09 3210.53 3211.53 14.29 3210.53 3211.35 13.36 urbain_22 3356.74 3361.97 17.35 3356.74 3361.97 11.83 3356.74 3361.97 17.19 3356.74 3361.97 18.47 3214.70 3853.71 18.44 3214.70 3853.35 15.23

	% of EVs in the fleet			0%			20%				40%			60%				80%	100%
	Instance	Best		Avg.	t(min) Best	Avg.	t(min) Best	Avg.	t(min) Best	Avg.	t(min) Best	Avg.	t(min) Best	Avg.	t(min)
	rural_18	7309.13 7347.27	6.13 7062.45 7104.56 10.36 6856.54 6905.22 10.79 6791.60 6817.33 13.60 6639.04 6685.71	N.R.	N.R.	N.R.	N.R.
	rural_19	4568.71 4577.93 10.07 4337.86 4354.07 17.17 4243.67 4249.34	8.72 4243.14 4249.58 10.79 4243.14 4247.26	N.R.	N.R.	N.R.	N.R.
	rural_20	6533.76 6552.38	9.74 6263.10 6277.22 19.66 6066.01 6079.87 21.25 5976.84 6013.09 18.53 5885.27 5934.11	N.R.	N.R.	N.R.	N.R.
	rural_21	8465.34 8508.08 18.73	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.	N.R.
	rural_22	6611.81 6700.23 20.53 6356.17 6475.38 20.17 6131.63 6181.34 20.36 6033.61 6068.19 19.41 5863.80 5871.67	N.R.	N.R.	N.R.	N.R.
	semi_urbain_18 4362.85 4388Avg. Gap 1.76							4.73			3.03			4.18	3.90
	Avg. Time				16.80			18.12			17.85				17.34		18.54	18.57
	Avg. Gap above BKS			15.97%			1.38%			0.17%	0.92%				4.50%	4.97%
	NBKS				0			0				4					0	0
	Cum. Number of Veh.		508			461									466		459
	Avg. Time (min)							159.58				6.20				157.03	156.05
				MSH(1k)				MSH(5k)				MSH(10k)
	Instance	n	v	Best	Avg.	t	n	v	Best		Avg.	t	n	v	Best	Avg.	t
	111c_21s	109	17 4798.52 4861.35 0.91 109	17	4780.01 4793.07	2.40	109	17 4777.91 4781.85	4.94
	111c_22s	109	17 4795.26 4858.09 1.04 109	17	4776.75 4789.81	2.13	109	17 4774.65 4778.80	4.69
	111c_24s	109	17 4795.19 4855.89 1.16 109	17	4776.68 4789.77	2.74	109	17 4773.67 4778.62	5.64
	111c_26s	109	17 4795.19 4855.89 0.79 109	17	4776.68 4789.77	2.56	109	17 4773.67 4778.62	5.23
	111c_28s	109	17 4793.57 4854.86 0.93 109	17	4775.06 4788.22	2.80	109	17 4772.46 4777.03	5.54
	200c_21s	192	32 9005.58 9039.97 2.96 192	31	8894.56 8923.18 19.80	192	31 8839.62 8879.98 19.96
	250c_21s	237	39 10702.76 10755.30 13.43 237	38	10534.52 10579.12 27.59	237	37 10482.52 10518.32 21.58
	300c_21s	283	45 12663.49 12719.25 44.28 283	44	12444.48 12548.70 37.43	283	44 12367.60 12421.75 47.53
	350c_21s	329	52 14431.27 14470.40 53.01 329	50	14146.67 14253.55 50.96	329	50 14073.34 14226.03 63.01
	400c_21s	378	60 16873.51 17143.75 64.29 378	59	16745.24 17128.76 67.91	378	59 16660.20 17119.89 71.70
	450c_21s	424	66 18569.58 18831.00 72.11 424	65	18351.72 19009.22 76.35	424	65 18241.48 18902.03 80.75
	500c_21s	471	75 20960.18 21678.33 80.11 471	73	20610.36 21297.07 84.95	471	73 20496.50 20997.04 89.95
	Avg. Gap above BKS		1.42%	2.64%				0.40%	1.48%			0.05%	1.02%
	NBKS			0					0						8
	Cum. Number of Veh.						445						444
	Avg. Time (min)				27.92					31.47					35.04

Table A .

 A 8 -Results of PMa on large instances of Hiermann et al. (2016).

	Type		A					B					C
		ALNS		PMa		ALNS		PMa		ALNS		PMa
	Instance	BKS				BKS					BKS
		Avg.	t (min)	Best	Avg.	t (min)	Avg.	t (min)	Best	Avg.	t (min)	Avg.	t (min)	Best	Avg.	t (min)
	c101	7162.01 7190.21	17.53 7162.01 7165.15	25.88 2495.00 2505.73	14.43 2495.00 2502.09	25.56 1809.93 1816.06	14.08 1810.12 1820.45	8.84
	c102	7139.45 7162.24	17.95 7139.45 7142.80	26.02 2445.99 2450.73	14.41 2448.76 2451.39	25.86 1759.73 1766.14	14.36 1763.97 1770.44	9.20
	c103	7121.12 7149.31	18.30 7121.12 7124.26	31.22 2430.72 2452.40	15.03 2430.72 2437.99	31.16 1746.80 1759.20	15.09 1746.80 1757.68	20.76
	c104	7099.88 7110.43	17.75 7099.88 7103.79	32.93 2404.73 2428.95	15.07 2404.73 2410.12	32.85 1719.67 1735.86	15.60 1721.53 1727.72	29.04
	c105	7140.31 7182.56	17.60 7140.31 7152.21	27.24 2472.93 2475.95	14.45 2473.79 2476.46	27.16 1783.25 1785.43	14.51 1785.31 1790.59	10.66
	c106	7136.90 7168.94	17.69 7136.90 7141.66	27.28 2462.54 2468.13	14.71 2462.90 2466.25	27.20 1774.77 1777.67	14.64 1775.08 1780.39	11.10
	c107	7139.12 7171.04	17.42 7139.12 7145.54	30.49 2458.37 2461.32	14.60 2458.37 2459.31	30.25 1764.02 1768.33	14.57 1765.45 1775.74	14.68
	c108	7133.30 7153.77	17.53 7133.30 7142.39	30.22 2450.17 2463.02	14.33 2452.58 2455.74	30.17 1761.41 1769.76	14.80 1768.05 1772.78	20.58
	c109	7120.33 7132.19	17.32 7123.15 7132.66	33.32 2431.87 2452.57	15.07 2431.87 2434.02	33.22 1740.18 1749.07	15.30 1744.09 1749.05	29.92
	c201	5736.35 5757.53	39.28 5736.35 5741.28	19.36 1730.41 1739.26	35.76 1730.41 1730.50	14.14 1210.41 1213.63	31.30 1210.41 1210.50	13.50
	c202	5740.43 5765.52	42.22 5740.43 5753.56	21.50 1729.73 1745.24	38.14 1729.73 1730.12	17.01 1209.73 1220.97	34.17 1209.73 1210.12	16.58
	c203	5726.08 5751.99	47.69 5744.25 5783.33	28.27 1716.29 1742.76	39.38 1730.26 1743.23	23.73 1212.34 1227.69	33.07 1216.18 1226.22	30.01
	c204	5705.82 5727.18	50.91 5724.75 5772.72	41.49 1699.07 1709.43	37.02 1724.75 1753.55	36.43 1179.25 1199.37	32.15 1204.68 1232.00	35.62
	c205	5694.58 5725.41	46.16 5694.58 5701.83	24.43 1694.58 1715.09	37.83 1694.58 1701.83	17.33 1188.92 1195.24	31.64 1192.54 1195.35	16.70
	c206	5689.08 5714.39	31.20 5689.08 5696.99	26.40 1689.08 1712.38	36.78 1689.08 1696.99	20.69 1183.42 1192.30	31.01 1188.81 1191.74	20.16
	c207	5696.54 5713.44	32.58 5696.54 5716.37	27.61 1694.61 1710.70	35.05 1696.54 1715.61	22.19 1183.42 1190.37	31.44 1187.49 1190.65	21.52
	c208	5682.60 5707.65	50.23 5682.60 5700.58	30.48 1681.47 1707.11	36.46 1682.60 1701.29	24.49 1181.47 1192.96	30.82 1182.60 1192.42	23.78
	r101	4366.21 4465.51	14.47 4366.21 4379.64	21.41 2249.14 2281.28	13.68 2253.04 2253.93	10.71 1954.00 1977.89	14.36 1959.04 1959.07	9.45
	r102	4176.82 4270.92	14.97 4176.82 4180.61	13.11 2047.89 2095.87	14.45 2053.42 2057.97	8.01 1757.13 1791.03	14.67 1761.48 1767.00	7.54
	r103	4043.29 4130.86	16.43 4043.29 4049.01	15.53 1892.84 1927.52	15.11 1892.84 1897.06	13.13 1587.32 1618.81	15.52 1587.32 1593.13	12.54
	r104	3969.12 4025.60	15.27 3969.12 3970.77	21.45 1747.65 1775.33	15.28 1749.48 1761.84	22.78 1424.30 1448.31	16.20 1427.76 1437.37	19.66
	r105	4138.01 4215.34	15.37 4138.01 4140.08	12.94 1997.75 2030.12	15.08 1998.80 2007.01	9.17 1699.34 1728.12	14.68 1699.34 1705.56	8.25
	r106	4068.15 4155.24	15.54 4068.15 4075.44	16.46 1925.83 1963.88	14.51 1925.83 1929.13	12.02 1604.55 1635.42	14.67 1609.78 1612.20	8.55
	r107	4003.84 4093.59	15.30 4003.84 4008.71	24.07 1818.03 1844.20	15.50 1818.03 1820.75	22.10 1490.04 1514.01	15.99 1493.12 1499.00	17.35
	r108	3955.48 4025.75	15.77 3955.48 3961.93	23.96 1720.78 1753.09	16.27 1720.78 1733.87	23.31 1397.86 1417.39	16.42 1397.86 1404.05	18.13
	r109	4016.33 4110.98	15.58 4016.33 4024.49	19.92 1863.48 1904.12	15.37 1863.48 1867.32	25.36 1550.40 1580.14	15.65 1551.69 1554.22	16.98
	r110	3969.02 4045.96	15.73 3969.02 3975.95	26.51 1752.40 1793.48	15.54 1752.40 1755.92	12.89 1420.13 1471.66	15.64 1420.13 1429.38	12.11
	r111	3972.84 4048.42	15.93 3972.84 3976.91	25.10 1765.08 1808.36	15.66 1765.08 1775.04	22.45 1438.81 1479.75	16.16 1439.02 1443.94	12.62
	r112	3942.68 4023.01	15.80 3942.68 3947.93	21.85 1714.93 1746.02	15.77 1714.93 1722.34	26.75 1388.22 1403.82	16.10 1388.22 1394.19	25.74
	r201	3413.93 3432.83	42.20 3515.45 3543.81	17.91 1594.58 1618.25	31.21 1697.28 1714.30	40.11 1366.63 1378.77	29.69 1474.75 1480.55	39.50
	r202	3270.49 3295.26	44.95 3366.13 3397.33	18.23 1468.05 1479.42	29.61 1541.45 1551.73	17.87 1236.97 1249.65	29.13 1316.45 1328.35	42.06
	r203	3136.47 3169.97	49.40 3224.34 3252.78	23.86 1340.00 1354.24	30.70 1424.34 1438.89	23.68 1104.85 1124.07	30.23 1183.71 1205.32	56.20
	r204	3008.01 3026.09	46.32 3098.94 3114.40	26.84 1203.89 1211.63	27.35 1278.59 1286.75	27.33 977.72 983.97	26.81 1032.48 1050.51	65.96
	r205	3234.26 3261.16	40.89 3326.29 3356.01	21.28 1430.70 1455.08	30.16 1524.53 1534.24	20.82 1217.77 1232.63	29.15 1293.54 1308.90	47.77
	r206	3172.50 3194.12	47.73 3261.05 3296.06	26.25 1361.69 1376.34	31.35 1461.05 1470.36	27.07 1136.83 1155.47	30.95 1230.67 1239.18	60.69
	r207	3079.39 3099.52	46.87 3159.76 3175.13	23.25 1256.22 1268.66	28.18 1359.76 1364.12	23.60 1031.22 1057.22	26.31 1110.93 1136.77	58.58
	r208	3010.51 3026.57	51.26 3088.29 3104.77	33.89 1198.39 1208.89	29.02 1281.92 1289.03	34.43 971.15 984.87	28.21 1032.20 1049.20	85.72
	r209	3142.72 3161.57	45.06 3224.28 3255.62	25.29 1333.33 1345.50	30.30 1381.34 1407.39	25.18 1099.24 1117.68	29.62 1156.34 1170.58	57.01
	r210	3110.90 3143.79	45.94 3181.95 3235.09	25.02 1314.16 1324.07	30.07 1359.10 1385.87	25.44 1087.21 1100.27	29.88 1134.10 1143.99	57.33
	r211	3041.93 3079.24	44.29 3133.26 3147.32	20.83 1231.38 1244.73	26.30 1319.15 1326.35	21.33 1006.38 1026.07	25.93 1087.95 1098.37	51.81
	rc101	5247.39 5346.49	14.13 5247.39 5255.03	13.21 2504.72 2560.33	13.71 2504.72 2508.65	9.42 2121.02 2153.24	13.80 2121.02 2123.26	8.96
	rc102	5114.15 5180.03	14.63 5114.15 5118.29	24.81 2330.50 2359.92	14.35 2334.08 2335.72	9.68 1947.54 1972.85	14.68 1947.54 1950.16	9.40
	rc103	4916.98 5007.37	14.53 4916.98 4924.38	18.60 2105.84 2136.78	14.14 2108.82 2113.31	10.78 1726.85 1764.22	14.54 1726.85 1736.49	10.72
	rc104	4801.06 4862.65	16.03 4801.06 4815.22	24.28 1979.16 2002.33	15.56 1979.16 1984.66	14.39 1595.44 1614.09	15.80 1596.18 1605.91	15.20
	rc105	5060.96 5117.09	14.48 5060.96 5068.85	14.74 2252.05 2287.95	13.87 2252.05 2264.49	9.81 1875.91 1900.42	14.14 1875.91 1880.03	9.93
														Continued on next page

Thèse de Doctorat Alejandro MONTOYA Problèmes de tournées de véhicules électriques : modèles et méthodes de résolution Electric Vehicle Routing Problems: models and solution approaches

These costs include energy (fuel or electricity), maintenance, and taxes

Although in the green vehicle routing problem literature the most popular acronym is G-VRP, we prefer to use Green VRP

This graph resembles that used by[START_REF] Villegas | Grasp/vnd and multi-start evolutionary local search for the single truck and trailer routing problem with satellite depots[END_REF] in a route-first cluster-second heuristic for the single truck and trailer routing problem with satellite depots.

This metric is the sum of the number of vehicles of all best solutions[START_REF] Bräysy | Vehicle routing problem with time windows, part II: Metaheuristics[END_REF]

This comparison tool was introduced by Vidal et al. (2015)

G * : Acyclic graph for the split procedure G * = (V * , A). V * : Set of vertices of G * . A: Set of arcs of G * . K: Set of TSP-like tours. R: Set of all feasible routes. Ω: Subset of R. S: Set of all feasible solutions to Green VRP. H: Set of sampling heuristics. K: Number of iterations in sampling phase. tsp: TSP tour (tsp ∈ K). s: Green VRP solution (s ∈ S). W: Set of customers visited by tsp. N : Ordered set of nonrouted customers. L: Randomization factor of each heuristic. f (s): Objective function of solution s. s * : Best solution found during execution of method. λ: Positive parameter. c r : Total distance of route r. tp r : Total time of route r. f r : Binary variable equal to 1 if route is feasible. label i : Shortest path label i. pred i : Predecessor labels i.

In our implementation we use intra-route and inter-route versions with best-improvement selection.

We randomly picked the route from the best solution found to a randomly picked instance

The best gap is the gap between the best solution found over 10 runs and the optimal solution

According to(Lahyani et al.

2015), if a VRP has at least nine physical characteristics, this problem can be classified as a RVRP

In our implementation we use intra-route and inter-route versions with best-improvement selection.

The best gap is the gap between the best solution found over 10 runs and the optimal solution

For confidentiality reasons, we do not report detailed results for each instance in this testbed

Hiermann et al. (2016) report the BKSs obtained using BnP or any run of their experiments

Notice that, because of Hiermann et al. (2016) just report the average solution, we could not compute the number of BKSs and the average best gap for the ALNS

This values were calculated by multiplying the distance travelled by the CVs by a constant emission rate available in their technical information

The increment on the cost is explained by two factors: i) the solutions employed more vehicles, and ii) CVs have to travel more distance (which leads to higher travel costs). Figure5.4a shows the comparison of the maximum number of routes over the 10 runs with and without the option of visit the CSs. Although in most of cases the number of vehicles is the same, for instances rural_18 and rural_22 the solutions employ an additional route when EVs are 80% of the fleet. This explains the large gap shown in Figure5.3 for these instances. Figure5.4b shows the gap in the average travelled distance of the CVs on the instances without mid-route charge with respect to the original solutions reported in Section 5.6.4. Since the EVs can not use the CSs, their autonomy is limited to the battery level when they leave the depot. Then, the CVs have to visit more request, and thus they have to travel more distance.8 Note that re-running the algorithm on instances with no EVs in the fleet would lead to the same solutions and that running the algorithm on instances with only EVs on the fleet would lead again to unfeasible solutions 9 G(%) = (cost withoutCSs -cost withCSs )/cost withCSs × 100

Hiermann et al. (2016) report the BKSs obtained using BnP or any run of their experiments

List of Figures

Cumulative Frequency

Battery level after charging (%)

e ij : Energy consumption between vertices i and j (i, j ∈ V ). t ij : Travel time between vertices i and j (i, j ∈ V ).

T max : Tour duration constraint. Q: Tank capacity.

Notation for hybrid metaheuristic

G * : Acyclic graph for the split procedure G * = (V * , A * ).

V * : Set of vertices of G * . A * : Set of arcs in G * . Ω: Pool of feasible eVRP-NL routes.

4.6.3 Notation for the fixed-route vehicle-charging problem Π: Energy-infeasible route. t: Total time of the route Π, which is the sum of the travel times plus the service times. ρ jk : Slope of the segment defined between the breakpoints k -1 and k ∈ B. ε π(i)j : Binary variable, equal to 1 if the EV charges at CS j ∈ F before visiting vertex π(i) ∈ Π. φ π(i) : Battery level. If ε π(i)j = 0, φ π(i) is the battery level when the EV arrives at vertex π(i). On the other hand, if ε π(i)j = 1, φ π(i) is the battery level when the EV arrives at CS j ∈ F right before visiting vertex π(i). θ π(i)jk : Binary variable, equal to 1 if the EV charges on the segment defined by breakpoints k -1 and k ∈ B at CS j ∈ F before visiting vertex π(i) ∈ Π. δ π(i)jk : Amount of energy charged when the charging finishes on the segment between breakpoints k -1 and k ∈ B at CS j ∈ F before the visit to vertex π(i) ∈ Π. µ π(i)jk : Battery level when the charging finishes on the segment between breakpoints k -1 and k ∈ B at CS j ∈ F before the visit to vertex π(i) ∈ Π. τ : Minimum time needed to visit CS j ∈ F between two vertices i and h ∈ I ∪ {0}. u: Lower bound on the travel time. v: Lower bound on the charging time. ec: Charge needed to recover the energy consumed to make the detour to j ∈ F . o: Lower bound on the energy remaining when an EV arrives at customer i ∈ I right after charging at CS j ∈ F . c(i): Closest station to the customer i ∈ I (including the depot). t : Lower bound on the new duration of a given fixed route if CS j ∈ F is inserted between vertices π(i) and π(i + 1) ∈ Π, i = n r .

Y i : Battery level at each vertex i ∈ Π, assuming that the EV fully charges its battery at each visited CS. s: Sum of the battery levels with negative values. t: Total time of the route (supposing a full charging policy).

Chapter 5

The technician routing problem with conventional and electric vehicles

The work of this chapter was presented at the following conferences and seminar:

-Montoya, A., Guéret, C., Mendoza solutions, we conducted a small study. We ran our PMa on the rural instances with 20% to 80% 8 of EVs in the fleet, considering that the EVs departs from the depot with the battery full. We compared these results to those obtained in the experiments reported in Section 5.6.4. We were able to obtain feasible solutions for all the tested instances. However, our results suggest that disallowing mid-route charging has a impact on the quality of the solution. Figure 5.3 shows the average gap between the cost of the solutions with and without the option of visiting CSs 9 . between technicians and requests; incompatibility between CSs and EVS; nonlinear battery charging functions; and a limited number of vehicles for each type. Decisions that must be taken simultaneously are the vehicle-to-technician assignment, the sequencing of the routes, and the battery charging program.

To solve the TRP-CEV we proposed a parallel matheuristic (PMa). This approach decomposes the TRP-CEV into a set of VRP-TW and eVRP-TW with lunch break (one for each couple technician-vehicle), and solves each sub-problem in parallel using a greedy randomized adaptive search procedure (GRASP). Later, PMa uses a set covering model to assemble a TRP-CEV solution. We built a set of real instances based on the ENEDIS operation. The aim of these instances is to evaluate the performance of the PMa, and to analyse the features of the solutions under the light of some important metrics for ENEDIS. The proposed PMa provides competitive results when compared against the current routing software used by the company and state-of-the-art methods for a related problem (the EFSMTW). Finally, we analyzed the solutions delivered by PMa for the TRP-CEV aiming to provide ENEDIS team with some insight about the behaviour of costs, emissions, and visits to CSs for different compositions of the fleet. Our analysis concluded that the cost decreases in function of the number of EVs in the fleet, and that the visits to the CSs are relevant only for the rural instances. t ij : Travel time between vertices i and j (i, j ∈ N ). e ijk : Energy consumption between vertices i and j (i, j ∈ V ) of the type of vehicle k ∈ K e .

Notation for the MILP formulation

F : Set containing the set F and the copies of CSs. N : Set of vertex containing the set F . V : Set of vehicles. V e : Set of electric vehicles. V c : Set of conventional vehicles. n il : Binary parameter, equal to 1 if request i ∈ C requires a technician with skill l ∈ S, and 0 otherwise. m ul : Binary parameter, equal to 1 if technician u ∈ T has skill l ∈ S, and 0 otherwise. h ik : Binary parameter, equal to 1 if the EV k ∈ V e can charge at CS i ∈ F . S max : Maximum charging time of any charging function.

x ijku : Binary variable, equal to 1 if technician u ∈ T travels from vertex i to j in vehicle k ∈ V , and 0 otherwise. τ iu : Arrival time of technician u ∈ T at vertex i ∈ N . y ik : Battery level of EV k ∈ V e upon departure from vertex i ∈ N . q ik : Battery level when EV k ∈ V e arrives at CS i ∈ F . o ik : Battery level when EV k ∈ V e departs from CS i ∈ F . α ikl : Binary variable, equal to 1 if the battery level is l ∈ B when vehicle k ∈ V e arrives at CS i ∈ F . β ikl : Binary variable, equal to 1 if the battery level is l ∈ B when vehicle k ∈ V e departs from CS i ∈ F . ∆ ik : Time spent by vehicle k ∈ V e at CS i ∈ F . w iju : Binary variable, equal to 1 if technician u ∈ T takes the lunch break between the vertex i ∈ N and j ∈ N and 0 otherwise. z iju : Binary variable, equal to 1 if the technician u ∈ T takes the lunch break after travelling from vertex i to j ∈ N and 0 otherwise. z iju : Binary variable, equal to 1 if the technician u ∈ T takes the lunch break before travelling from vertex i to j ∈ N , and 0 otherwise. Erdogan & Miller-Hooks (2012) reported the best distance (Best), and the number of vehicles of the best found solution (v) of multiple runs with different parameters for their MCWS and DBCA; they did not give the exact number of runs. [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF] reported the best distance, the number of vehicles of the best found solution, and the average computing time (t in minutes) over ten runs of their VNS/TS, [START_REF] Schneider | An adaptive VNS algorithm for vehicle routing problems with intermediate stops[END_REF] reported the best distance, the average distance (Avg.), and the average computing time over ten runs of their AVNS. Finally, [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF] reported the best distance, the number of vehicles of the best found solution, and the computing time over a single run of their 48A and SA. For all the algorithms, we provide the number of customers served (n). The last rows of the table summarize the average BKS gap, the cumulative number of vehicles, the number of times each method found the BKS, and the average running time. Values in bold indicate that a method found the BKS. A.2 Detailed results of the eVRP-NL Tables A. 3 and A.4 show the results of our two ILS versions (i.e, ILS(S)+HC and MILP(H)+HC) for the small and large eVRP-NL instances. In Table A.3, we compare our results with the optimal solutions found by Gurobi using the MILP formulation. In Table A.4, we compare our results with the best results obtained with Gurobi. For each instance, we report the problem name 1 , and the best solution (BKS) taken from the results of Gurobi, ILS(H)+HC or ILS(S)+HC.

Notation for the PMa

For the results obtained with Gurobi, we report the best solution (Best), and the gap with respect to the BKS (G) 2 . For the results obtained with the ILS(H)+HC and ILS(S)+HC, we report the best solution, the average solution (Avg.), and the average computing time (t in seconds) over ten runs. For the ILS(S)+HC, we also report the average computing time using the preprocessing procedure (t * ). For the two ILS+HC versions, we provide the gap of the average solution and best solution with reference to the BKS. The last rows of the table summarize the average and maximum BKS gap, the number of times each method found the BKS, and the average and maximum running time. Values in bold indicate that a method found the BKS. 1 tcα α αcβ β βsµ µ µc # # #, where α α α is the type of the location of the customers (i.e., 0:randomize, 1:cluster, 2: mixture of both), β β β is the number of customers, µ µ µ is the number of the CSs, is 't' if we use a p-median heuristic to locate the CSs and 'f' otherwise, and # # # is the number of the instance for each combination of parameters (i.e., # = 0, 1, 2, 3, 4)

2 G(%) = (of -of BKS )/of BKS × 100 A.3 Detailed results for the TRP-CEV Tables A. 5 shows the results of our PMa for the small TRP-CEV instances. We compare our results with the optimal solutions found by Gurobi using the MILP formulation presented in Section 3.3.2. For the results obtained with the PMa, we report the best solution, the average solution (Avg.), and the average computing time (t in seconds) over ten runs. The last rows of the table summarize the average gap, the number of times PMa found the BKS, and the average running time. Values in bold indicate that a method found the BKS. Tables A .6 shows the results of our PMa for the large TRP-CEV instances. For each instance and percentage of EVs in the fleet, we report the best solution, the average solution (Avg.), and the average computing time (t in seconds) over ten runs. The last rows of the table summarize the average gap, and the average running time. 3 and updated with some new BKSs found by our PMa.

Hiermann et al. ( 2016) reported the best solution (Best) only for the small instances, the average cost (Avg.), and the average computing time (t in minutes) over ten runs of their AVNS. The last rows of the table summarize the average BKS gap, the number of times each method found the BKS, and the average running time. Values in bold indicate that a method found the BKS.