Ciro Amaya

Cyril Briand

Véronique

Pascale,S Isabelle Ingrid Bernadette

Majo, Camilo, Sebas Juan Jose

Juan Pablo

Jorge Caro

Guille, Adri, Elyn, Syrine Marie Karen

Andrea Felipe

Aleja, Samir, So, Nacef Mourad David

It

Table 3.4: VRPSD Proposed testbed results -Costs: Continued

First of all I will like to thank You. After all you are the purpose and reason of my life. Then, I would like to thank my advisers Nacima Labadie and Nubia Velasco

Chapter 1 Introduction

The world of science lives fairly comfortably with paradox. We know that light is a wave, and also that light is a particle. The discoveries made in the innitely small world of particle physics indicate randomness and chance, and I do not nd it any more dicult to live with the paradox of a universe of randomness and chance and a universe of pattern and purpose than I do with light as a wave and light as a particle. Living with contradiction is nothing new to the human being" Madeleine L'Engle

The transport activities are one of the more important drivers in many of the economic activities.

In fact, according with the International Trade Administrator, logistics and transportation activities 1 accounted for 8% of the Gross Domestic Product in the United States in 2015. Additionally, only considering the logistic activities, transportation can represent up to 60% of the total costs [16], making this topic an important area to study.

In the Operations Research context, the transportation problems and specically the Vehicle Routing Problems (VRPs) has been one of the most studied problems with more than a thousand published papers between 1954 to 2006 [6]. Overall, the VRP consist of nding a set of minimal cost routes performed by a set of vehicles, satisfying a set of clients, and respecting specic constraints, according to the particular context. Furthermore, most of the works on this eld consider that all problem parameters, travel and service times, and demands, are known in advance. Actually, Braekers et al. [3] identify that more than 80% of 277 published articles between 2009 and 2015 in the VRPs context are deterministic.

Many factors could aect the certitude of the information on problems parameters. For example, the cities population concentration and their consequent trac congestion make the travel times uncertain [11]. Furthermore, the accelerated use of Information Technologies could generate incertitude about the requests because clients could made them more frequently and are subject to random external factors. Moreover, the time spent at each location visited by the vehicles can change because of the complexity of the task to satisfy at the clients or due to the environmental conditions.

The consequence of neglecting the variability is that the solutions obtained could perform badly in real uncertain environment [17]. Actually, based on a previously literature results, ignoring the incertitude of the information could increase up to 10% the costs for the stochastic demands case [8] or 4% for stochastic times [1]. Moreover, the amount of the objective (costs or utility) variation between deterministic and stochastic solutions depends on the randomness of the problem and can be as signicant as 20% when travel and service are uncertain [4].

In recent years, some authors have been working on stochastic VRPs, one of the most comprehensive study about this topic is given by Gendreau et al. [10]. The authors focus their attention on the stochastic programming modeling which is the predominantly approach in the eld. Thus, only parameters which are dened by random variables or scenarios are considered 2 . Furthermore, to classify the dierent problems the authors consider two main characteristics: the solution paradigm and the stochastic parameters (customers, demands, and times). Solution paradigms are divided in two: a priori and reoptimization paradigms. The former, consider that solutions are created before any information is revealed (rst phase) and they are barely modied during the execution (second phase). Meanwhile, reoptimization approach aims to modify the solution as new information becomes available to improve. The latter approach is more related with the Dynamic problems (see Psaraftis et al. [18] for a recent review). A priori paradigms can be further divided into Stochastic Problems with Recourse (SPR) and Chance Constrained Problems (CCP). The SPR use recourse to react to failures, or constraints violations during the second phase. Meanwhile, the CCP bounds the probability of possible failures appearing during the second phase.

The a priori paradigm has been the predominant approach to solve stochastic VRPs, particularly using SPRs. This can be explained by the fact that this type of models allow stable tactical routes which are operationally desirable [9]. Moreover, according Bekta et al. [2], a priori paradigm is suitable when anticipating uncertainty is crucial to nd feasible solutions and avoid penalties (economic and reputation). Nevertheless, the a priori paradigm adds complexity of dealing with probability calculus overhead. This makes that the size of stochastic problems that can be solved (exactly and approximately) is rather small. For example, Gauvin et al. [8] are able to optimally solve only one instance with 100 customers in the context of the VRP with stochastic demands. Therefore, there is a need to develop solution methods able to tackle closer to real problems settings, in terms of size, multiple uncertainties, and assumptions.

The literature review reveals that there is a lack of comparisons among dierent methods and works between stochastic VRPs, which might be caused by the lack of standardized benchmarks.

Although each problem has its own characteristics (type of distribution, amount of variance, etc.) base tests serve to prove the usefulness of new models and strategies to solve the stochastic VRPs.

Moreover, these test cases need to evolve, particularly in terms of size, to prove the capacity of the new methods to deal with closer to life real problems. Being able to solve this type of problems will demonstrate their capacity to fully exploit the benets of stochastic solutions over deterministic ones.

Last but not least, new stochastic models need to incorporate real life constraints, such as hard time windows, a characteristic that has been largely studied in deterministic VRPs but not so much in the stochastic ones. This thesis addresses two kinds of VRPs using the stochastic programming framework and the a priori paradigm: the rst one considers the demand as a random variable and is reported in chapter 3 where results for middle and large instances are reported, second it tackles in chapter 4 the stochasticity on the travel and service time. The latter under the presence of hard time windows which can conduce to unserviced customers and using dierent types of continuous distribution for the stochastic parameters. Both problems are presented in the context of maintenance operations for which anticipating uncertainty is imperative. The last problem studied in this thesis considers maintenance planning on wind farms. It extends the stochastic VRPs in which technicians are to be scheduled to perform their tasks under the appearance of uncertain new tasks (demands), random weather conditions, and with stochastic service times. The present thesis is developed as follow:

Chapter 2 introduces an extensive review on the Vehicle Routing Problems (VRPs). Starting with the description of deterministic VRPs, it makes its path to uncertain VRPs as their natural evolution. Then, the attention is focused in the three main paradigms to model uncertain VRPs, namely, Stochastic Optimization, Interval Optimization and Fuzzy Logic. Special consideration is given to the Static Stochastic VRPs with a comprehensive review of the solution approaches and dierent problems variants tackled in the literature. This revision shows that albeit the increase of CHAPTER 1. INTRODUCTION research in the eld, there is still a lack of detailed results for big instances. Moreover, a lack of research for the stochastic VRPs with hard time windows is established. Especially for the case with stochastic travel and service times modelled by continuous random variables, so it is likely they are in real applications.

Chapter 3 is devoted to the VRP with stochastic demands (VRPSD). In the VRPSD, customers demands are modeled by random variables and their realization value are only known when the vehicles arrive at the customers. A simple classical recourse action is used to model the VRPSD as a stochastic problem with resource. To tackle the VRPSD a Greedy Randomized Adaptive Search Procedure (GRASP) is used to restart a Memetic Algorithm (MA) and eciently solve the problem at hand. In this chapter it is shown that large instances (up to 385 customers) can be eectively handled by the MA+GRASP. A comparison with the state-of-the-art algorithms for the VRPSD shows that the MA+GRASP provides better and more accurate solutions in very competitive computational times. Moreover, the chapter establishes a new testbed of instances (based on instances already used in deterministic context) with a higher number of customers than the traditional Christiansen and Lysgard benchmark [5]. These results are important to open the space to discussion and further design of other methods to VRPSD. The work presented in this chapter is under minor revision in the Computers & Operations Research journal and an earlier version was presented at the CIE45 conference [12].

Chapter 4 focuses on a VRP with uncertain times. It presents a VRP considering stochastic travel and service times with hard time windows. The problem is thought-out in a maintenance activities context and the uncertainty in times are modeled through continuous probability distributions. A model is proposed to enable the control of customers service levels but also considers the implications of missing the time windows. To overcome the problem of modeling the arrival times, it is shown that they can be fairly approximated using a log normal distribution. To solve the problem, a Multi-population Memetic Algorithm (MPMA) exploiting dierent characteristics in each population (running in parallel) is proposed. Results are presented for instances with up to 100 customers derived from the the Solomon [19] benchmark. Additionally, the MPMA is compared against state of the art methods although these allow late services, and the proposed approach shows very good performance.

The results of this third chapter are gathered in a paper which is accepted for publication in the Computers & Industrial Engineering journal. Preliminary results were presented at MIM2016 [14] and CLAIO2016 [13] conferences.

Chapter 5 introduces a general review of the wind farms maintenance activities. Two related problems are further explored in this context. First, a multi-objective approach to deal with maintenance scheduling of wind farms is addressed. In this problem the operator of the wind farm needs to decide the order, time, and resources assignation to execute a set of maintenance tasks in a short term horizon. Moreover, the operator aims to minimize its costs while the investors want to maximize the energy production. A linear integer model is used to model the problem and the epsilon-constraint method is designed to approximate the optimal Pareto front. Tests are performed on the set of instances proposed by Froger et al. [7] pointing out that objectives are in conict. Also it is shown that the variation of energy production in the short term can be highly aected by the scheduling of the activities. The second problem, extends the scheduling problem and explores the selection of maintenance strategies for wind farms. Dierent strategies are evaluated within a event discrete simulation approach to compare them on a long-term horizon. Furthermore, dierent ways of solving the scheduling the maintenance tasks in the short term are compared within the simulation. Failures appearances as well as maintenance times are considered as random variables. The proposed shows that even simple heuristic rules are used to tackle the scheduling of technicians, they can have important eects on both the costs and the energy production. The results for the scheduling maintenance activities problem considering multiple objectives were presented at IEOM 2017 conference held at Bogota [15].

Finally the thesis ends with chapter 6 drawing the conclusions of the thesis and research clues for future research on the uncertain vehicle routing problems, scheduling of maintenance activities, and strategy selection in the wind farm context.

Chapter 2

Literature Review -Vehicle Routing Problems Nearly 60 years have passed since the introduction of the vehicle routing problem (VRP). The rst literature appearance of the VRP can be traced back to the seminal work of Dantzig and Ramser [41] named The Truck Dispatching Problem. Since 1959, the number of articles and applications have grown tremendously. Since that year, a general search for the Vehicle Routing Problem in a searching engine such as Google Scholar gives more than 20.000 results that can be reduced to over 4600 if the words are present in the document title. Eksioglu et al. [49] reviewed nearly 1500 VRP references from 1954 to 2006 1 , which included journal articles, books, book chapters, technical reports, and conference articles. The authors proposed a taxonomy to classify the vast literature and asserted that it grew exponentially at a rate of approximately six percent per year. A more recent classication work can be found in Braekers et al. [30] where 277 VRP journal articles from 2009 to mid-2015 were arranged using a taxonomy similar to the one used in Eksioglu et al. [49].

The massive amount of research related to the VRP can be twofold explained. First, transportation plays a central part in many human activities, economics, and the environment. According to Hesse and Rodrigue [82] transportation accounted for nearly 6% of the Gross Domestic Product (GDP) of the United States in the year 2000. Moreover, transportation transcends the purely economic trend. In fact, the subject has been studied in the context of disaster relief and humanitarian logistics ( [74,34]), and services delivery (health care [55], technicians [START_REF] Ribeiro | A simple and robust simulated annealing algorithm for scheduling workover rigs on onshore oil elds[END_REF], among other). Second, the eld has been the seed for many developments of several exact and heuristic methods for combinatorial optimization problems [START_REF] Laporte | Fifty years of vehicle routing[END_REF]. These developments have an impact in other elds in the Operational Research community [84] and therefore make the VRP research an active and important part of the scientic development.

Besides, the nearly sexagenarian problem has seen a myriad of variants and extensions of its basic version. Either by the addition of more characteristics, constraints, or changes on the objective function, new problems have risen to adapt the VRP to many contexts. Within this variety, the last few decades have seen an increment in the study of problems where the parameters information is not certain [63]. Beyond the pure theoretical value of the works, applications deal with a reality in which information is far from perfect and stochastic (weather, accidents, drivers skills, etc.). Moreover, despite the usual necessary eort to solve problems with uncertain information, its value is not trivial [6,33,60]. Therefore, the uncertain VRPs are an important eld to make both theoretical and practical research. This chapter presents a review of the VRP. It starts by introducing one of the most basic and known version of the problem, the Capacitated Vehicle Routing Problem (CVRP). The CVRP serves 2.1. DETERMINISTIC VRP as proxy to introduce the mathematical formulations and a very special generalization called the VRP with Time Windows (VRPTW). Then, the uncertain VRPs modeling and solution approaches are explored. The chapter ends with concluding remarks on the importance of stochastic VRPs.

Deterministic VRP

Capacitated Vehicle Routing Problem

In its basic form, the Capacitated Vehicle Routing Problem (CVRP) can be dened by a complete undirected graph G = (V, E) where V = {0, 1, . . . , i, . . . , n} and E = {[i, j] ∀i, j ∈ V | i < j} are the vertex and the edge sets respectively. Moreover, let V c = V \ {0} be the customers subset. Each customer has a non-negative demand q i . Vertex 0 is a depot where is located a set of homogeneous vehicles with limited capacity Q. Furthermore, each edge [i, j] ∈ E has a non-negative cost c ij . The objective of the CVRP is to build a set of routes with minimum cost considering that each route starts and ends at the depot, the vehicle capacity Q must be respected, and no split deliveries are allowed. Moreover, a generic route r is dened as an ordered sequence of nodes r = {r 0 , r 1 , . . . , r j , . . . , r k , r k+1 } where r j represents the j th visited node. Each vehicle starts and ends its route at the depot, therefore, r 0 = r k+1 = 0 for every route. Even more, each route r has an associated cost C r = k j=0 c rj ,rj+1 .

Several other extensions and variants for the CVRP have been proposed aiming to bring the models closer to real life applications. Among these, one can nd the Distance Constrained VRP (DVRP) which limits the total distance traveled by each vehicle to a threshold [19,5]; the Heterogeneous VRP (HVRP) where the eet of vehicles is, as its name states, heterogeneous (in terms of capacity or costs) [7,[START_REF] Li | A multistart adaptive memory-based tabu search algorithm for the heterogeneous xed eet open vehicle routing problem[END_REF][START_REF] Penna | An iterated local search heuristic for the heterogeneous eet vehicle routing problem[END_REF]; the Multi Depot VRP (MDVRP) which considers multiple depots where the vehicles start and end their routes [START_REF] Renaud | A tabu search heuristic for the multi-depot vehicle routing problem[END_REF][START_REF] Vidal | A hybrid genetic algorithm for multidepot and periodic vehicle routing problems[END_REF]; the Periodic VRP (PVRP) that requires repeated visits to customers [57]; the open VRP (OVRP) which does not require vehicles to return to depot after serving the last customer [START_REF] Li | The open vehicle routing problem: Algorithms, large-scale test problems, and computational results[END_REF][START_REF] Li | A multistart adaptive memory-based tabu search algorithm for the heterogeneous xed eet open vehicle routing problem[END_REF]; the Orienteering Problem (OP) where customers have an associated prot (or score) collected by a xed size eet not necessarily sucient to visit all the customers, and the objective is to maximize the total prot [START_REF] Vansteenwegen | The orienteering problem: A survey[END_REF][START_REF] Labadie | Hybridized evolutionary local search algorithm for the team orienteering problem with time windows[END_REF]80]. Other variants include more additional constraints such as the Pickup and Delivery Problems (PDP) where people or goods must be transported from dierent origins to dierent destinations [START_REF] Parragh | A survey on pickup and delivery problems[END_REF][START_REF] Parragh | A survey on pickup and delivery problems[END_REF]17,18]. The reader is referred to the mentioned bibliography and to Toth and Vigo [START_REF] Toth | Vehicle Routing[END_REF] for further details on variants of the CVRP.

CVRP Mathematical formulations

Three formulations are mainly used to model the CVRP [START_REF] Semet | Chapter 2: Classical Exact Algorithms for the Capacitated Vehicle Routing Problem[END_REF][START_REF] Laporte | What you should know about the vehicle routing problem[END_REF], the vehicle ow, the commodity ow and the set partitioning formulations. The vehicle ow formulation uses integer variables

x ij ∀i, j ∈ V to represent the number of times an edge is used in the optimal solution [START_REF] Laporte | A branch and bound algorithm for the capacitated vehicle routing problem[END_REF][START_REF] Laporte | Optimal routing under capacity and distance restrictions[END_REF].

Model M1CVRP presents a classical two-index network formulation.

M 1CV RP : min

i,j∈V

x ij c ij (2.1)

j∈V c x 0j = 2m (2.2) i<p|i∈V x ip + j>p|j∈V x pj = 2 ∀p ∈ V c (2.3) i∈S,j / ∈S or i / ∈S,j∈S x ij ≥ 2b(S) ∀S ⊂ V c (2.4) CHAPTER 2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS
x ij ∈ {0, 1} ∀i, j ∈ V c (2.5) x 0j ∈ {0, 1, 2} ∀j ∈ V c (2.6) In M1CVRP the objective (2.1) minimizes the total costs associated with the used edges. Constraint (2.2) determines the degree of the depot, using m as the number of vehicles (m can be a variable). Constraint (2.3) ensures that any vehicle visiting a customer must leave to another node. Constraint (2.4) serves to guarantee capacity restrictions as well as to prevent subtours formation, that is, ensemble of connected customers without being linked to the depot. Practically, the term b(S) can be set to i∈S qi Q , therefore, b(S) is a lower bound on the number of vehicles needed to satisfy the demand of the subset of customers S. Furthermore, constraints (2.5) and (2.6) stand for variables nature.

The second formulation, called commodity ow formulation makes use of continuous variables to model the amount of vehicle load and empty space on the vehicle, when an edge is used. The reader is referred to Baldacci et al. [9] for a complete formulation. The third formulation is the set partitioning one. This one relies on the enumeration of all feasible routes which are then selected through a set partitioning problem [11]. Associated to each route a binary variable serves to decide if it is included within the solution or not. The reader is referred to Laporte [START_REF] Laporte | What you should know about the vehicle routing problem[END_REF] for the complete model.

Other formulations beside the three described can be used. For instance, three-index formulations add an index to identify each vehicle separately. This types of models are useful when particular characteristics of each route aect its feasibility or cost. One of the most common VRP extension modeled by three-index formulation is the one with Time Windows (VRPTW), where time constraints are imposed for the customers visits. It is now explored in more depth.

VRP with Time Windows

The Vehicle Routing Problem with Time Windows (VRPTW) is one of the most important and wellstudied VRPs. Usually the VRPTW uses an extended graph G = V , A . The vertex set V includes an exact copy of the depot node called n + 1, therefore, V = V ∪ n + 1. Moreover A stands for the arcs set dened as as A = {(i, j) ∀i, j ∈ V | i = j}. In addition, each arc (i, j) takes a time t ij ∀ (i, j) ∈ A to be traversed. Also, each customer requires a time t i ∀ V c to be served.

The VRPTW extends the CVRP by dening a time window [e i , l i ] ∀i ∈ V c . The vehicle must start to service the customer during this lapse of time. Time windows constraints can be dened as hard or soft [43]. In the hard version, the vehicles cannot start their services outside the time windows. Nevertheless, early arrivals, i.e. arriving before e i , are possible but vehicles must wait until the opening of the time window. Soft version allows services outside the time window at the expense of a penalization cost. Besides, the depot often has a time window [e 0 , l 0 ] representing the earliest departure time and the latest arrival time for vehicles to the depot. The time window is the same for the depot copy n + 1, i.e. [e n+1 , l n+1 ] = [e 0 , l 0 ]. These additional constraints on service start times add another layer of complexity when compared to the classical CVRP. When the number of vehicles is xed, even computing a feasible solution is NP-Hard [START_REF] Savelsbergh | Local search in routing problems with time windows[END_REF]. The objective of the VRPTW can dier from that of the CVRP. According to Desaulniers et al. [43] exact approaches to the VRPTW usually consider the same objective function as in the CVRP, i.e. the total cost of the routes. Meanwhile, heuristic and metaheuristic methods are often designed to rst minimize the number of required vehicles then the total cost in a hierarchical way.

Mathematical formulations

Similar to the CVRP there exist many formulations for the VRPTW. M2VRPTW presents a threeindex formulation for the version with hard time windows. In this, the binary variable x ijl ∀ i, j ∈ V , l ∈ L takes value one if the vehicle l uses the arc (i, j). Besides, T il stands for the time when x ijl = 1 ∀i ∈ V c

(2.8)

j∈V x 0jl = 1 ∀l ∈ L (2.9) i∈V x ijl - i∈V x jil = 0 ∀j ∈ V c , l ∈ L
(2.10)

i∈V c ∪0
x i,n+1,l = 1 ∀l ∈ L

(2.11)

T il + t i + t ij -T jl ≤ M (1 -x ijl ) ∀l ∈ L, (i, j) ∈ A (2.
12)

e i ≤ T il ≤ l i ∀l ∈ L, i ∈ V (2.13) i∈V c q i j∈V x ijl ≤ Q ∀l ∈ L (2.14)
x ijl ∈ {0, 1} ∀ (i, j) ∈ A, l ∈ L

(2.15)

T il ∈ + ∀i ∈ V , l ∈ L (2.16)
Constraint (2.8) guarantees that each customer is served by only one route. Meanwhile, constraint (2.9) species that all vehicles must leave the depot. In addition, constraint (2.10) guarantees that a vehicle visiting a customer must leave to another node. Moreover, constraint (2.11) states that all vehicles nish their route at node n + 1. Note that in this formulation, the number of eectively used vehicles can be less than |L| as far as variable x 0,n+1,l can take value 1. Furthermore, constraints (2.12) to (2.13) ensure that the times when services start, respect the nodes time windows. In constraint (2.12) term M is a large value that let the equation holds when x ijl takes value zero. Vehicles capacity constraint is guaranteed by (2.14). Last but not least, constraints (2.15) and (2.16) stand for the variables nature.

Solution Methods

Solving the VRPs is not an easy task, since the CVRP is an NP-Hard problem and most of its variants are NP-Hard as well, including the uncertain VRPs. Solution approaches can be classied according to the nature of the solution. In this vein, exact methods guarantee that the optimal solution will be found but at the expense of a prohibitive running time even for medium size instances. Approximate methods on the other hand usually provide quickly a solution but this last can be not optimal. Figure 2.1 provides a simple scheme to VRP solutions approaches. The scheme is not exhaustive but allows to navigate through the extensive amount of methods.

CHAPTER

Exact methods

Exact approaches for the VRPs are highly related to the way the problem is modeled. Dierent methods have been used for solving the CVRP and the VRPTW such as the Branch-and-price (BP), Branch-and-cut (BC) and Branch-and-price-and-cut (BPC) [58,43]. BC has been mainly based on the two-index formulation (M1CVRP) or an analogous version for the VRPTW [12,90]. Overall, the BC works by relaxing integrality constraints and discarding the set of constraints represented by (2.4) for the CVRP and (2.12) to (2.14) for the VRPTW. The BC solves the relaxed problem and identies any subset of variables that violates the removed constraints. If this set is found, it generates the violated constraints, add them to the problem and reiterates. Moreover, when no constraints are identied, the BC branches on a fractional variable and creates problems that are solved with the same approach.

BPC works similarly to BC methods. The dierence relies in the fact that each subproblem (usually the Shortest Path Problem with Resource Constraints) relaxation is solved by means of a column-generation approach. This last approach exploits the set covering formulation for both the CVRP and the VRPTW. Dierent versions of BPC algorithms exist since dierent approaches can be used to solve the subproblems (e.g. ng-routes, q-routes, bidirectional search), or the types of cuts (constraints) added during the iterations. BPC has shown to be the state-of-the-art to solve the CVRP [8,58,[START_REF] Pecin | Improved Branch-Cut-and-Price for Capacitated Vehicle Routing[END_REF] as well as the VRPTW [START_REF] Kohl | 2-path cuts for the vehicle routing problem with time windows[END_REF]42,87]. Nevertheless, Baldacci et al. [10] has proposed the best method based on a reduced set partitioning for the VRPTW. Further analysis and description of the methods are available in [89,43].

Approximate methods -Heuristic and Metaheuristics

While exact methods have seen an incredible development in the last years, the combinatorial nature of the VRPs limit their use to relatively small instances. Nowadays, for example the CVRP can be consistently solved for problems with up to 200 customers [START_REF] Pecin | Improved Branch-Cut-and-Price for Capacitated Vehicle Routing[END_REF]. Meanwhile, VRPTW is consistently solved for instances with up to 100 customers [10]. Still, since applications can easily overpass this size, heuristics and metaheuristics are omnipresent in the literature. Heuristics are approximate algorithms which try to nd good solutions in competitive running times.

Laporte and Semet [START_REF] Laporte | The vehicle routing problem. chapter Classical Heuristics for the Capacitated VRP[END_REF] classify heuristics under constructive and two-phase methods. Labadie et al. [START_REF] Labadie | Metaheuristics for Vehicle Routing Problems[END_REF] also follow this classication. In general words, constructive heuristics work by creating an initial solution that can be further improved [START_REF] Laporte | Chapter 4: Heuristics for the Vehicle Routing Problem[END_REF]. Clarke and Wright savings [39] is by far the most 2.1. DETERMINISTIC VRP known heuristic due to its simplicity [START_REF] Laporte | The vehicle routing problem. chapter Classical Heuristics for the Capacitated VRP[END_REF]. The heuristic works by iteratively merging pairs of routes into single ones, provided that this implies a saving and guarantees feasibility. The process is repeated until no further merges are feasible or the maximum possible saving is negative. Further information on basic heuristics can be found in Toth and Vigo [START_REF] Toth | The vehicle routing problem[END_REF].

Several constructive heuristics tailored for the VRPTW are introduced by Solomon [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF]. By far the most important and successful is the insertion heuristic called I1. Iteratively I1 creates routes starting with some seed customers. These are selected from dierent rules, such as the farthest not visited customer or the one with the earliest initial time window. Once a seed has been selected, I1

calculates an insertion (of non visited customers) criteria based on the classical savings (distance) and the extra time required by adding the new customer. The best customer not yet visited is added in the best possible position. The algorithm iterates until no customers can be added to the current route, then it starts a new one until all customers are visited. The reader is referred to Bräysy and

Gendreau [31] for other constructive approaches.

Two-phase methods can be further divided into cluster-rst, route second and route-rst, cluster second. The rst one is based on the idea of creating groups or clusters of customers respecting the capacity constraint. Then, customers in the cluster are ordered to completely dene a route, by solving a Traveling Salesman Problem (TSP) for each group. Several approaches can be used to create clusters: the sweep algorithm [70] uses angular sectors from the depot to create the necessary clusters.

Fisher and Jaikumar heuristic [56] also uses the idea of clusters around seeds aiming to minimize the distance from customers to cluster seeds.

The route-rst, cluster second approach is mainly based on the idea of creating a giant tour (TSP tour) without considering capacity or other constraints, and then splitting it into feasible routes. This approach was introduced by Beasley [14] and has received more attention since the work of Prins [START_REF] Prins | A simple and eective evolutionary algorithm for the vehicle routing problem[END_REF]. Indeed, Prins showed that route-rst, cluster second algorithms could be as ecient as methods relying on classical methods at the date, such as the Tabu Search. Further examples on the routerst, cluster second can be found in Labadie et al. [START_REF] Labadi | A memetic algorithm for the vehicle routing problem with time windows[END_REF] with an application to the VRPTW, Prins et al. [START_REF] Prins | Tour splitting algorithms for vehicle routing problems[END_REF] addressing the Capacitated Arc Routing Problem and the CVRP, Mendoza et al. [START_REF] Mendoza | A memetic algorithm for the multi-compartment vehicle routing problem with stochastic demands[END_REF] in a multi-compartment vehicles with uncertain parameters problem, Velasco et al. [START_REF] Velasco | A non-dominated sorting genetic algorithm for a bi-objective pick-up and delivery problem[END_REF] with a multi objetive pick-up and delivery problem, and Mendoza et al. [START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF] in the context of a CVRP with stochastic demands. A recent review on on the route-rst, cluster second approach can be found in Prins et al. [START_REF] Prins | Order-rst split-second methods for vehicle routing problems: A review[END_REF].

Although heuristics commonly provide a good trade-o between eciency and quality, they are usually coupled with local search or improvement procedures. The underlying concept of local search is the denition of neighborhoods [31,[START_REF] Labadie | Metaheuristics for Vehicle Routing Problems[END_REF]. These lasts are considered as close related solutions to a generic solution s. Neighborhoods are structured in a way such that movements can be performed on s to achieve a new solution s . Usually, this type of procedures contain two types of movements, namely intra-route and inter-route ones [START_REF] Laporte | Fifty years of vehicle routing[END_REF]. The rst one aects only one route trying to improve it while the second considers and changes more than one route. Moreover, neighborhoods can be wholly explored to select the best improving movement (best acceptance) or partially explored until a movement improves the current solution (rst acceptance) [START_REF] Osman | Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problem[END_REF]. The exploration of neighborhoods is performed in an iterative way, until reaching a stopping condition or achieving a local (global) optimal solution. Among the most used neighborhoods one can nd the k-opt movements [START_REF] Lin | An eective heuristic algorithm for the traveling-salesman problem[END_REF] for which 2-Opt and 3-Opt are the most popular cases, the b-cyclic, k-transfer scheme [START_REF] Thompson | Cyclic transfer algorithm for multivehicle routing and scheduling problems[END_REF], Or-opt movements [START_REF] Or | Traveling salesman-type combinatorial problems and their relation to the logistics of regional blood banking[END_REF] and λ-interchange movements [START_REF] Osman | Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problem[END_REF]. More complex movements can also be found in the literature, e.g. GENI exchange [62] and ejection chains [START_REF] Rego | A Parallel Tabu Search Algorithm Using Ejection Chains for the Vehicle Routing Problem[END_REF]73]. For further details on these neighborhoods and their characteristics, the interested reader is referred to [59,31]. Furthermore, details on ecient implementations of evaluation tests allowing to know either these movements are feasible or not, and to compute the extra cost generated, can be found in [START_REF] Kindervater | Vehicle routing: handling edge exchanges[END_REF][START_REF] Vidal | Timing problems and algorithms: Time decisions for sequences of activities[END_REF].

Despite the success of local search procedures, they are often trapped in local optima when they CHAPTER 2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS are applied to only one initial solution obtained with a constructive heuristic. To overcome these problems, metaheuristics are a good option. According to Bianchi et al. [24] metaheuristics are high level procedures that combine heuristics in a more general framework. Osman and Laporte [START_REF] Osman | Metaheuristics: A bibliography[END_REF] dened them as A metaheuristic is formally dened as an iterative generation process which guides a subordinate heuristic by combining intelligently dierent concepts for exploring and exploiting the search space, learning strategies are used to structure information in order to nd eciently nearoptimal solutions. Originally, metaheuristics were easily identied and dierentiated, however, the increasing hybridization of such methods have blurred the lines between them [START_REF] Laporte | Chapter 4: Heuristics for the Vehicle Routing Problem[END_REF]29].

Still, according to Laporte et al. [START_REF] Laporte | Chapter 4: Heuristics for the Vehicle Routing Problem[END_REF] metaheuristics can be classied into trajectory and populationbased methods. In the rst one, a solution move to another by searching in a neighborhood. Meanwhile, population-based methods use a set of solutions that interact to improve the solution quality.

Within the rst category classication one can nd methods such as the Simulated Annealing (SA) [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], Tabu Search (TS) [72,68], Variable Neighborhood Search (VNS) and Variable Neighborhood Descent (VND) [START_REF] Mladenovi¢ | Variable neighborhood search[END_REF], Deterministic Annealing (DA) [47,[START_REF] Li | Very large-scale vehicle routing: new test problems, algorithms, and results[END_REF], Iterated Local Search (ILS) [13,[START_REF] Lourenço | Iterated Local Search: Framework and Applications[END_REF].

Population-based metaheuristics include Genetic Algorithms (GAs) [START_REF] Reeves | Genetic Algorithms[END_REF][START_REF] Prins | A simple and eective evolutionary algorithm for the vehicle routing problem[END_REF], Ant Colony Optimization (ACO) [START_REF] Reimann | D-ants: Savings based ants divide and conquer the vehicle routing problem[END_REF] and Scatter Search (SS) [71].

Another important trend in solution methods for the VRP are the matheuristics [44,[START_REF] Maniezzo | Matheuristics: Hybridizing Metaheuristics and Mathematical Programming[END_REF]. Matheuristics work by hybridizing heuristics (or metaheuristics) and exact algorithms such as Integer Programming (IP). These methods cooperate and share information to improve the solution quality.

An interesting approach is the Petal heuristics which make use of the set partitioning model. In this one, promising routes (often called petals [START_REF] Laporte | Fifty years of vehicle routing[END_REF][START_REF] Labadie | Metaheuristics for Vehicle Routing Problems[END_REF] in this context) are added to a set partitioning problem as it is done in a column generation scheme [11,[START_REF] Ryan | Extensions of the petal method for vehicle routeing[END_REF][START_REF] Renaud | A tabu search heuristic for the multi-depot vehicle routing problem[END_REF]. A similar approach is used by Mendoza et al. [START_REF] Mendoza | A multi-space sampling heuristic for the vehicle routing problem with stochastic demands[END_REF][START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF], by constructing several solutions from constructive heuristics and then solving a set partitioning problem. Further information on solution methods for dierent variants of VRP can be found at [28,67,32,69,[START_REF] Labadie | Metaheuristics for Vehicle Routing Problems[END_REF].

VRP under uncertainties

A majority of the studies on the dierent VRP variants have been carried out under the assumption that all relevant information is known with certitude when solving the problems. That is, problems are solved in a deterministic static way [START_REF] Pillac | A review of dynamic vehicle routing problems[END_REF]30]. Nevertheless, dealing with real life applications implies the appearance of uncertainties. These lasts arise from many reasons, e.g. weather conditions, accidents, customer presence, etc. Neglecting the uncertainties is not always an option since the eects of variability can have important consequences on solutions quality. Actually, it has been shown that deterministic solutions 2 can lead to systematically bad solutions in an uncertain environment [START_REF] Louveaux | An Introduction to Stochastic Transportation Models[END_REF][START_REF] Sörensen | A practical approach for robust and exible vehicle routing using metaheuristics and monte carlo sampling[END_REF].

Talking about uncertainties is a discussion highly related with information. Certainly, the quality of information and the times when it is available play a major role in both, the models and solution approaches for this type of problems. Moreover, what can be known about the uncertain parameters denes the framework in which an uncertain VRP can be tackled. Pillac et al. [START_REF] Pillac | A review of dynamic vehicle routing problems[END_REF] dene four categories due to sub levels of information characteristics, i.e. information evolution and information quality. Their taxonomy is presented in table 2.1 by considering general uncertain inputs. 2 Based for example on the expected value of the variable parameters The top-left box stands for Static and deterministic problems. This category represents the problems where all parameters are known beforehand, and their values are certain. Moreover, once the solution is deducted it remains unchanged as far as no new information arise. The CVRP and VRPTW lie within this category. In Dynamic and deterministic problems, part or all the inputs is revealed dynamically during the execution of the routes [START_REF] Pillac | A review of dynamic vehicle routing problems[END_REF]. More importantly, there is no exploitable information about the dynamic parameters. Therefore, solutions are construct in an online way.

Static and dynamic uncertain problems (top-right and bottom-right boxes) share the fact that at least one input (parameter) is not known for sure. Same as the dynamic deterministic case, parameters true value is revealed at some specic moments, e.g. when the vehicle arrives at the customer.

Nevertheless, in uncertain problems there exist exploitable information about the parameter such as its probability distribution, the interval set for its value, moments of the distribution, etc. Therefore, static and dynamic uncertain problems can exploit this information to devise their solutions. There are two main dierences between static and dynamic uncertain problems. First, in static problems all inputs (even uncertain ones) are dened before solving the problem while this is not the case in dynamic ones. Second, both categories dier in how solutions are created and treated. In the static uncertain problems, solutions are created before the realization of the uncertain parameters and are barely modied. Meanwhile, dynamic problems construct and change the solution as new information arise.

The framework to handle the uncertainties (static or dynamic) clearly depends on the type of problem tackled and how the information evolves. Additionally, how the uncertainties are modeled plays a major role in the approaches to solve uncertain VRPs. Three main approaches have been used

to deal with uncertainties in VRPs, namely, stochastic, robust optimization, and Fuzzy logic. They are now explored in further detail.

Stochastic VRP

Stochastic programming has been the main paradigm to deal with uncertain VRPs, driving the development of what is known as Stochastic Vehicle Routing Problems (SVRP) [64]. The reader is referred to the work of Birge and Louveaux [27] as a good entry point to the stochastic programming eld.

The main characteristics of SVRPs is the fact that uncertainty alters the condence on the problems parameters, which are modeled as random variables. Thus, in static SVRPs, information is assumed to be prior available to characterize the random parameters of the problem i.e. the probability distributions of the random parameters. As long as parameters are not certain, some constraints might not be fully satised. Whenever a constraint is violated (due to the variability of the parameters) a failure occurs.

SVRPs can be tracked down to the pioneer work of Tillman [START_REF] Tillman | The multiple terminal delivery problem with probabilistic demands[END_REF], who deals with the VRP with stochastic demands (VRPSD) considering multiple depots. The VRPSD is an extension of the classical CVRP where demands are modeled as random variables. Due to this fact, failures can take place whenever a customer' demand is higher than the available remaining capacity in the vehicle. Aiming to overcome and solve the SVRPs, two main paradigms have been mainly used [63]: the a priori (static) and reoptimization (dynamic) paradigms.

CHAPTER 2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS

The a priori optimization works by constructing a solution before the realization of the random variables [86,63]. Therefore, a priori paradigm can be viewed as a two-stage approach. After the prior solution is built (rst-stage) the random parameters are revealed (second-stage). The way in which this information is available depends on the problem. For example, many VRPSD formulations [38,76,60,[START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF] assume that the customer demand is only known when the vehicle arrives at the customer location. Besides, two main approaches are commonly used to optimize the a priori SVRPs, the Stochastic Problem with Recourse (SPR) and the Chance Constraint Programming (CCP). SPR use recourse which can be dened as policies or rules and actions that adjust the prior solution to deal with specic situations such as failures in the second stage. Indeed, recourse can be used to respond to failures. For example, the so called classical VRPSD recourse consists in a return to the depot to load (unload) when a customer demand exceeds the current capacity of the vehicle and restarts the route from this customer. It shall be noticed that recourse actions usually generate a cost which is properly considered in the objective function. Thus, the SPR minimizes the rst-stage cost of the solution plus the expected cost of the recourse (second-stage). CCP works by bounding the probability of failure to a threshold. CCP can be used when recourse can be hardly dened [START_REF] Stewart | Stochastic vehicle routing: A comprehensive approach[END_REF] or to guarantee a service level [START_REF] Li | Vehicle routing problems with time windows and stochastic travel and service times: Models and algorithm[END_REF]52,[START_REF] Zhang | A stochastic vehicle routing problem with travel time uncertainty: Trade-o between cost and customer service[END_REF]. A CCP and SPR formulations for the VRPSD are introduced in the next sub-section 2.2.1.

The reoptimization paradigm on the other hand does not rely on a prior (static) solution. Conversely, it constructs and changes the solution as new information arises (dynamic). The increasing amount and development of Information and Communication Technologies (ICT) has enabled to tackle and evaluate problems in such a dynamic way [START_REF] Pillac | A review of dynamic vehicle routing problems[END_REF]. Although this paradigm usually improves the solutions quality when compared to its static counterpart, it also poses new challenges in terms of computational eciency at operational level. Indeed, the speed in which solutions must be computed is a limitation on such methods. Besides, a priori (static) approaches conduce to stable tactical routes, which are operationally desirable [63]. Moreover, the a priori approach is preferable when anticipating uncertainties is important for routes feasibility and to avoid costs (economic and reputational) [15]. The main focus of the SVRPs is the static a priori approach, further information on dynamic stochastic problems can be found at Pillac et al. [START_REF] Pillac | A review of dynamic vehicle routing problems[END_REF], Bekta et al. [15], and Psaraftis et al. [START_REF] Min | Dynamic vehicle routing problems: Three decades and counting[END_REF].

Given the static and stochastic context, the SVRPs literature is frequently classied with respect to the parameters stochasticity [64]. The most common studied versions are the VRP with stochastic demands, with stochastic customers, and with stochastic times. They are now further explored.

VRP with stochastic demands

Among the SVRP, the VRP with stochastic demands (VRPSD) is the furthermost studied problem.

The rst reported solution method for the VRPSD is proposed by Tillman [START_REF] Tillman | The multiple terminal delivery problem with probabilistic demands[END_REF] who used a modication of the Clarke and Wright heuristic [39]. The general VRPSD with recourse is described in the model M3VRPSD while the CCP formulation is presented in model M5VRPSD. Both models are proxies for the SVRPs. The notation of Gendreau et al. [64] is used. Dierences among M3VRPSD and M1CVRP are the inclusion of the expected recourse cost (Q (x)) in objective function (2.1) and constraint (2.4) which need to consider the demand expected value. Since demands are modeled as random variables, they become q i ∀i ∈ V c with expected value E [ q i ] = µ i ∀i ∈ V c . The expected recourse cost is determined by Q (x) [START_REF] Laporte | An integer l-shaped algorithm for the capacitated vehicle routing problem with stochastic demands[END_REF].

M 3V RP SD min

i,j∈V

x ij c ij + Q (x, q i ) (2.17) i,j∈S x ij ≤ |S| - i∈S µ i Q ∀S ⊂ V c , 3 ≤ |S| ≤ n -1 (2.18) (2.2) -(2.

M 4V RP SD min

i,j∈V

x ij c ij (2.20) i∈S,j / ∈S or i / ∈S,j∈S x ij ≥ 2V β (S) ∀S ⊂ V c , 3 ≤ |S| ≤ n -1 (2.21) (2.2) -(2.3) , (2.5) -(2.6) (2.22)
SPR has been the dominant trend in the VRPSD (Bertsimas [22], Gendreau et al. [65,66], Hjorring and Holt [83], Christiansen and Lysgaard [38], Goodson et al. [76], Gauvin et al. [60], Mendoza et al. [START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF]) when compared to the CCP formulation (Stewart and Golden [159,[START_REF] Stewart | Stochastic vehicle routing: A comprehensive approach[END_REF], Dror et al. [46]). The most used recourse policy, the classical recourse is dened as follows. When the load of the vehicle is fullled it returns to the depot to unload the charge, and resumes its assigned route from the failure point [22]. However, other recourse policies have been studied and implemented through several studies: preventive restocking policies are extensions of the classical recourse where return trips to the depot are performed even if the vehicle is not empty to avoid future failures [START_REF] Yang | Stochastic vehicle routing problem with restocking[END_REF]23,25,[START_REF] Luo | Adaptive large neighborhood search heuristics for the vehicle routing problem with stochastic demands and weight-related cost[END_REF][START_REF] Zhang | On-time delivery probabilistic models for the vehicle routing problem with stochastic demands and time windows[END_REF]; pairing strategies allowing the cooperation of multiple vehicles [3]; split deliveries between paired routes [START_REF] Lei | The vehicle routing problem with stochastic demands and split deliveries[END_REF] in which some customers are served by two vehicles; and backup routes [50] that receive customers from primary routes. Although the use of more complex recourse policies can generate a signicant saving relative to simpler ones [3], the latter have been preferred since they allow more tractable models and stable tactical routes [63]. 2.2 shows a summary of the literature of VRPSD. Exact methods have been used to solve the VRPSD, for instance Gendreau et al. [65] employed an integer formulation mixed with the L-shaped method [START_REF] Laporte | The integer l-shaped method for stochastic integer programs with complete recourse[END_REF] to optimally solve instances with up to 70 customers given that all customers are present. Still, their method is restricted to discrete probability distributions [65]. Hjorring and Holt [83] tackle the VRPSD with only one route. The classical recourse is used considering two scenarios: the rst one is the normal stockout which represents the case when the vehicle returns to the depot and back to the customer where the failure takes place. The second one is the exact stockout, representing the case when the vehicle has just enough capacity to serve the current customer. In that case, the vehicle returns to the depot and continues towards the next customer. The authors use an L-shaped approach to solve the problem. Furthermore, the authors propose new optimality cuts that improve the performance of the algorithm.

CHAPTER 2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS

An important development in the eld is presented by Christiansen and Lysgaard [38]. The authors proposed the rst branch-and-price algorithm to solve the VRPSD with classical recourse.

Moreover, they determined a testbed set of 40 instances based on Augerat sets A and P. The results allowed to have a common benchmark. Moreover, although the method solved problems with up to 60 customers, it showed to be more powerful with a higher number of vehicles when compared to the L-shaped method [65,83,[START_REF] Laporte | An integer l-shaped algorithm for the capacitated vehicle routing problem with stochastic demands[END_REF]. Gauvin et al. [60] developed an improvement of Christiansen and Lysgaard work introducing a Branch-and-cut-and-price. The column generation is accelerated using a TS heuristic and a bidirectional labeling algorithm. The BCP outperforms the BP of Christiansen and Lysgaard [38]. In fact, 20 more instances are solved to optimally, so 38 of 40 instances are closed.

Furthermore, the method achieves to solve problems with up to 100 customers in less than 20 minutes.

Both works assume independent Poisson distributions. More recently, Biesinger et al. [25] propose to use the L-shaped algorithm for the Generalized VRP with stochastic demands. More interesting is that the recourse action corresponds to a preventive restocking policy [START_REF] Yang | Stochastic vehicle routing problem with restocking[END_REF] in contrast to the classical one. Therefore, this is is the rst exact method that considers such a dierent recourse. Results are presented for small size instances with dierent types of variations. The solution method is proven to be eective for instances involving up to 40 nodes and with less than three expected restocking actions.

Due to the complexity of the VRPSD, approximate methods have been mostly designed to solve real size instances. Stewart and Golden [START_REF] Stewart | Stochastic vehicle routing: A comprehensive approach[END_REF] proposed a CCP and two SPR for the problem. The SPR are dierentiated by the penalty induced by the recourse. In the rst one, a penalty is taken into account disregarding the amount of the violation (lack of capacity). The second one induces a penalty proportional to the amount of the violation. Using adaptations of the Clarke and Wright [39] and the Generalized Lagrange Multipliers [START_REF] Stewart | A vehicle routing algorithm based on generalized lagrange multipliers[END_REF] heuristics, the authors solve some instances assuming the demands are independent normal variables. Moreover, the authors present an approach to transform the CCP formulation into an equivalent deterministic CVRP model. Also, solutions are provided for the correlated-demands case. Dror and Trudeau [46] consider the models presented in [START_REF] Stewart | Stochastic vehicle routing: A comprehensive approach[END_REF] and prove that in the VRPSD framework, the direction of the routes can have an important impact on the objective function. This important fact shows the dierent structure that SVRPs take when compared to deterministic VRPs. The Clarke and Wright [39] heuristic is adapted to incorporate the expected recourse cost. More important, when two routes are being considered to merge, the direction of the resulting route must be evaluated since this impacts the objective function. A comparison to the adaptation of [START_REF] Stewart | Stochastic vehicle routing: A comprehensive approach[END_REF] shows that their method performs better in terms of the number of vehicles, deterministic length and expected cost.

Bertsimas [22] considers the VRPSD as an SPR. The authors propose two strategies (recourse) depending on the available information. Strategy a on the one hand, is the classical recourse policy previously described. Strategy b on the other hand, assumes that demands are revealed before the tours (routes) start, thus, customers with zero demand are omitted. Although, the last scenario seems CHAPTER 2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS more suitable for a reoptimization paradigm, strategy b is interesting when resources do not allow to perform a reoptimization scheme e.g when computational resources are scarce. Furthermore, this study proposes closed-form expressions to calculate the expected costs of routes under general probabilistic assumptions. Besides, under certain conditions (customer locations) the a priori strategies perform very closely to the reoptimization paradigm. Further works of Gendreau et al. [66] led to the development of a tabu search in the context of the VRPSD where additionally, customers might or not be present. To calculate the expected cost, a proxy function is used to approximate the expected solution costs.

Yang et al. [START_REF] Yang | Stochastic vehicle routing problem with restocking[END_REF] extended the classical recourse policy in their work: they assume that vehicles return to depot when capacity is fullled but they can return before this event happens. When leaving a customer, the capacity of the vehicle is compared to a threshold to determine if a visit to the depot is valuable before attending the next customer. One threshold is associated with each customer given a xed route. Moreover, this restocking policy is proven to be optimal given a xed route.

Two heuristics are proposed for the VRPSD with route durations limits for instances with discrete triangular distributed demands. Bianchi et al. [23] consider the same restocking policy proposed by Yang et al. [START_REF] Yang | Stochastic vehicle routing problem with restocking[END_REF] for the VRPSD. Their objective is to compare multiple metaheuristics, each of them under two approaches for the local search: using the solution routes representation and a TSP representation. The rst one evaluates the pertinence of local search movements given the cost changes in the routes. The second one uses the deterministic tour length of an underlying TSP representation of the solution to estimate if a movement is performed or not. Overall, the TSP approach can be seen as an acceleration technique to improve the execution of the local search since recourse costs are not calculated. Several metaheuristics are tested on instances with up to 200 customers.

Goodson et al. [76] approach the VRPSD using a cyclic-order representation of the solutions. The classical recourse action is considered and demands are assumed to be Poisson distributed [38,60]. An ecient way to deal with neighborhoods under a solution cyclic-order representation is introduced. This type of representation is based on a cyclic permutation of numbers, which is then decoded into detailed routes. The authors use this type of representation since small changes in the permutation can conduce to important changes in the decoded routes. Simulated Annealing is selected by the authors to tackle the problem. The method achieves 16 out the 18 optimal solutions presented by Christiansen and Lysgaard [38]. Later, Mendoza et al. [START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF] investigated the VRPSD with route durations constraint. A hybrid metaheuristic composed of a GRASP plus Heuristic Concentration (HC) are used. The solution method operates several simple heuristics to construct a pool of routes, then, a set partitioning problem is solved giving a nal solution. The method solves the VRPSD with classical recourse and achieves to nd all the best-known solutions of the Christiansen and Lysgaard [38] testbed with average gaps of 0.02%. For the version with duration constraints two models are considered, a CCP and an SPR. Ak and Erera [3] proposed an alternative recourse policy called Paired locally coordinated (PLC).

The PLC extends the idea of usual recourse actions which consider vehicles separately to paired vehicles. When the rst vehicle (type I) achieves its capacity it returns to the depot nishing its service. Then, the second vehicle (type II) incorporates the unserved customers at the end of its route. Besides, if the second vehicle faces a failure, it uses the classical recourse. The aim of this type of coordination is to improve the expected total cost. A Tabu Search is used to solve instances with up to 150 customers where customers have homogeneous discrete demands. Results show that comparing to the classical recourse, the PLC allows signicant savings, ranging from 3% to 25% when instances with 50 or more customers are considered. Erera et al. [50] worked on the VRPSD with hard time windows. A limit on route duration is also considered, so drivers return to the depot respecting the working hours. Although deliveries and analyzed data comprises a whole week, routes are created for specic days. One of the main characteristics of the work is the fact that routes must if possible visit the same customers. The motivation for this idea is to create long-term-relationships with the 2.2. VRP UNDER UNCERTAINTIES customers along the horizon. Customers are classied into two groups, one that can be visited by any driver and another group that must be visited at most by two drivers within a week. A CCP formulation is used to constraint the probability of missing a customer' time window. The problem is solved by means of a heuristic which uses Monte Carlo simulation to verify constraints satisfaction.

A novel recourse approach is used by creating primary and backup routes. The latter can receive customers to guarantee feasibility or improve costs. A case study is presented and results are shown for it.

Erera et al. [51] address the version of the VRPSD with duration constraints. Two recourses are considered, the classical recourse and a variant where the vehicle performs a return to the depot without servicing any demand if its current load is less than the demand met at the last customer, then it comes back to this customer and satises its entire demand. To guarantee the route duration constraint the authors propose to solve an Adversarial Problem. Since recourse actions generate additional time depending on demands realization, the Adversarial Problem seeks to maximize the additional time needed over all the possible demands values. Thus, solving this problem for a given route allows to guarantee that in the worst case scenario, the duration constraint holds. A TS is developed to solve the problem with independent demands which follow a discrete uniform distribution.

Results are presented for instances with up to 100 customers. Moreover, the inuence on the eet size for the VRPSD with duration constraint is discussed. Although no general conclusions can be drawn, the results indicate that duration constraints increase the required eet size. Mendoza et al. [START_REF] Mendoza | A memetic algorithm for the multi-compartment vehicle routing problem with stochastic demands[END_REF] used a Memetic Algorithm (MA) to solve the VRPSD with multi-compartments. As far as the calculus of failure probabilities is expensive, an approximation called Take-all-policy (TAP) is employed. Furthermore, two heuristics are used to compare the eciency of the MA with TPA to solve the problem. The rst one is based on an adaptation of the Clarke and Wright heuristic proposed by Dror and Trudeau [46], the second one solves the deterministic problem with the MA but reduces the capacity of the vehicles by a percentage. Demands are assumed to be independent normally distributed and results are presented for instances with up to 484 customers. The MA mixed with TAP approximation gives the best results in terms of number of vehicles and operational costs.

Other related variants of the VRPSD can also be found in the literature. A weight-related cost version with dynamic recourse strategy is introduced in [START_REF] Luo | Adaptive large neighborhood search heuristics for the vehicle routing problem with stochastic demands and weight-related cost[END_REF] while a multi-objective VRPSD is proposed in [61]. Time dimension has been also considered, either in the form of time windows [START_REF] Lei | The capacitated vehicle routing problem with stochastic demands and time windows[END_REF][START_REF] Nguyen | Satiscing measure approach for vehicle routing problem with time windows under uncertainty[END_REF] or route duration constraints [START_REF] Yang | Stochastic vehicle routing problem with restocking[END_REF]50,51,77,[START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF][START_REF] Nguyen | Satiscing measure approach for vehicle routing problem with time windows under uncertainty[END_REF]. Furthermore, an important version of the VRPSD considering stochastic customers is discussed by Gendreau et al. [65,66].

Regarding the reoptimization paradigm, the work of Dror et al. [45] introduced this approach using a Markov Decision Process (MDP) to solve the single vehicle case. In this vein, single-vehicle cases have been investigated in [START_REF] Secomandi | Comparing neuro-dynamic programming algorithms for the vehicle routing problem with stochastic demands[END_REF][START_REF] Secomandi | A rollout policy for the vehicle routing problem with stochastic demands[END_REF][START_REF] Secomandi | Reoptimization approaches for the vehicle-routing problem with stochastic demands[END_REF]. A multi-vehicle version was rst proposed by Goodson [77] aiming to maximize the expected demand served.

Summing up, heuristics for the VRPSD have been mainly based on the adaptation of heuristics initially designed for the CVRP [START_REF] Tillman | The multiple terminal delivery problem with probabilistic demands[END_REF][START_REF] Stewart | Stochastic vehicle routing: A comprehensive approach[END_REF]46,[START_REF] Mendoza | A memetic algorithm for the multi-compartment vehicle routing problem with stochastic demands[END_REF]. Similarly to the deterministic VRP, metaheuristics have played an important role to solve the VRPSD. Tabu Search (TS) is the most used metaheuristic in dierent variants of the VRPSD [66,3,51,[START_REF] Nguyen | Satiscing measure approach for vehicle routing problem with time windows under uncertainty[END_REF]. Other methods such as the LNS [START_REF] Lei | The vehicle routing problem with stochastic demands and split deliveries[END_REF], SA [76],

Memetic Algorithms (MA) [START_REF] Mendoza | A memetic algorithm for the multi-compartment vehicle routing problem with stochastic demands[END_REF], and hybridized metaheuristics e.g. Multi-Space Sampling [START_REF] Mendoza | A multi-space sampling heuristic for the vehicle routing problem with stochastic demands[END_REF][START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF] are addressed in the literature. Exact methods have been concentrated on L-shaped method [65,83,[START_REF] Laporte | An integer l-shaped algorithm for the capacitated vehicle routing problem with stochastic demands[END_REF]25], branch-and-price [38] and branch-and-cut-and-price [60]. L-shaped methods have shown to be able to solve instances with up to 100 customers but few number of vehicles [START_REF] Laporte | An integer l-shaped algorithm for the capacitated vehicle routing problem with stochastic demands[END_REF]. Branch-andprice and branch-and-cut-and-price on the other hand have been able to solve one instance with up to 100 customers [60]. Nevertheless, the former methods show a better performance when dealing with multiple vehicles in the solutions [38,60]. It shall be noticed that only the work of Biesinger et al. [25] on the Generalized Vehicle Routing Problem with Stochastic Demands (GVRPSD) makes use of exact methods considering a restocking policy, while the others remain on the classical recourse. CHAPTER 2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS Moreover, it is evidenced that many of the works use their own set of instances [START_REF] Laporte | An integer l-shaped algorithm for the capacitated vehicle routing problem with stochastic demands[END_REF][START_REF] Yang | Stochastic vehicle routing problem with restocking[END_REF]23,3,50,51,[START_REF] Mendoza | A memetic algorithm for the multi-compartment vehicle routing problem with stochastic demands[END_REF][START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF][START_REF] Luo | Adaptive large neighborhood search heuristics for the vehicle routing problem with stochastic demands and weight-related cost[END_REF] since no standard benchmark has been proposed. Meanwhile, other works have rely on original or modied well-known data sets such as Nguyen et al. [START_REF] Nguyen | Satiscing measure approach for vehicle routing problem with time windows under uncertainty[END_REF] and [START_REF] Lei | The capacitated vehicle routing problem with stochastic demands and time windows[END_REF][START_REF] Lei | The vehicle routing problem with stochastic demands and split deliveries[END_REF] using the Solomon [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF] instances. Likewise, since the introduction of the Christiansen and Lysgaard [38] benchmark other works [76,60,[START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints[END_REF][START_REF] Luo | Adaptive large neighborhood search heuristics for the vehicle routing problem with stochastic demands and weight-related cost[END_REF][START_REF] Marinaki | A glowworm swarm optimization algorithm for the vehicle routing problem with stochastic demands[END_REF] have used it as a base to compare their results (or improve them).

VRP with stochastic customers

The VRP with stochastic customers (VRPSC) refers to the variant where the presence of the customers is stochastic, i.e. they might be present or not. Among the classical stochastic parameters, this one is the far less studied. The roots of the VRPSC are derived from the Traveling Salesman Problem with Stochastic Customers (TSPSC). TSPSC was rst introduced in the thesis of Jaillet [85]. Both the VRPSC and the TSPSC are modeled as a two-stage problem. In the rst stage, a solution to the TSP is computed, while in the second stage, absent customers are revealed so the routes (tour) follow the original order but skip the absent customers. Bertsimas and Howell [21] reviewed the main results of Jaillet and proposed several heuristics based on classical neighborhoods ideas, e.g. 2-opt and 3-opt, and others on angular sorting and space lling curves. A branch and cut approach is proposed in [START_REF] Laporte | A priori optimization of the probabilistic traveling salesman problem[END_REF] for the TSPSC. Bertsimas [22] also derives heuristics and bounds for the case when demand is binary, i.e. either one or zero, and shows that the TSPSC is a special case of this problem. Further development on the VRPSC is due to a more general problem, the VRPSC with stochastic demands (VRPSCD). A tabu search is proposed in [66] while Gendreau et al. [65] use the L-shaped algorithm to optimally solve the VRPSCD.

VRP with stochastic times

The Vehicle Routing Problems with Stochastic Times (VRPST) are versions of the VRP in which the travel or/and service times are random variables. Table 2.3 presents a summary of some of the most important works in the eld. Compared to the VRPSD, the VRPST is a more recent studied version. The rst literature reference is presented by Laporte et al. [START_REF] Laporte | The vehicle routing problem with stochastic travel times[END_REF] who considered the case with stochastic travel and service times, with both CCP and SPR formulations. The authors use a two-stage approach, designing an a priori solution before travel and service times are revealed. Moreover, in the second stage the vehicles follow the a priori solution and a penalization for late arrivals proportional to the length of the delay is considered. Optimal solutions are found for small instances with up to 20 customers and with at most ve travel time scenarios. An adaptation of the Clarke and Wright [39] is used by Lambert et al. [START_REF] Lambert | Designing collection routes through bank branches[END_REF] for the version with stochastic travel times in a bank money collection context. Similarly to Laporte et al. [START_REF] Laporte | The vehicle routing problem with stochastic travel times[END_REF], a penalization is incurred on late arrivals, in this work, due to the money lost of interest. More recently, Kenyon and Morton [92] tackled the version with stochastic service and travel times under two dierent objectives. The rst one minimizes the maximum completion times of the whole set of routes, while the second one maximizes the probability of completing the routes within a pre-specied deadline. A Generalized Variable Neighbourhood Search (GVNS) for the capacitated VRP with stochastic service times is designed in [START_REF] Lei | A generalized variable neighborhood search heuristic for the capacitated vehicle routing problem with stochastic service times[END_REF] using a penalization on route duration. Service times are assumed to be normally distributed with a coecient of variation (CV) of 0.25. A closed form to calculate the expected delay is presented under normality assumption. GVNS shows better results at the expense of higher computational times in comparison to a VNS and a VND elaborated by the same authors. An interesting study is conducted by Lecluyse et al. [START_REF] Lecluyse | Vehicle routing with stochastic time-dependent travel times[END_REF] for the VRPST with stochastic time-dependent travel times. The authors consider four scenarios depending on two relevant dimensions: speed proles and weather/road conditions. The proposed objective function includes the expected duration of the routes but also its variability. This last is added through the standard deviation, which is multiplied by a parameter β(β ≥ 0). Parameter β is intended to represent the risk preferences of the decision maker. Results are obtained with a TS heuristic on 27

VRP UNDER UNCERTAINTIES

Augerat datasets. The authors conclude that their solution approach gives more robust and reliable routes. In fact, their solutions compensate the increase in the expected route durations by diminishing their variability, especially, when roads are congested for the most part of the day or when weather conditions are bad.

When the VRPST is considered, one fundamental dimension to take into account is the time windows. VRPST with soft time windows have been mainly studied [START_REF] Ta³ | Vehicle routing problem with stochastic travel times including soft time windows and service costs[END_REF][START_REF] Ta³ | The time-dependent vehicle routing problem with soft time windows and stochastic travel times[END_REF][START_REF] Ta³ | Vehicle routing with soft time windows and stochastic travel times: A column generation and branch-and-price solution approach[END_REF] when compared to hard time windows [53,26]. However, many other works enforce the respect of the early time window constraint but allow late services [6,[START_REF] Li | Vehicle routing problems with time windows and stochastic travel and service times: Models and algorithm[END_REF][START_REF] Zhang | A stochastic vehicle routing problem with travel time uncertainty: Trade-o between cost and customer service[END_REF]48,[START_REF] Miranda | The vehicle routing problem with hard time windows and stochastic travel and service time[END_REF][START_REF] Nguyen | Satiscing measure approach for vehicle routing problem with time windows under uncertainty[END_REF], thus mixed soft and hard time windows are used. One of the reasons for the use of soft time windows can be related to the convolution property.

The use of convolution property enables the deduction of closed form expressions to calculate the recourse cost (SPRs) or the probability of constraints violations (CCPs). Hard time windows have a truncation eect on the cumulative times, thus, convolution rules are no longer useful. Aiming to overcome this problem some authors have used dierent techniques to model the arrival times when hard time windows are considered: assuming that they can be modeled a well-known probability distribution [36,48], using discrete approximations [START_REF] Zhang | A stochastic vehicle routing problem with travel time uncertainty: Trade-o between cost and customer service[END_REF][START_REF] Miranda | The vehicle routing problem with hard time windows and stochastic travel and service time[END_REF] to estimate their distribution or modeling travel and service times by phase type distributions [75] and then approximate the distribution arrival times.

Wang and Regan [START_REF] Wang | Assignment models for local truckload trucking problems with stochastic service times and time window constraints[END_REF] A similar problem is tackled in the article of Li et al. [START_REF] Li | Vehicle routing problems with time windows and stochastic travel and service times: Models and algorithm[END_REF]. Both travel and service times are assumed to be independent random variables following a normal distribution. Time windows are considered in the same way as in [6], thus services are not allowed before the opening time window.

Two models are presented: a CCP and a SPR. The rst one bounds the probability of time windows violation on both customers and the depot (route duration). The SPR, penalizes late arrivals regarding time windows and route duration. To solve both formulations a modied TS is used. TS embeds a Monte Carlo simulation to estimate probabilities of failure or the recourse cost. Results for own generated instances with up to 100 customers are discussed. One of the main conclusions of the authors is that CCP is overly constrained and harder to solve than the SPR. As a matter of fact, when constraints on customers and route duration are simultaneously active, their method fails to nd a feasible solution. A generalization of the models presented by Li et al. [START_REF] Li | Vehicle routing problems with time windows and stochastic travel and service times: Models and algorithm[END_REF] is introduced by Zhang et al. [START_REF] Zhang | A stochastic vehicle routing problem with travel time uncertainty: Trade-o between cost and customer service[END_REF] proposing a combined SPR and CCP. The model allows to compare and make trade-os between customers service and the operator costs by changing the associated constraints.

A computational study based on six Solomon [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF] instances restricted to 20 customers is performed.

Travel times are assumed to be log-normal distributed while service times follow a normal distribution.

Coecients of variations range from 0.2 to 0.6 for the travel times and are xed to 0.4 for the service times. The authors use a TS to solve the set of instances. Furthermore, a discrete approximation of the arrival and service start time distributions, called α Discrete Approximation Method (αDAM), is validated and used within the TS. αDAM considers a discrete set of possible values and the associated probability mass function to represent the travel times, as well as the service times. A discrete approximation of the arrival and service start times distributions is recursively derived. A comparison with the normality assumption made by Chang et al. [36] shows the better performance of αDAM to estimate the service levels, expected earliness and tardiness, and expected route costs. Nevertheless, for the sake of fairness it must be mentioned that only one partial route involving 5 customers is tested.

Ta³ et al. [START_REF] Ta³ | Vehicle routing problem with stochastic travel times including soft time windows and service costs[END_REF] deal with the pure soft time windows version of the problem under stochastic travel times. In their work, services can start before or after the time window at the expense of a penalization. The cost incurred by early or late services is dened as the service costs, while the transportation costs account for vehicle travel costs and driver overtime. The two types of costs are weighted in the objective function. Authors use the Gamma distribution to model the travel times with various coecients of variation. To solve the problem a three-phase algorithm based on TS is used. In the rst phase, the authors construct a solution by means of I1 heuristic proposed by Solomon [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF], without considering service costs. The initial solution is further improved by the TS during the second phase. In the third phase, the best solution found by TS is improved by changing the departure time of the vehicles. Since time windows are soft, the departure time can have a signicant eect on total costs, by applying this post-optimization procedure the total costs are reduced on average by 1.3%. Ta³ et al. [START_REF] Ta³ | Vehicle routing with soft time windows and stochastic travel times: A column generation and branch-and-price solution approach[END_REF] propose a column generation within a B&P procedure to optimally solve the problem introduced in [START_REF] Ta³ | Vehicle routing problem with stochastic travel times including soft time windows and service costs[END_REF]. Using a set partitioning formulation, the authors use a labeling algorithm to solve the associated Elementary Shortest Path Problem with Resource Constraints (ESPPRC) [54]. Moreover, since time windows are soft, the only resource constraint is the capacity of the vehicle. Due to this, the authors reduce the capacity of the vehicles for type 1 instances from the classical Solomon instances [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF] to make them tractable. Results are provided for 20, 25, 50 some 100 customers instances. An extension of the work presented in [START_REF] Ta³ | Vehicle routing problem with stochastic travel times including soft time windows and service costs[END_REF] is developed by Ta³ et al. [START_REF] Ta³ | The time-dependent vehicle routing problem with soft time windows and stochastic travel times[END_REF] including stochastic time-dependent travel times (see also [START_REF] Lecluyse | Vehicle routing with stochastic time-dependent travel times[END_REF]). In this work, the authors considered 5 intervals of time within a day in which travel speeds vary due to trac conditions. Assuming that the time to travel a unit distance can be modeled by a Gamma random [START_REF] Miranda | The vehicle routing problem with hard time windows and stochastic travel and service time[END_REF].

A CCP formulation is used so late arrivals are constrained by a desired service level. Early services are not allowed but late services are permitted. The truncated random variables are modeled using a discrete approximation of the arrival and service start times distributions similar to the one of Zhang et al. [START_REF] Zhang | A stochastic vehicle routing problem with travel time uncertainty: Trade-o between cost and customer service[END_REF]. An exhaustive comparison to αDAM is presented showing an improvement in various metrics. Moreover, time consumption is greatly improved. The discrete approximation is embedded within an ILS algorithm to solve some Solomon [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF] instances.

Nguyen et al. [START_REF] Nguyen | Satiscing measure approach for vehicle routing problem with time windows under uncertainty[END_REF] present an interesting approach for the VRPST with stochastic travel times.

A Satiscing Measure Approach (SMA) is derived to evaluate the dissatisfaction of the customers for late arrivals. Late services are allowed and early arrivals are forbidden. Travel times are assumed to follow an ambiguous distribution, though, for simulation purposes it is assumed that they follow a Gamma distribution. The SMA is incorporated into the objective function as main objective (among a hierarchy). By construction of the SMA, this implies a minimization of the expected tardiness at customers. Their approach is embedded within a TS algorithm showing a remarkable performance; moreover, the results show that by incrementing the number of vehicles used, as well as the distance of the routes, the expected tardiness at customers can be signicantly reduced. Errico et al. [53] extended their work in [52] dealing with stochastic travel times and hard time windows. In their more recent work, the authors used a combined SPR and CCP to solve the problem. Penalties are incurred each time a customer is not provided with a service. Furthermore, customers that are not served are picked up based on the two proposed recourse actions. Moreover, customers may be avoided if they induce the route to be operationally infeasible. In such a case, the current or the next customer is selected to be omitted. Also, two constraints are added to ensure that the probability of having no failure at a route reaches a required threshold and limit the number of recourses to at most one.

Another interesting point in this work is the way in which information is revealed. When a vehicle arrives at a customer, it spends some xed time to determine the actual service time. Only then the choice to avoid current or next customer is made. Service times are assumed to be discrete random variables with triangular distribution and two cases regarding the range of possible values are considered for experiments. The problem is modeled as a set partitioning problem and solved through a branch-and-price.

Binart et al. [26] also studied the problem with hard time windows. In fact, their model considers two types of customers, mandatory and optional. Only mandatory customers have an associated time window. Uncertainty aects both the service and travel times which are modeled by discrete triangular distributed random variables. Moreover, multiple depots are considered. A two-phase approach is used to solve the problem. During the rst phase a skeleton of routes serving mandatory customers is created. During the second phase, optional customers are added to the planned routes which are modied in real time to enable time windows to be respected. Furthermore, in the second phase, the routes are evaluated in an SPR and CCP context. Penalties are included if the time window is missed, while the probability of success is also constrained.

VRP UNDER UNCERTAINTIES

In summary, VRPST have been solved mostly with the use of metaheuristic methods because of the complexity of the problem. The solutions approach usually rely on a static approach. Simple heuristic methods have been rarely used [START_REF] Lambert | Designing collection routes through bank branches[END_REF]. Within metaheuristics, approaches such as TS [START_REF] Russell | Vehicle routing with soft time windows and erlang travel times[END_REF][START_REF] Li | Vehicle routing problems with time windows and stochastic travel and service times: Models and algorithm[END_REF]48,[START_REF] Ta³ | Vehicle routing problem with stochastic travel times including soft time windows and service costs[END_REF][START_REF] Zhang | A stochastic vehicle routing problem with travel time uncertainty: Trade-o between cost and customer service[END_REF][START_REF] Ta³ | The time-dependent vehicle routing problem with soft time windows and stochastic travel times[END_REF][START_REF] Nguyen | Satiscing measure approach for vehicle routing problem with time windows under uncertainty[END_REF] are the foremost used, other metaheuristics such as ALNS [START_REF] Ta³ | The time-dependent vehicle routing problem with soft time windows and stochastic travel times[END_REF], GNVS [START_REF] Lei | A generalized variable neighborhood search heuristic for the capacitated vehicle routing problem with stochastic service times[END_REF], GA [6], and approximate methods with Branch-and-cut based frameworks [92,1] have also been used. Exact methods are limited to branch-and-cut [START_REF] Laporte | The vehicle routing problem with stochastic travel times[END_REF], branch-and-price [START_REF] Ta³ | Vehicle routing with soft time windows and stochastic travel times: A column generation and branch-and-price solution approach[END_REF] and branch-and-cutand-price [53]. The use of these approaches is limited to cases where few scenarios are considered [START_REF] Laporte | The vehicle routing problem with stochastic travel times[END_REF], when convolution properties enable closed forms evaluations [START_REF] Ta³ | Vehicle routing with soft time windows and stochastic travel times: A column generation and branch-and-price solution approach[END_REF] or when discrete variables are considered [53].

Robust Optimization

While the SVRPs have been mostly solved using the Stochastic Programming paradigm other approaches in the literature can be found. In fact, stochastic models have been criticized on three main points: the underlying probability distributions must be known in advance, convolutions must be computationally tractable and provided solutions might be infeasible for some realizations of the random variables [START_REF] Solano-Charris | Local search based metaheuristics for the robust vehicle routing problem with discrete scenarios[END_REF]. Actually, if the underlying probability distributions are unknown or information is very limited this precludes the use of SVRPs. Robust Optimization (RO) oers a dierent approach to deal with uncertainties. In fact, uncertainties become deterministic and set-based [20]. The objective is to create solutions that are feasible for any realization of the uncertain parameters rather than being immunized in a probabilistic sense [20]. The interested reader is referred to the book of Ben-Tal et al. [16] and the paper of Bertsimas et al. [20] as excellent introductory points.

Although less studied than SVRPs, robust optimization has also been applied to VRP problems. [1] propose to solve the RVRP with deadlines and travel times uncertainty. Moreover, the authors also tackle the SVRP variant under the same type of uncertainties. For the RVRP a performance measure called lateness index is dened to minimize the risk of violating deadlines, while the SVRP formulation uses a sampling approach similar to the one in [92]. Furthermore, an extension to consider soft time windows is described, in the case of the RVRP this is achieved by adding an earliness index. To solve both formulations, the author use a branch-and-cut approach. Results are reported for instances considering directed graphs with up to 80 customers and 240 arcs. A comparison between the SVRP and the RVRP shows that SVRP heavily resides in the prior knowledge of the underlying probabilistic distributions, while the RVRP outperforms the SVRP if this information is incorrect.

Fuzzy Logic

Fuzzy logic has also been applied to deal with VRP under uncertainties. In fact, given the inherent variability of the problems parameters, fuzzy numbers and variables can be used to model these parameters using the fuzzy logic sets. Fuzzy sets work dierently than deterministic sets. In the former, a number has a degree of pertaining to the set with a grade of membership which ranges between zero and one [4]. Deterministic sets are thus a particular case of the fuzzy ones where partial membership to a set is not allowed, therefore, the grade of membership is fully dened by either one or zero. Moreover, fuzzy numbers can be dened as convex and normalized fuzzy sets. Using fuzzy numbers, the parameters of the problem can be represented given that the underlying probability distribution is unknown or hard to be computed. For further details on Fuzzy variables the interested reader is referred to [91,[START_REF] Shapiro | Fuzzy random variables[END_REF].

Teodorovi¢ and Pavkovi¢ [START_REF] Teodorovi¢ | The fuzzy set theory approach to the vehicle routing problem when demand at nodes is uncertain[END_REF] consider the VRP under demands uncertainties. Using the fuzzy set theory, the authors model the demand at each customer as a triangular fuzzy number. Moreover, by using fuzzy arithmetic, it is shown that the remaining capacity after visiting a customer can also be represented by a triangular fuzzy number. Two approximate reasoning algorithms are then proposed.

The rst one, considers only the remaining capacity to decide if another customer should be serviced by a partial route. The second one, considers both the remaining capacity and the customer demand to decide if it should be included in a route or not. To solve the problem, the sweeping algorithm [70] is modied to consider the fuzzy logic. The solution is further evaluated within a simulation procedure to calculate the expected distance traveled due to failures. The classical recourse action is considered. Results are reported for a single instance with 100 customers with an improvement when the second reasoning algorithm is selected. A similar problem is considered by Cao and Lai [35] addressing the Open VRP with fuzzy demands (OVRPFD). The problem is modeled as a CCP, making use of fuzzy credibility measure to guarantee customers satisfaction. Demands are assumed to behave as a triangular fuzzy number and the decision maker can select a preference index representing is risk attitude towards the possible failures. Moreover, the classical recourse for the VRPSD is used when failures arise. The OVRPFD is solved with a hybrid method combining Monte Carlo simulation and a dierential evolution algorithm. Results are presented for instances with up to 100 customers.

Kuo et al. [START_REF] Kuo | Integration of fuzzy theory and ant algorithm for vehicle routing problem with time window[END_REF] worked on the VRP with time windows in a context where both travel times and 2.3. CONCLUSION time windows are modeled through fuzzy numbers. Both parameters make use of the triangular membership function. To solve the problem, a two phase algorithm is proposed. In the rst part, an Ant Colony optimization procedure is used to solve the associated TSP. Then in a second phase, routes are constructed using as input the results of the AC procedure. Routes are considered feasible if they respect the capacity and if the time windows respect a required service level. This one is considered as the largest membership value found in the intersection of the fuzzy arrival time and time window.

Some results are presented for a case study with 50 customers and dierent types of time windows.

He and Xu [81] deal with the VRP with uncertainties in customers demands and travel times. In this study, the authors consider travel times as independent normal random variables, and the demands as fuzzy variables. In fact, demands follow a normal distribution with a mean value represented by a trapezoidal fuzzy number. The variance of the demand is assumed to be crispy

Conclusion

The VRP is an active research eld. The last sixty years have shown a massive amount of research within the overwhelming dierent types of VRPs. Although most of the research has been concentrated on deterministic problems, in the past three decades an important amount of time was devoted to variants that consider the inherently uncertainty of information in many applications. How the uncertainties are handled is hardly dependent on the information available, e.g. if the underlying probability distribution of the parameters is known or not.

Stochastic optimization has been the main approach to solve uncertain VRPs. Even though SVRP eld is still at an early stage of development [63] it has an active community. Still some issues can be detected in the reviewed bibliography. For instance, many of the works deal with small to at most medium instances. That is the case for the VRPSD and the Cristiansen a Lysgaard [38] benchmark whose instances contain at most 60 customers. Incrementing the instances size of the testbeds can push knowledge boundaries by challenging the creation or adaptation of new methods (exact and heuristic)

to solve these problems. Moreover, this implies to design strategies that exploit characteristics of the SVRPs, while overcoming the inherent computational challenge of solving the SVRPs.

In SVRPs with stochastic times problems involving soft time windows (or mixed soft and hard time windows) are more common while the hard time windows case is very rare. Additionally, when hard time windows are considered, usually discrete probability distributions are used to handle the problems complexity. Thus, SVRP with hard time windows remains a fairly unexplored area in the eld. Besides, the largest amount of articles in the literature rely on convolution properties (also for the VRPSD) to deal with stochastic parameters. Even if convolution can arise naturally when dealing with random parameters, there is no proof that it is always (or recurrent) the case for SVRPs.

Therefore, solution approaches should include mechanism to be exible when convolution properties CHAPTER 2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS do not hold. This without forgetting the tradeo between accuracy and exibility.

Given the before, this rst part of the thesis focuses on solving the SVRPs for large instances as a way to create new comparable results. This is achieved by adapting solution methods with tailored strategies for the SVRPs, particularly for the VRPSD in chapter 3. Then a SVRP with stochastic travel and service times and hard time windows is considered in chapter 4.

Chapter 3

Hybrid metaheuristic for the VRPSD

Introduction

Since its introduction by [9], the Vehicle Routing Problem (VRP) has been widely studied in the literature, becoming a classical combinatorial problem. While the VRP comprise a broad family of problems, it is commonly used to refer to the Capacitated Vehicle Routing Problem (CVRP) version.

The CVRP objective is to build the set of vehicle routes with minimum cost satisfying customers demands. Stochastic Vehicle Routing Problems (SVRPs) are generalizations of the VRP where one or more parameters of the problem are associated with random variables. Recently, more attention has been paid to SVRPs. Their importance relies on their closeness to reality and their ability to take into account variability of data. A recent review of the main SVRPs variants studied in the literature can be found in [18].

To model SVRPs two stochastic approaches have been widely used in the literature: Stochastic

Programming with Recourse (SPR) and Chance Constraint Programming, using probabilistic constraints (CCP). Contrarily to CCP, the rst approach considers actions (called recourse) to overcome or react to possible violations of the constraints. In fact, since parameters are random variables sometimes constraints might not hold leading to failures. SPR takes into account the cost associated with recourse within the objective function. The CCP introduces constraints to limit the probability of failures to a threshold aiming to guarantee a quality level of the solution. It should be noted that SPRs and CCPs are not exclusive and can be used in mixed formulations (see for example [12]).

The Vehicle routing problem with stochastic demands (VRPSD) is an extension of the VRP in which the demand of each customer is a random variable. VRPSD was originally proposed by [40] who solve the problem through a modication of the well-known Clarke and Wright heuristic ( [8]). Exact methods can be found in the VRPSD literature and are divided in two approaches: the L-shaped algorithm (see [23]) and branch-and-price based. L-shaped methods ([?], [21], [24], [22], [4], [35]) have been the preferred approach to solve the VRPSD and are able to optimally solve instances with up to 100 customers and few vehicles when discrete distributions are considered ( [24]). In [22] normal distributions are used to model the demands attaining optimal solutions to instances with 60 to 80 nodes and two to four vehicles.

Branch-and-price based methods have shown to solve problems with larger number of vehicles. The rst branch-and-a-price algorithm for the VRPSD is presented in [7] to deal with instances with at most 60 customers under a Poisson demands assumption. More recently, a branch-and-cut-and-price algorithm is implemented by [15] and tested on instances with up to 101 customers and 15 vehicles.

Heuristics and metaheuristics methods are more often used to solve the VRPSD. In [39] a CCP and two SPR models are provided for the problem at hand, which is then solved by means of the Clarke and Wright heuristic [8] and a Lagrangian Relaxation based heuristic for instances with normally independent and correlated demands. In this study, a transformation of the VRPSD into CVRP

INTRODUCTION

under certain conditions is also discussed.

[?] propose a Tabu Search called tabustoch to tackle the extension of VRPSD where additionally, customers are present or not with a given probability.

[41] consider the VRPSD under restocking possibilities, thus adopting strategies for preventive restocking. That is, return trips to the depot to restock even if the vehicle is not empty to avoid future failures. The authors develop and embed an optimal restocking policy in the route design. In [3],

hybridization of ve metaheuristics with dierent objective function approximations and preventing restocking is devised. More recently, [20] propose a Simulated Annealing procedure to solve the VRPSD that uses a cyclic-order encoding to represent solutions. The authors designed a two-phase technique embedded in the solution method. In the rst phase they consider only deterministic costs, while the recourse cost is explicitly integrated into the second phase.

[28] use a Memetic Algorithm to solve the multi compartment VRPSD which is a generalization of the classic VRPSD. In the multi compartment VRPSD, the customers are associated to several products demands that cannot be mixed. This constraint imposes to load each product in a dierent compartment. In a more recent study, [29] designed a Greedy Randomized Adaptive Search Procedure (GRASP) enhanced with Heuristic Concentration (HC) to solve the VRPSD with maximum route duration constraints. To the best of our knowledge, the method of [29] reports the best overall results for [7] testbed.

Other variants of the VRPSD have also been considered in the literature. [27] addressed the VRPSD with weight-related costs and solved it by means of an Adaptive Large Neighborhood heuristic using several approximate methods. [37] dealt with the VRPSD with dynamic requests meaning that previously unknown customers can be received and scheduled over time. A Variable Neighborhood Search (VNS) based approach was proposed by the authors to solve both stochastic and dynamic cases.

[32] study the VRPSD with time windows using a Satiscing Measure Approach (SMA). The SMA is embedded in a tabu search showing very competitive results in small computational time. A multiobjective version of the VRPSD considering total traveling distance, total driver remuneration, number of vehicles and drivers remuneration balance is proposed in [16] and solved using a multi-objective evolutionary algorithm.

Concerning the stochastic models, SPR formulations have been dominant in the VRPSD literature compared with CCP formulations ( [38,39,10]). The most used recourse policy, called from now onward in this paper the classical recourse, is dened as follows. When the charge of the vehicle is emptied (fullled) it returns to the depot to replenish (unload) the charge, and resumes its assigned route from the failure point ( [2]). However, other recourse policies have been studied and implemented through several studies: preventive restocking policies are extensions of the classical recourse where return trips to the depot are performed even if the vehicle is not empty to avoid future failures ( [41,3,4,27,42,35,36]); pairing strategies allow the cooperation of multiple vehicles ( [1]); split deliveries between paired routes ( [26]) in which some customers are served by two vehicles; and backup routes ( [11]) that receive customers from primary routes.

Although the use of more complex recourse policies can represent a signicant saving relative to simpler ones ( [1]), the latter have been preferred since they allow more tractable models and stable tactical routes ( [17]). For this reason only the works of [4] and [35] deal with exact methods for the VRPSD using recourse actions dierent than the classical one. [4] consider a restocking policy for the generalized VRPSD using a single vehicle. In [35] optimal restocking policies allow vehicles to decide between a visit to the depot to replenish or proceeding to the next customer. Decisions are made using an optimal remaining capacity threshold for each customer within the route. The authors show that under arbitrary discrete probability distributions, instances with up to 60 customers and four vehicles can be solved. Moreover, it is shown that restocking policies can reduce the recourse cost to half of those achieved by the classical recourse. This chapter is dedicated to the vehicle routing problem with stochastic demands using the clas-CHAPTER 3. HYBRID METAHEURISTIC FOR THE VRPSD sical recourse. We propose a Memetic Algorithm hybridized with a GRASP (MA+GRASP) to solve eciently the VRPSD. The GRASP framework is embedded within a MA and is used as a way of restarting the algorithm. In this chapter, it is shown that the MA+GRASP is a valid and ecient method to solve the VRPSD. Moreover, a new testbed built from instances originally designed for the CVRP with up to 385 customers is proposed. The obtained results on new large instances serve as benchmark for testing future methods in real life scale problems. The remainder of this chapter is organized as follows. The problem formulation is introduced in section 3.2. Section 3.3 presents the developed solution approach for the problem. Numerical results are given and discussed in section 3.4. Finally, section 3.5 concludes the chapter.

Problem formulation

The traditional CVRP can be described as follows. Let G = (V, E) be a complete undirected graph where V = {0, 1, . . . , i, . . . , n} and E = {[i, j] ∀i, j ∈ V | i < j} are the vertex and the edge sets respectively. Moreover, V c = V \ {0} is the customers subset, each customer has a non-negative demand q i ∀i ∈ V c . Vertex 0 stands for a central depot where a homogeneous set of vehicles with a limited capacity Q each are initially located. Furthermore, each edge [i, j] ∈ E has a non-negative cost c ij . The objective is to build a set of routes with minimum cost considering that each route must start and end at the depot, the maximal capacity Q must be respected and, no split deliveries are allowed. This last constraint means that each customer is serviced once by a vehicle which deliver (or pickup) its whole demand.

The problem tackled in here presumes that the demand q i of each customer i follows a probability distribution ψ, with expected value and variance noted E [q i ] > 0 and V ar [q i ] > 0 respectively. We assume that probability function ψ is known and demands are mutually independent. Furthermore, we consider as other authors ( [7,15,20,29]), that ψ distribution has a cumulative property, i.e. the sum of demands probability functions is also ψ distributed. Furthermore, the vehicles eet is assumed to be unlimited and no xed cost per vehicle is involved. This problem is formulated in this work as an SPR in which the rst recourse policy of Bertsimas [2] is employed. This recourse (classical ) assumes that whenever a route achieves its maximum capacity Q, the vehicle returns to the depot to load/unload, then it returns to the customer where the capacity was fullled to complete the unserviced demand and then continues its route from the failure point. Let r be a route dened as a sequence of nodes r = {r 0 = 0, r 1 , . . . , r i , . . . , r k , r k+1 = 0}, the cumulative demand up to a client r i in a route r can be

dened as D ri = i j=1 q rj with E [D ri ] = i j=1 E [q r l ], V ar [D ri ] = i j=1 V ar [q r l ].
The expected recourse for a given customer r i ∈ V c in a route r is estimated by equation 4.1, as done by [15].

ERC ri = 2 • c 0ri • ∞ u=1 P D ri-1 ≤ uQ -P (D ri ≤ uQ) (3.1)
Given a route r, the term P D ri-1 ≤ uQ stands for the probability of cumulative demand up to customer r i being less than or equal to a multiple of the capacity of the vehicle. Therefore, the probability part of the expression represents the sum of having the u th failure of the route at client r i .

It shall be noticed that this expression does not consider the case when the remaining capacity equals the demand of a customer (exact stock-out). In such case the vehicle can return to the depot and then, continue towards the next customer in the route. Indeed, Hjorring and Holt [21] have already addressed the exact stock-out recourse. Nevertheless, we keep expression (4.1) since probability of such events might be rather low and its consideration can make the calculation computationally inecient.

Hence, the expected cost of a given route r can be calculated using equation 4.2.

E [C r ] = k j=0 c rj rj+1 + k j=1 ERC rj (3.2)

SOLUTION APPROACH: HYBRID METAHEURISTIC

In order to avoid multiple failures in a route r, the expected demand is limited to be at most equal to the maximum capacity Q, this assumption has been already used by many authors such as Laporte et al. [24], Christiansen and Lysgaard [7], and Gauvin et al. [15]. Indeed, if constraint (4.3) is not included, the optimal solution tends to be composed by only one route which visit all the customers.

k i=1 E [q ri ] ≤ Q ∀r (3.3)
The objective is thus to create a set of routes with minimum expected cost calculated by means of equation 4.2 and respecting constraint 4.3.

Solution approach: hybrid metaheuristic

A hybrid metaheuristic, combining a Memetic Algorithm (MA) and a GRASP, is proposed to solve large VRPSD instances. The proposed method is described in Algorithm 1. The MA basis is borrowed from the ideas of Prins [33] and it works with a xed population size. MA starts by creating an ordered initial population called P op (line 1 Algorithm 1). New individuals are created from the crossover of two chromosomes selected from P op (line 5 Algorithm 1). Moreover, a mutation procedure can be performed on the new individual with probability p mp (line 6 Algorithm 1). A local search is then executed with an associated probability of p ls (line 8 Algorithm 1) on the resulting solution. A procedure called Split is used to evaluate an individual tness (lines 7 and 10 Algorithm 1) and is presented in section 3.3.1. It allows also to convert chromosomes to VRPSD solutions by computing the detailed routes. Furthermore, when the local search is carried out, the routes are concatenated before using Split (line 10 Algorithm 1). To ensure a diversity on the population, clones are not allowed (line 12 Algorithm 1). Indeed, a new chromosome is kept for the next iteration only when its distance to the current population is not null. The distance measure used is the broken pairs [6] that counts the number of times a pair of consecutive customers in a rst individual is broken in a second one. When the new individual is accepted to enter the population, it is added to this last in a position that keeps the population ordered (line 13 Algorithm 1).

Since the MA works with a xed size population, when a new chromosome is entering the population, another one already in P op is removed. This last is randomly selected among those with tness superior to the median (line 13 Algorithm 1). Moreover, the MA uses a restart procedure. Indeed, after each φ iterations without improving the best solution the MA discards all the individuals except the best one (line 17 Algorithm 1). The population is completed using a GRASP procedure as explained in section 3.3.2. The algorithm stops when a time limit τ is achieved or if ρ iterations have been performed without improving the best solution found so far (line 4 Algorithm 1). end if 20: end while

Chromosomes

The MA + GRASP uses a twofold representation for each individual. The rst one follows the idea of Prins [33] and consists in representing a solution by a permutation of the V c customers. This representation have been already used by authors such as Mendoza et al. [28] for the multicompartment VRP with stochastic demands, Mendoza et al. [30] and Goodson et al. [20] for the VRPSD, and Mendoza et al. [29] for the VRPSD with time time duration constraints. The second one gives the detailed routes composing the solution. In order to decode the permutation of customers into a set of routes, the Split procedure presented in [33] is employed.

Split works by constructing a directed graph H = (W, Y ) composed by its vertex set W =

{W 0 = 0, W 1 , • • • , W i , • • • , W n }
, where W 0 serves as a dummy auxiliary vertex while the rest of

vertex {W 1 , • • • , W i , • • • , W n } ∀i ∈ V c are the permutation of the V c customers.
The arcs set Y is built in way that every arc (W i , W j ) | j > i represents a feasible route starting at the depot, visiting customers W i+1 , • • • , W j and returning to the depot. Since each arc is associated with a feasible route r, the arcs have an associated weight equal to E [C r ]. Therefore, it is during the construction of routes (arcs) composing graph H that stochasticity is considered, by properly calculating the costs using equation (4.2). After constructing the set of all feasible routes (respecting constraint 4.3) the goal is to nd the shortest path from vertex W 0 to W n and thus the arcs composing the shortest path are the optimal partition for the permutation of customers in the ordered sequence to routes. The set of routes is, in fact, the second representation of the individual, and is used when the local search is performed. For more details about Split method the reader is referred to [33]. shows the optimal decoding of the permutation into routes, with the cost of the whole set of routes which is the tness or cost of the individual. Concatenation procedure allows to pass from a set of routes representation to a permutation one.

It is achieved by removing the depot node from the start and end of the routes. Then, the customers sequences of the routes are added one after another. See for example the part (d) of gure 3.1, where routes 0 -2 -0, 0 -4 -1 -0, and 0 -5 -3 -0 derive to the customers permutation 2 -4 -1 -5 -3.

Besides, the calculus of the broken pair distance [6] is performed using the routes representation.

Initial population and Restart

Population P op is initially lled with three individuals created with the well-known heuristics: Clarke

and Wright [8], Gillet and Miller [19] and best insertion. Each of the aforementioned heuristics is executed twice with slightly dierences in the capacity Q. In the rst call to the heuristics, the full capacity of the vehicles (Q) is considered, while in the second one the capacity of vehicles is reduced

to Q = 0.9 • Q. That is, constraint (4.
3) right side is changed to 0.9 • Q. By doing so, it is expected that constructed routes will have a lower probability of failures (see [28]), introducing important information to the population. These heuristics are modied to consider the associated recourse costs.

The routes obtained by each heuristic are then concatenated to obtain a chromosome which is then evaluated by the splitting procedure described in 3.3.1. Only the best three individuals among the six generated ones are kept in the initial population, the other three are discarded. To complete the population P op, the remaining chromosomes are created from completely random permutations of customers which are next evaluated with the Split procedure. The population size is denoted P opsize afterward in the paper.

Restart

Restart procedure is executed after performing φ iterations without improving the best solution found so far. Algorithm 2 shows the pseudo-code of the procedure. It starts by deleting all the individuals of P op except the best one. In order to generate P opsize -1 missing individuals the next strategy based on a GRASP procedure is used. The GRASP generates

(P opsize-1) 2
individuals. Meanwhile, the remaining ones required to achieve P opsize are created from completely random permutations of customers, which are next evaluated with the Split procedure. This approach is used to ensure diversication within the restart.

Greedy Randomized Adaptive Procedure

Greedy Randomized Adaptive Search Procedure (GRASP) is a metaheuristic introduced by Feo and Rasende [13]. GRASP consists in using a randomized constructive procedure to build solutions that are improved after by a local search approach. This process is repeated through a number of iteration.

The behavior of the constructive procedure is controlled thanks to a parameter called greediness which permit to balance between greediness and randomness.

The proposed GRASP works as follows. Iteratively, the Greedy procedure of the GRASP, which is based on the Nearest Neighbor (RNN) heuristic is called. The RNN works by creating a permutation of the customers. At each iteration, the RNN picks randomly a client among the k nearest neighbors to the last visited customer (initially it starts from the depot), and add it to the permutation. After several preliminary tests we picked k with the expression M ax(2, |V | 60 ) and iterations to a value of 10. The RNN works without considering recourse costs as in [3,20]. After generating a permutation by using the RNN, it is decoded using the Split procedure. The resulting solution undergoes just after the local search procedure (see section 3.3.4) considering only deterministic costs to avoid the overhead of considering recourse costs. The best solution found by this steps is selected and kept.

Then, the local search procedure considering the recourse cost is called to improve the individual.

This last is added to the population and the process is repeated until

Crossover

Crossover procedure is performed in order to create new individuals, this is done by means of the OX crossover. Two individuals p 1 and p 2 from the population are selected from binary tournaments as well as two random positions i, j | i = j, j > i. Using the permutation representation of solutions, the information from position i to position j (included) are copied from p 1 to the new individual in the same positions. Additionally, p 2 is circularly traversed from position j + 1 to j completing the ospring from j + 1 to i, with clients in p 2 not included yet. By changing the roles of p 1 and p 2 , another individual can be created using the same procedure. Among the two new chromosomes one is randomly retained.

Mutation and Local Search

After being created, a new chromosome can be modied by the mutation procedure (line 6 Algorithm 1). It consists in moving nm random selected customers from their current positions to new ones randomly selected. In our implementation, after preliminary tests the value of nm is set to two. This operation is performed on the chromosome. The proposed local search (LS) for the MA is organized as a Variable Neighborhood Descent (VND), a variant of the Variable Neighborhood Search proposed in [31]. The LS is performed on the routes representation. The neighborhoods used for the VND (line 8 Algorithm 1) are the Or-opt and 2-opt movements in their intra and inter route versions, and the inter-CROSS movement. The order in which neighborhoods are explored is randomly selected each time the LS procedure is called. For a further related review on neighborhoods structures and types, the reader is referred to [5]. The LS procedure starts by exploring the neighborhoods trying to improve a solution. Whenever an improving movement is found in a current neighborhood it is executed, and the procedure passes to the next one when no improving movement is possible. Each time a solution has been improved, the procedure restarts from the rst neighborhood, running so on until any of them is able to enhance the solution. Inter route Or-opt movements are limited to sequences of at most three customers, and inter-cross movements exchange to at most two costumers

Numerical results

The tests are carried on two groups of instances: the rst group is already used in the VRPSD literature and is due to Christiansen and Lysgaard [7]. The second one, with larger graphs, is elaborated in this work to asses the performance of the developed approach on real size benchmarks which can be used for future works. The detailed results are provided in the following sections.

Classical testbed from Christiansen and Lysgaard

In order to assess the performance of the hybrid MA+GRASP, the tests are carried on 40 existing instances proposed in [7]. In this benchmark, the number of customers varies from 16 to 60 and the minimum number of vehicles needed to satisfy the customers demands is comprised between two and fteen. Christiansen and Lysgaard [7] testbed is based on Augerat test sets A and P and Christodes and Elion test set E. Demands are assumed to be Poisson distributed as in [7,20,30,15,29], with expected values equal to the deterministic demand values. Moreover, travel costs are calculated as the Euclidean distance between two nodes rounded to the nearest integer as done in [15]. For Christiansen and Lysgaard [7] benchmark, the best Known Solutions (BKS) are either taken from [15] in which 38 solutions are proven to be optimal, or from [29] and [20] which use heuristic approaches. All tests were conducted on a Dell Latitude E6420 personal computer with Intel Core i7-2760QM 2.4 GHz, running under Windows 7 Professional 64 bits. The algorithms were coded on Java and compiled with JavaSE-1.8 4 5 with maximum allocated memory of 1 Gb. Random variables and computation probabilities were generated by the library of Stochastic Simulation in Java ( [25]). Preliminary tests were performed to select the parameters used by the heuristic, local search rate and mutation probability. These two last parameters are set to 0.15 and 0.2 respectively. The hybrid MA+GRASP stops after 5000 iterations without improving the best solution or after running for 10 seconds.

Table 4.1 summarizes the results of the proposed MA+GRASP compared to those of Mendoza et al. [29] (GRASP-HC), and Goodson et al. [20] obtained with a Simulated Annealing (SA) procedure.

For each method are reported: the average gap on 10 runs (Avg. Gap), the maximum average gap across the 40 instances (Max Gap), the average gap of the best solution found over the 10 runs (Avg.

Best Gap), the number of best known solutions found considering the 10 runs (NBKS), the average time over 10 runs for the 40 instances (Avg. CPU), and the maximal and minimal average time over the 40 instances (Max. CPU Min. CPU).

The MA+GRASP and the GRASP-HC achieve to nd all the BKS (40 out of 40) whereas SA nd 33 BKS. Moreover, our method reaches an average gap below 0.01% (0.004%) which is nearly four times lower than the GRASP-HC. The low average gap also shows the stability of the method when dealing with dierent types of instances. The Max. Gap of 0.14% conrms the ability of the MA+GRASP method to regularly nd near optimal solutions, this metric shows also that our results addressed by Gauvin et al. [15]. The proposed testbed has on average nearly 2.4 more customers than the benchmark of [7] which oers the opportunity for future comparisons of new methods dealing with the VRPSD. Moreover, travel costs are calculated as the Euclidean distance between two nodes rounded to the nearest integer as done in [15] Some adjustments were made for the stopping condition of our hybrid method, for instances with Min veh•Q , the best known solution (BKS) provided either by [15] or among the several preliminary runs and the reported results, the expected cost of the optimal deterministic solution (BDS) 1 in the presence of uncertain demands 2 , and the value of the stochastic solution (VSS) which stands for the percentage of improvement among the BDS and the BKS. Furthermore, for the BKS is presented: the number of vehicles in the BKS (Veh), the total expected cost (Total), the deterministic cost (Det), and the recourse cost (Rec). The BDS also presents its total expected cost (Total) and its recourse cost (Rec).

As shown in table 3.3 most of the instances (38 out 39) show a positive VSS. Instance E-n23-k3 present a null VSS since the optimal deterministic solution is also optimal in the stochastic scenario.

Overall the VSS rounds the 5.37%, showing the importance of considering the stochasticity. VSS seems to grow with the number of nodes within each instance, for instances with less than 100 nodes its average value is 4.67%. Meanwhile, this value increases to 6.39% for instances with more than 100 nodes. As well, instances with more than 150 nodes achieve a VSS of 7.93%. Even if these results are not conclusive in a direct relation between the VSS and the number of nodes, it seems that the larger are the instances, the higher is the VSS. Additionally, the BKS cost (Total) is mainly composed by the deterministic cost averaging 97.38%, while the recourse only achieves 2.66%. Indeed, the positive VSS can be explained as follows. The BDS is optimal in the deterministic component of the cost. However, the recourse value of BDS solutions becomes far more important than it is in the BKS. Indeed, the recourse cost averages 10.88% of the Total BDS in contrast of the 2.55% of the Total BKS. Therefore, the trade-o between the deterministic and recourse cost is better in the BKS, generating a positive VSS.

Tables 3.4 and 3.5 illustrate the results on the new set of instances using the MA+GRASP and two alternative versions, namely MA+RANDOM and NR-MA. MA+RANDOM only changes the way new individuals are created during the restart procedure, using completely random customers permutations while NR-MA is a simpler version of the algorithm without considering the restart procedure. For each instance on this testbed are provided in table 3.4: the best known solution (BKS) with the same considerations as in table 3.3, the best cost found out on the 10 runs (Best Cost) and the average cost 3.5 are reported: the computational time (Time) in seconds, the gap between the average cost and the BKS (Avg. Gap) and the gap between the best solution found and the BKS (Best Gap). Moreover, ten out of the 39 instances have a proven optimal solution given by Gauvin et al. [15] 3 , these values are marked with an asterisk in table 3.4.

Overall, the MA+GRASP shows a good performance presenting 33 out of 39 instances with an average gap bellow one percent. Furthermore, the average gap is as low as 0.380%. MA-RANDOM ranks second in this metric with 0.418% and the NR-MA achieves 0.434%. MA+GRASP presents an average gap of 0.64% for instances with 100 customers or more, while instances with fewer customers show an average gap of 0.20%. The dierence can be explained by the inherent combinatorial nature of the problems at hand. MA+GRASP accomplishes to nd nine out of ten proven optimal solutions, still instance P-n65-k10 is the only one which cannot be found by this method. Furthermore, MA+GRASP

General performance discussion

Authors like Laporte et al. [24] or Jabali et al. [22] had reported the diculty of their exact approaches to solve the VRPSD when the number of vehicles increases or when the lling coecient (FC) approaches one. Although in our approach the eet is considered unlimited instead of xed, we test if such parameters had an impact on our solution method. We use as proxy the average gap retrieved by the MA+GRASP to compare it against both the FC and the number of vehicles. Figure 3.3 presents the relation between the FC and the average gap. In the left gure the FC reported in table 3.3 is used while the right part represents a FC calculated using the actual number of vehicles used in the BKS. Instance E-n23-k3 was discarded for the graphs as far as this point is a very extreme point 4 . The graph with the original FC shows a slightly positive trend between the variables while the modied FC presents no relation at all.

The same approach to FC is used to compare the avg gap of MA+GRASP with the number of vehicles used in the BKS.

Addressing the eect of GRASP restart

In order to test the impact of the GRASP and the restart procedure included in the MA, the variance analysis method of Friedman is used to compare the MA+GRASP against the MA+RANDOM. The results presented in table 3.4 are intended to perform two comparisons in terms of the best solution found, and the average cost. The Friedman tests works as follows. For each methods comparison, results are ranked within each instance giving one to the best value to two to the worst one, average values are used in case of ties (table 3.7). Let R (X ij ) be the rank of the heuristic j in instance i, b the number of instances, k the number of tested heuristics and

R j = b i=1 R (X ij ) ∀j). Moreover, A = b i=1 k j=1 (R (X ij )) 2 , C = bk(k+1) 2 4 , T 1 = (k-1) k j=1 (Rj - b(k+1) 2 ) 2 A-C , T 2 = (b-1)T 1 b(k+1)-T 1 . In addition t 1-( α 2 ) stands for the
|R l -R m | > t 1-( α 2 ) 2 bA - k j=1 R j 2 (b -1) (k -1) (3.5)
Let R 1 , R 2 be the R j metrics for the MA+GRASP, and MA+RANDOM respectively. Table 4.5

summarizes the necessary values for performing the tow two-pair comparisons. Moreover, the right side of equation 4.5 has values of 14.33 and 7.77 provided that α = 1% for the average cost and α = 5% for the best cost metric. As far as equation 4.5 holds for every |R l -R m | we can conclude the statistical dierence the heuristics within each metric and therefore, MA+GRASP can be determined as the best alternative. Indeed, results show that the MA+GRASP presents the best performance and it is statistically signicant on the average cost at a condence level of α = 1%, and for the best solution found at The mttt-plot method is an extension of the time-to-target plots (ttt-plots) introduced in [14]. According to Reyes and Ribeiro [34] the plots can be used to compare the running times of stochastic algorithms or dierent strategies for solving a given problem. The approach is based on the construction by simulation presented in [34]. The resulting graphs represents the probability of nding the proposed target values (one per instance) for the set of instances in a specic amount of time.

A-n61-k9 1 1 2 A-n62-k8 1 1.5 1.5 A-n63-k10 1 1.5 1.5 A-n63-k9 1 1 2 A-n64-k9 1.5 1.5 1.5 1.5 A-n65-k9 1 1.5 1.5 A-n69-k9 1 1.5 1.5 A-n80-k10 1 1.5 1.5 CMT12 1.5 1.5 1.5 1.5 CMT11 2 1 2 CMT4 1 1 2 CMT5 2 2 1 E-n101-k14 1 1.5 1.5 E-n101-k8 2 1.5 1.5 E-n23-k3 1.5 1.5 1.5 1.5 E-n30-k3 1.5 1.5 1.5 1.5 E-n76-k10 1 1.5 1.5 E-n76-k14 1 1.5 1.5 E-n76-k7 2 1.5 1.5 E-n76-k8 2 1.5 1.5 P-n101-k4 1 1 2 P-n55-k8 1.5 1.5 1.5 1.5 P-n65-k10 1 1 2 P-n70-k10 1 1.5 1.5 P-n76-k4 2 1 2 P-n76-k5 1 1.5 1.5 Tai100a 1 1 2 Tai100b 2 1.5 1.5 Tai100c 1 1 2 Tai100d 2 1.5 1.5 Tai150a 1 2 1 Tai150b 2 2 1 Tai150c 1 2 1 Tai150d 1 1 2 Tai75a 1 1 2
To extract the necessary information, we run again the three methods for each of the new testbed instances. The stopping criteria is set to a maximum time of half an hour or when a BKS is achieved.

During the execution time, for each instance run and method, the best solution is always surveyed, so when a new best solution is found the necessary time to nd it is recorded. Ten runs are performed per instance and per method. Using this information we estimate the probability distribution of the time required to nd dierent target values. Indeed, we vary the target value for each instance as a function of a gap to the BKS reported in table 3.3 6 . This allows to see how the methods perform over dierent targets giving a more comprehensive story of the performance of the algorithms.

As far as we imposed a limit time, it may happen that a method cannot meet the required target.

In such case we let the maximum time as record. The rationale behind this is that even if we do not know the actual time required to nd the target value, it stills contains valuable information.

Indeed, more stable methods should be able to nd the BKS (or very near BKS) targets more often.

Therefore, letting the maximum time as record when the target value is not found accomplishes the aim of characterizing the solution method times. Moreover, this allows to limit the computational times required to perform a mttt-plot analysis, as far as some methods can require huge amount of time to nd certain targets for some instances.

Figure 3.6 shows the mttt-plots for target values not farther than 5% and 1% of the BKS. The 5% gure concludes that the NR-MA is better than MA+RANDOM and MA+GRASP, and with less time it has higher probabilities of achieving a 5% gap to BKS for the whole benchmark of instances.

Nevertheless, the 1% gure displays a completely dierent behavior. In this last, the MA+RANDOM and MA+GRASP present a very similar behavior which is far better than NR-MA. For example, the probability of nding a target value within 1% of the BKS (for the whole set of instances) within 2500 seconds is one for MA+GRASP and MA+RANDOM, for NR-MA is around 70%.

Going further to determine the best method, gure 3.7 presents the mttt-plots for target values of 0.5% and 0% of the BKS. In both targets it is shown that MA+GRASP outperforms the MA+RANDOM and the NR-MA showing a higher probability of nding the matches for a xed amount of time. It is interesting that as the target approaches the BKS the performance of MA+RANDOM and NR-MA become very similar. This shows that the random restart might have a very little impact on the solution approach in long runs. However, the restart based on the GRASP does show an important impact, giving the edge to MA+GRASP as the best option among the tested ones to solve the VRPSD.

Conclusions

In this chapter the Vehicle Routing Problem with Stochastic Demands (VRPSD) was studied. In order to solve this problem a MA+GRASP method is proposed. The obtained results on a classical testbed from Christiansen and Lysgaard [7] show that our method outperforms state-of-the-art algorithms in This chapter deals with maintenance scheduling and routing problems where groups of technicians must be assigned to visit a set of customers to execute repairing tasks within given time windows.

The goal is to build a minimum cost set of routes subject to the following constraints: every route starts and ends at the depot; each customer is visited and serviced once by only one route within its time window, and the sum of customer demands on each route does not exceed the vehicle capacity.

With regard to time window constraints, there are two cases that are encountered in the literature:

hard time windows exist when the service cannot start outside the time interval; early arrivals are permitted but the vehicle must wait until the opening of the time window. In contrast, soft time windows services are allowed outside the time windows usually with a penalty cost. In this study, the is considered. In fact, travel times modeled as random variables consider the uncertainty in time due to trac, weather, accidents, driving skills, etc. Furthermore, the stochastic service time accounts for the complexity of maintenance tasks, which may only be known when technicians arrive at customer locations.

A recent survey on SVRPs is presented in Gendreau et al. [16] reviewing dierent variants, e.g. stochastic demands, customers, travel and service times. The version with stochastic demands (VRPSD) was rst proposed by Tillman [45] with later works on exact methods in [28,8,15] and heuristic methods in [42,19,33]. A VRPSD where the presence of customers is also a random variable is presented by Gendreau et al. [17] and solved with a Tabu Search.

Earlier related works concerning stochastic travel and/or times problems include those of Laporte et al. [27] for the Vehicle Routing Problem with Stochastic Travel and Service times, in which the authors propose three models, including a CCP and an SPR, solving instances with up to 20 nodes where travel times are restricted to ve discrete states. Lambert et al. [26] propose an adaptation of the Clarke and Wright [9] saving heuristic to solve the VRP with stochastic travel times for money collection. Later, two models for the VRP with stochastic travel and service times are presented in Kenyon and Morton [24]. The rst minimizes the expected completion time, while the second model maximizes the probability of completion within a pre-specied deadline. To solve small instances, the authors used a branch-and-cut algorithm; for larger ones, the branch-and-cut scheme was embedded in a Monte Carlo sampling approach.

A Genetic Algorithm reporting signicant improvements in terms of costs and delay penalties is referred to in [2] for the SVRP with time windows and stochastic travel times. More recently, the SVRP with stochastic travel and service times is modeled by Li et al. [31] as a CCP and an SPR.

The problem is solved using a Tabu Search, and Monte Carlo simulation is used to check stochastic constraint feasibility and estimate the expected recourse value. The authors test this approach on their own instances with up to 100 clients and soft time windows, concluding that CCP models are harder to solve than SPR. Zhang et al. [49] also address the SVRP with stochastic travel and service times with soft time windows. The authors propose a CCP model to guarantee a service level to customers services, arrival at the depot and on the total duration of the routes. To estimate the vehicles arrival times, a discrete approximation method is embedded into a Tabu search heuristic to solve some Solomon [41] instances with up to 20 customers. A Tabu Search is also proposed in Ta³ et al. [43] for the SVRP with stochastic travel times and soft time windows. The objective is to minimize a weighted cost composed of two parts: the rst relates to the service (early and late arrivals) and the second to transportation (distance, xed vehicle costs and overtime). Ta³ et al. [44] analyze the same problem but deal with stochastic time-dependent travel times. The authors deduced a way of calculating exactly the mean and variance of arrival times provided that service times for customers equal zero. Conversely, when service times are non-zero, an approximation is used. Both, exact and approximation schemes are derived if travel times are Gamma distributed.

A Satiscing Measure Approach (SMA) to mitigate the dissatisfaction experienced by customers is proposed in Nguyen et al. [39]. The SMA is realized for both the SVRP with stochastic demands and time windows and the SVRP with stochastic travel times and time windows. The latter model considers the early time window (ready time) as a strict requirement and the late time window (due time) as a soft requirement. The results for Solomon [41] instances show small computational times.

Errico et al. [13] eectively solve the SVRP with hard time windows and stochastic service times with a branch-cut-and-price algorithm for instances with up to 50 customers. The model proposed by the authors is an SPR, nevertheless two conditions (constraints) are imposed: a service level on being operationally feasible for each route is required and a maximum of one recourse per route is allowed. A variant of the problem that considers multiple depots and client priorities is addressed by Binart et al. [4]. The problem incorporates two types of customers: optional and mandatory. The TIMES former customers have an associated hard time window that must be respected. The problem is solved with a two-stage method and is proved eective in instances containing up to 50 customers and three vehicles.

Jula et al. [22] develop approximations for the mean and variance of the arrival times for the Traveling Salesman Problem with time windows (TSPTW) which is solved via dynamic programming. The arrival times mean and variance are estimated using a rst-order Taylors series expansion.

The authors guarantee a service level requirement on clients time windows using the Chebyshev and Cherno bounds. In the same vein, Ehmke et al. [11] estimate both the arrival and service times for the VRP with stochastic travel times and time windows using a normal approximation. The authors embed their approximation within a Tabu Search method to solve a CCP formulation on Solomon [41] instances where late services are allowed. The work of Gomez et al. [18] population to enhance its overall performance. Moreover, the contributions found in this chapter are: it is shown that despite the eects of hard time windows and the sum of random variables with dierent probability distributions, the proposed approximation of the arrival times is a valid and fast approach to guide a solution algorithm. Moreover, the experiments made with the MPMA using dierent populations reveal that tackling the problem with dierent assumptions on the parameters, while the populations share their solutions, greatly improves the best and average solutions found.

Furthermore, the use of a combined model (CCP + SPR) allows it to easily adapt the problem to comprise dierent objectives, such as costs and customers satisfaction.

This chapter is organized as follows. The SVRP with hard time windows and stochastic service and travel times is introduced in section 4.2. The approach to estimate arrival times is assessed in section 4.3. In section 4.4, the Multi-Population Memetic Algorithm developed to solve the problem is described. Numerical results are reported in section 4.5, and lastly a conclusion is put forward in section 4.6.

Problem Denition

The Vehicle Routing Problem with hard Time Windows and Stochastic Travel and Service Times (SVRPTW) is a generalization of the VRPTW where travel and service times are modeled by random variables. The VRPTW is dened by a vertex set V = {0, 1, . . . , i, . . . , n} and an edge set

E = {[i, j]∀i, j ∈ V | i < j} that composes of a complete graph G = (V, E). Each vertex i ∈ V is
characterized by a coordinate (x i , y i ), and a time window [e i , l i ] in which the service must start; e i being the opening time and l i the closure time. Furthermore, V c = V\ {0} is the customers subset, each of which has a non-negative demand q i and a specied service time s i . Vertex 0 represents a depot where a eet of homogenous vehicles with a limited capacity Q is located. In addition, to each edge [i, j] ∈ E there is an associated non-negative cost c ij and a travel time t ij .

Let M be the xed cost associated to each used vehicle, and the vertex n + 1 a dummy copy 4.2. PROBLEM DEFINITION of the depot 0 with the same location and time window, moreover allow

V = V ∪ {n + 1} and E = E ∪ {[i, n + 1] ∀i ∈ V c } with V c = V\ {0}.
Service times and travel times are random variables denoted by s i (∀i ∈ V c ) and t ij (∀ [i, j] ∈ E ) respectively.

To model the SVRPTW, a combined CCP and SPR formulation is used to deal with dierent settings. On the one hand, the CCP part lets managerial decisions to be considered when solving the problem (by controlling the condence levels). Conversely, the SPR adds the recourse cost component, hence ensuring a more reliable measure of the quality of the solution. Our model integrates stochastic constraints (CCP) to guarantee the condence levels α, β and γ for hard time windows constraints at clients, depot time window, and success for the whole set of routes respectively. For simplicity, it is assumed that every client in V c has the same required condence level α. Let r be a route dened as an ordered sequence of clients r = {r 0 = 0, r 1 , . . . , r j , . . . , r k , r k+1 = n + 1} where r j represents the j th visited customer. Furthermore, AT rj stands for the arrival time to client r j ∈ V c on the route r and P (A) the probability of an event A. AT rj is a random variable because it depends on travel and service times which are dened as random variables. Aiming to guarantee the condence levels to the above-mentioned, the following constraints are imposed on every route r:

P AT rj ≤ l rj ≥ α ∀r j ∈ V c (4.1) P AT r k+1 ≤ l 0 ≥ β (4.2)
Constraint (4.1) enforces a service level for every customer, stating that the customer service must start within its time window a with probability of at least α. Also by (4.2) it is guaranteed that vehicles will return within the depot time window with a probability of at least β. Even when a route r meets equations (4.1) and (4.2), sometimes it can miss the customer or depot time window closure;

these events are called failures. Let U r be the probability of having no failures in a route r. A solution for the SVRPTW will normally be composed of more than one route, therefore let s be a solution for the SVRPTW composed by a set K of routes; the following constraint is imposed for every solution:

r∈s U r ≥ γ (4.3) Equation (4.
3) guarantees a service level on the whole route plan rather than only on customers [12], i.e. none of the routes in the solution will miss time a window constraint with customers or at the depot, with a probability γ. This formulation is valid if, and only if, the set of routes in the solution are independent. This assumption holds if travel times are independent and also if every client is visited by only one vehicle, i.e. the probability of no failures is not related from one route to others.

The recourse action is considered as follows: if a vehicle arrives at a customer i ∈ V c later than the closure of its time window, the vehicle will continue its route (without performing the service) towards the next client. The recourse is founded in the idea of rescheduling the visit of customer i later. The associated recourse cost can be seen as a penalization for missing the customer time window, this penalty is represented by the cost of a vehicle visiting exclusively customer i. The recourse albeit simple, ts the maintenance scheduling problem. As far as service and travel times can only be known until they are completed, technicians cannot anticipate failures. Therefore, a failure is only known at the customer location where it takes place. Furthermore, this simple recourse lets us introduce the cost of missing a service while keeping the computations of probabilities tractable.

The proposed recourse has already been addressed in the literature: Nguyen et al. [39] described this recourse action, nevertheless, it is not used in their approach; Wang and Regan [48] also proposed this recourse but the authors did not associate a cost to it. In this work we dene the expected cost of a route r as:

c r = M + k j=0 c rj rj+1 + k j=1 P AT rj > l rj • 2 • c 0,rj + M (4.4) CHAPTER 4.

VEHICLE ROUTING PROBLEM WITH STOCHASTIC TRAVEL AND SERVICE TIMES

The rst part of (4.4) stands for the xed vehicle costs M ; the second part is the cost of edges traversed by the route and the third represents the expected recourse cost. It is important to note that the recourse cost does not depend on the length of the delay for the late arrival as it only uses the probability of missing the time window.

In here, the eet size is considered unlimited and homogenous. The probability density functions ψ i (∀i ∈ V c ) associated with service times are known, and service and travel times are assumed to be mutually independent. The probability density function φ ij ∀ [i, j] ∈ E of every travel time is known.

We set the capacity of the vehicles to innity (Q = ∞) considering that it does not limit the technicians capacity to provide their services. We assume that service times are identically Gamma distributed and travel times are identically distributed with a Log-normal distribution. Kaparias et al. [23] and Lecluyse et al. [29] have already used log-normal distributions, recognizing the importance of skewed distributions to model travel times. Furthermore, Gamma distribution is selected for the service times as long as it respects the principle of increasing repair rate.

Estimation of arrival times 4.3.1 Arrival and starting service times denition

Since the vehicles arrival times to customers and depot depends on travel and service times which are by denition random variables, thus arrival times become random variables too. The same condition applies to the starting time when the service is performed at a customer, as long as it depends on the arrival times and time windows. Let ST rj denote the random starting time of the service at client r j ∈ V c in a route r. Because of the problems dened, the service is performed only if AT rj ≤ l rj .

Otherwise, if a failure takes place ( AT rj > l rj ) the vehicle continues towards the next node in its route, i.e. the service is not performed. Also, because time windows are hard, the service for a customer can only start at or after the opening of its time window. Let 1 Z be an indicator function which takes value one if a condition Z holds or zero otherwise, the time when a service starts at a customer can thus be set as:

ST rj = e rj • 1 { AT r j <er j } + AT rj • 1 {er j ≤ AT r j ≤lr j } (4.5)

Aiming to dene arrival times and considering the used recourse, it is mandatory to check if a failure took place at the last customer visited or not. Therefore, arrival times can be dened by equation (4.6).

It should be noted that recursively using equation (4.6) implies that failures aect the distribution of the arrival and initial service times.

AT rj = ST rj-1 + s rj-1 + t rj-1rj , AT rj-1 ≤ l rj-1 AT rj-1 + t rj-1rj , AT rj-1 > l rj-1 (4.6)

Mean and Variance estimation

To estimate the mean and variance of the arrival times an approach similar to the one used by Ehmke et al. [11] is proposed, nevertheless, we do take into account the impact of possible service failures on AT rj as well as the fact that service times are random.

Let µ X = E X be the expected value of a variable X and σ 2 X = V ar X its variance. By denition the standard deviation is set to σ X = σ 2 X . For simplicity let P 1 = P e rj ≤ AT rj ≤ l rj | AT rj ≤ l rj and P 2 = 1 -P 1 = P AT rj < e rj | AT rj ≤ l rj . Applying the laws of total expectation and total variance to the equation (4.5), the mean and variance for the initial service times at customer j are presented in equations (4.8) and (4.9). µ ST r j | AT r j ≤lr j = µ AT r j |er j ≤ AT r j ≤lr j

• P 1 + e rj • P 2 (4.8) σ 2 ST r j | AT r j ≤lr j = σ 2
AT r j |er j ≤ AT r j ≤lr j

• P 1 + e 2 rj • P 1 • P 2 +µ 2
AT r j |er j ≤ AT r j ≤lr j

• P 1 • P 2 -2 • e rj • µ AT r j |er j ≤ AT r j ≤lr j • P 1 • P 2 (4.9)
Using the results from equations (4.8) and (4.9), we now apply the laws of total expectation and total variance to equation (4.6). Again for simplicity let P 3 = P AT rj-1 ≤ l rj-1 and P 4 = 1 -P 3 = P AT rj-1 > l rj-1 . The mean and variance for the arrival times at a node j are dened by equations (4.10) and (4.11).

µ AT r j = µ ST r j-1 | AT r j-1 ≤lr j-1 + µ sr j-1 + µ tr j-1 r j • P 3 + µ AT r j-1 | AT r j-1 >lr j-1 + µ tr j-1 r j • P 4 (4.10) σ 2 AT r j = σ 2 ST r j-1 | AT r j-1 ≤lr j-1 + σ 2 sr j-1 + σ 2 tr j-1 r j • P 3 + σ 2 AT r j-1 | AT r j-1 >lr j-1 + σ 2 tr j-1 r j • P 4 + (P 3 • P 4 ) • µ ST r j-1 | AT r j-1 ≤lr j-1 + µ sr j-1 -µ AT r j-1 | AT r j-1 >lr j-1 2 (4.11)
Although the equations (4.8) to (4.11) are computed for AT rj , they are valid for AT rj too, since no assumption is made on the distribution. Moreover, by using AT rj instead of AT rj two problems are solved. The rst one is that the calculus of probabilities can be easily made. The second is that the calculation of the mean and variance of the truncated variables, e.g. the mean of the upper truncated arrival time (µ AT r j | AT r j ≤lr j ), can be evaluated, for example using the closed forms expressions gathered in [21] for several distributions.

Summing up, by iteratively applying equations (4.8) to (4.11) using AT rj , one can calculate the parameters of the AT rj distribution for a given route. Assuming a technician starts its route from the depot at time zero with no variance, and using equations (4.8) and (4.9), the mean and variance of the starting service time at the rst customer are derived. This can be easily achieved since AT r1 depends only on the travel time from the depot to r 1 . Next, the parameters of the arrival time at the second customer can be calculated using equations (4.10) and (4.11), all the more, after performing this step, equations (4.8) and (4.9) are used to compute the parameters of the starting service time at the second customer. The process is repeated until the whole route is evaluated. Additionally, by replacing AT rj in constraints (4.1) and (4.2), the feasibility of the route can be checked.

Multiple authors have used dierent distributions and assumptions to estimate the arrival times (readers can refer to [7,11,18,34]). It is assumed that AT rj follows a log-normal probability distribution. The log-normal assumption is twofold motivated. First, an experiment using Monte Carlo simulation shows that statistically, the log-normal distribution tted the arrival times better than other distributions (e.g. normal, Gamma). Second, although normality assumption has been proven to be eective in previous works [7,11]; skewness is an important factor to be considered in stochastic vehicle routing algorithms to lead to reliable routing decisions [18]. The presence of left truncation due to early arrivals induces more asymmetry in the time when service starts which is transferred to arrival times, and because normal distribution has zero-skewness it does not appear to be the best distribution to approximate arrival times. Based on these reasons and the tractable computational times, the log-normal assumption is retained.

Validating the log-normality approximation

An experiment to validate the log-normality assumption is conducted as follows. A solution is rst created for each instance (see section 4.5.1 for further detail on the instances) using the Clarke and Wright [9] heuristic combined with the simulation procedure to verify constraints (4.1, and 4.2). For each tested instance, we considered the two routes in the solution with the greatest number of visited clients. The reason for this choice is that larger routes will be harder to t and to approximate as far as more truncation eects are summed up. Furthermore, for the whole set of selected routes a test is performed to compare the mean, standard deviation, and percentiles of the arrival times at each node. This is estimated through simulation (10000 trials) against the values obtained while considering the variables AT i log-normally distributed (estimated). The reported gaps are calculated as

|X simulation -X estimated | X simulation
where X simulation and X estimated are replaced by the mean, standard deviation and each percentile of the arrival times. Table 4.1 shows basic information of the experiments, including the number of evaluated routes, the total number of evaluated nodes, the average number of nodes per route and the percentage of arrival times that tted log-normal, gamma and normal distributions.

Although more classic probability distributions were tested e.g. exponential, Weibull, chi-squared, Poisson, Pareto, triangular, uniform, Cauchy, logistic, Laplace, and Erlang, we reported the results respectively. The standard Kolmogorov-Smirnov test was used for testing the distributions. The rst customer of each route was not considered because the arrival time distribution will be the same distribution of travel times between the depot and the customer. The test concludes that nearly one in three of the arrival times are log-normal distributed while this number decrease to around 15% when tting a normal distribution. Gamma distribution achieves to t around 28% of arrival times. Based on these results it appears that the log-normal approximation can be used as a compromise solution to model the arrival times while allowing tractable computation of the cost function and probabilities.

ESTIMATION OF ARRIVAL TIMES

Thus, we selected it to model AT rj and continued testing its pertinence. As shown in Table 4.2, the proposed approximation displays an average gap of the arrival times mean of at most one tenth of one percent, showing good performance. The lowest average gap is reached by families R2 and RC2

where the mean has an average gap of only 0.03% while the highest gap of 0.1% is achieved in family C1. On this, standard deviations (SD) have bigger gaps, yet, half of the studied families present an average standard deviation gap below one percent, and only one family of instances exceeds a gap of two percent.

Given that equations (4.1) and (4.2) make use of probabilities, it is mandatory to test how accurate the estimations of dierent percentiles are. Therefore, table 4.2 also presents the average absolute gap for some notable percentiles and the mean, minimum, and maximum, average absolute gap when the 1 st to 99 th percentiles are evaluated with increments of one percent.

Results in table 4.2 show the pertinence of the log-normality assumption. The 90 th percentile has an average gap under a quarter of one percent, while the 95 th percentile estimation is on average quite bellow one percent error rate for the whole set of instances. Moreover, the 99 th percentile gives an average gap less than one percent for all families except for RC1 for which the maximum is reached at 1.45%. It should be noted that the errors grow with the percentile value. The proposed approach is very exible and can easily modify the assumption that AT rj is lognormally distributed to other probability distributions, e.g. Normal, Gamma, Exponential. In fact, the mean and standard deviation of AT rj were also calculated assuming AT rj is normally and Gamma distributed. However, log-normal approach presents the best results, improving by more than 80% of the absolute gap error of arrival times' mean when compared to the normality assumption (see table 4.3). Indeed, these results with those reported in table 4.2 show that the log-normal distribution gives better results in terms of parameters and percentiles estimations. Therefore, we retained this distribution as the best option. To check the robustness of the log-normal approximation we conducted an additional experiment. During the simulation of the routes, we added white noise to the realization of the travel and service times. The aim was to inspect if the results held up under these scenarios, thus representing the uncertainties in the probability distributions of both the travel and service times.

Noise was added using a normal random variables X with parameters µ X = 0 and standard deviations as a percentage of the expected value in travel and service times. This percentage varied from 0 to 0.3 with steps of 0.05. Figure 4.2 shows the average absolute gap error of percentile 95% for dierent levels of noise in travel and service times. Percentile 95% is selected since it is commonly used as a standard service level. The results show that the maximum absolute error is just above 5% when noise of travel and service times is 30% of their means. However, if noise levels are at 15% or less, the error is below 2%. Indeed, even at noise values of 20% the error stays under the 3% mark. Therefore, even with the presence of noise in times, the log-normal approximation reaches small gaps at percentiles estimation.

Multi-population Memetic Algorithm

To solve the SVRPTW we propose to use a Multi-Population Memetic Algorithm (MPMA) framework.

The choice of a MA is due to its exibility and performance in a variety of VRPs (for example see [40,47]). For the sake of clarity this section is structured as follows: rst the MA general structure 4.4. MULTI-POPULATION MEMETIC ALGORITHM is presented, this is followed by details of the main components of the proposed MA. Finally we will present the Multi-Population MA general framework.

MA general structure

A Memetic Algorithm (MA) can be described as the hybridization of a Genetic Algorithm with local search procedures [36]. Our MA structure follows the ideas of Prins [40] and is summarized in algorithm

4:

C ← Mutate(C) with probability pm.

5:

D ← Decode(C).

6:

D ← LocalSearch(D) with probability p ls .

Decoding

To decode a chromosome, we use the Split method of Prins [40]. Split works by considering an auxiliary

directed graph H = (W, Y ) with vertex set W = {W 0 = 0, W 1 , • • • , W i , • • • , W n }. W 0 represents a
dummy vertex and vertices W 1 . . . W n ∈ V c characterizes an ordered sequence of customers dened by a chromosome. An arc (W i , W i+d ) represents a feasible route for visiting the customers from W i+1 to W i+d , with its associated cost. The Split procedure nds the shortest path from vertex W 0 to vertex W n and subsequently the optimal set of routes associated to the chromosome. Normally this can be done using a shortest path algorithm. Nevertheless, constraint (4.3) precludes this implementation for the SVRPTW. As far as constraint (4.3), it is dependent on the set of routes which comprises a solution, this constraint must be guaranteed when decoding a chromosome. Shortest Path Problem (CSPP). In this CSPP the scarce resource is θ and each edge consumes ζ r units of the resource, where r is the route associated with the edge. The CSPP is solved using a Labeling Algorithm as the one employed in [14]. Thus, Split allows to decode the optimal routes partition for a chromosome, giving a solution which is feasible, i.e. respect constraints (4.1) to (4.3).

An example of the split method using the log-normal approximation is reported in gure 4.3. For this example it is assumed that α, β, γ are set to 95%. Customer time windows are reported in part (a) of the gure along with the chromosome. Moreover, a coecient of variation of 0.2 is used for the travel and service times. Their mean values are reported in part (b). Distances are assumed to be equal to the mean of travel times. For simplicity, the xed cost of each vehicle is assumed to be zero.

Besides, the time windows for the depot are dened by e 0 = 0 and l 0 = 140. Part (c) of gure 4. 3 shows the auxiliary graph of the Split, and the bold arcs represent the optimal solution to the CSPP.

In this auxiliary graph the arc (2, 1) for example, represents a route starting at the depot, visiting customers 4 and 1, and then returning to the depot. This route has an expected cost (see equation 4.7) of 71.8 and a probability of 93% of having no failures. The reader shall notice that only the arcs representing feasible routes (respecting the service levels α and β) are considered in gure 4.3 part (c). Lastly, part (d) displays the individual, or the decoded routes with the cost and the probability of having no failures in the solution.

Population

The population (P op) is dened as an ensemble of chromosomes, thus, an ensemble of coded solutions. P op is composed of P opSize chromosomes which are ordered in a decreasing way with respect to the cost of their detailed corresponding solutions. P opSize is constant, so there are always the same number of chromosomes in P op. The diversity of the population is controlled by mean of Campos et al. [6] distance measure. The latter is computed at the detailed solution level. To enhance diversication, clones (distance zero) are discarded, and chromosomes with a positive distance are allowed to enter. It shall be noted that the although this makes the MA a version with population control, the distance measure is only used to avoid clones. Further improvements can be performed by dynamically adjusting the minimum distance to accept the chromosomes willing to enter the population 1 . Furthermore, as populations are of a xed size, a random chromosome among the worst half of the population is deleted before the new one is inserted.

The initial population is created as follows: Four individuals are computed by using heuristics. The rst is created with the Clarke and Wright [9] heuristic, the second and third from Solomon [41] insertion heuristic (using two sets of parameters), and the fourth from Algorithm D proposed by Nagata and Bräysy [37] using as a starting solution the best individual among the Clarke and Wright and Solomon heuristics. All heuristics are run for the deterministic problem using the mean travel and service times as the true values. The four heuristic solutions are converted to chromosomes by concatenating their routes and then added to P op. Meanwhile, to achieve the value P opSize, the population is lled with chromosomes created with random permutations of the customers. If clones are built during this procedure they are discarded.

Crossover

Crossover procedure is used to create new chromosomes from those already present in P op. The crossover used in this work is the well-known Ordered Crossover (OX). OX works as follows: two chromosomes π 1 , π 2 are selected from the population and two random positions i, j | i < j ≤ n are chosen. The information between positions i, j is copied from π 1 to the new ospring. To complete the latter, π 2 is circularly traversed from position j + 1 to position j. The customers which are not already with the new child are then copied circularly from position j + 1 to i. The roles of π 1 , π 2 are exchanged to produce another ospring. Before performing the crossover, the MA selects π 1 , π 2 with binary tournaments, and then one of the two children generated by OX is randomly picked.

Local search and mutation

The local search (LS) is designed to improve the objective function of a solution by performing dierent modications on the solution itself. The proposed local search is based on the Variable Neighborhood Search (VNS) [35]. In this context Neighborhoods are structured in such a way that movements can be performed on a solution s to achieve a new solution s . Further explanation on neighborhoods denitions for VRPs can be found in Labadie et al. [25] and Bräysy et al. [5]. The neighborhoods used for the local search are: Or-opt, 2-opt movements, both intra and inter routes, and the CROSS exchanges in their inter-routes version.

LS starts by searching for a movement in the rst neighborhood which improves the solution. If such movement does not exist, it passes to the next neighborhood, however, if a new best solution is found in the current neighborhood, the LS starts again from the rst neighborhood. The process is repeated until no neighborhood can improve the solution. The order in which neighborhoods are explored is randomized and it changes every time the LS is performed. A movement is executed immediately if it improves the solution, therefore a rst accept criterion is used.

Since LS can be expensive in terms of running time, the neighborhoods exploration is constrained, Besides the exibility and performance of MAs they also allow a multi-population approach, e.g. using distributed evolutionary algorithms. In the latter, several populations can be used to enhance the diversication [1]. Our MPMA aims to enhance the diversication but also to speed-up and improve the solutions through the introduction of dierences between the problems tackled in each MA.

Algorithm 4 presents the framework of the proposed MPMA, it uses a set of MAs with the same presented structure, but with dierences in the problem handled by each one. In the rst line of Moreover, since Split using simulation is a highly time consuming task, we take advantage of every route simulated during the Split of the best MA chromosomes. This is done by saving each route into a pool, which acts as a list of feasible routes that have been validated by Monte Carlo Simulation. After all selected chromosomes have been evaluated, the routes in the pool are used to solve a set partitioning problem (see for example [33]) which uses equation (4.12) to guarantee the feasibility of the solution. By solving the set partitioning using the list of Monte Carlo Simulation validated routes, we guarantee that the routes in the solution respect the problem constraints. The solution retrieved by the set partitioning problem is returned as the best solution found by the MPMA.

MA dierences

To allow the communication of the MAs some changes are performed to the base MA. Every τ seconds, each MA sends two copies of its own chromosomes. One is randomly selected from the best half of P op, and the other from the worst. The destination MA of each child is randomly picked. Moreover, 4.5. NUMERICAL RESULTS destination MAs use the received chromosome in the next iteration of algorithm 3 instead of creating a new one by a crossover procedure.

Our MPMA uses three MAs (K = 3), namely MA1, MA2, MA3 with the same general structure presented in section 4.4.1 but with its own particularities. MA1 tackles the SVRPTW by embedding the log-normal modeling presented in section 4.3 while MA2 and MA3 work on a deterministic VRPTW. The idea behind MA2 and MA3 is that even if deterministic solutions might not perform well on stochastic environments, they may serve as good seeds to potential good stochastic solutions ( [3]).

MA2 and MA3 dier since MA2 travel and service times are set to their respective mean values, whereas MA3 service and travel times are set to the 75 th percentile that is s i = ψ i -1 (75%) ∀i ∈ V c and

t ij = φ ij -1 (75%) ∀ [i, j] ∈ E where ψ i -1 ∀i ∈ V c and φ ij -1 ∀ [i, j] ∈ E represent the inverse probability
functions of the service and travel times respectively. The 75 th percentile is selected after performing some preliminary tests. Increasing the service and travel times provides solutions with low probability of failures. A similar idea has already been used in the context of the VRP with stochastic demands ( [32]). Therefore, MA3 is used to provide robust solutions which can be fast improved upon to perform better in a stochastic environment. Moreover, MA2 and MA3 use labels to enhance the LS procedure (see [47]) and allow unfeasible solutions with the returns in time (time wraps as explained in [38]).

Additionally, MA2 and MA3 use the Bellman algorithm during the Split procedure (section 4.4.3) as done by Prins [40].

MA1 presents its own important considerations. Log-normal modeling allows MA1 to check feasibility, estimate the costs, and the underlying probabilities when using Split (section 4.4.3) and while performing the LS. This increases the complexity of the LS since feasibility checks require the evaluation of the whole route. Therefore, MA1 uses an additional strategy. Indeed, MA1 alternates its p ls each seconds by changing the probability of LS from zero to a given p ls and from p ls to zero value.

When MA1 p ls is greater than zero, it uses deterministic labels to eciently evaluate the pertinence of a movement (see [47]). The idea is that if a movement is unfeasible in a deterministic environment (times equal to their mean values), it is very likely that it will violate constraint (4.1) or (4.2). Subsequently it will be a waste of time to reevaluate the route using the log-normal approximation. If the movement is feasible for the deterministic model, then the log-normal approach is used to evaluate its feasibility and pertinence.

Numerical Results

Instances

To test the MPMA, results are gathered for modied Solomon [41] instances. Customer service times are set as Gamma distributed variables with mean equal to the deterministic service time given in the original instances, and a standard deviation inferred for a coecient of variation of 0.2. Edges travel times are set as log-normally distributed variables with mean equal to the length of the edge and with standard deviation derived from a coecient of variation of 0.2. This value is used when considering authors like Turner et al. [46] (as cited in [11]) previously stated that it varies from 0.15 to 0.25 when dealing with freeways and from 0.20 to 0.25 for principal and secondary roads. As the original instances do not state the kind of roads represented by the edges, the value 0.2 works overall for any type of road. Solomon [41] If two condence interval solutions overlap, it is assumed that they are, on average, equal. Detailed results of the best an average solutions are presented in tables 4.11 to 4.12.

Type 2 instances with random located customers are harder to solve than type 1 instances when considering the percentage of times that the best solution is found. This is evident for family RC2 with 100 customers, where the best solution is only found 11.25% of the time among the 10 runs.

Still, 4 out of 6 families found the best solution almost more than 30% of the times in 100 customer instances. In contrast, C2 family (with 100 customers) achieves 70% of the time for the best solution.

Overall the MPMA shows a good performance, retrieving on average 49.1% and 35.5% of the time the best solution for instances with 50 and 100 clients respectively. In terms of computational time, the MPMA needs at most three minutes on average to solve any instance (100 customers). The algorithm achieves these low times thanks to the log-normal approach presented in section 4.3 and because simulation is used limitedly. Indeed, the Split with simulation plus the set partitioning model take on average around a third of the total time for 100 customer instances.

Moreover, during the whole set of experiments it was evidenced the importance of considering the stochasticity of the parameters for the proposed problem. When the best chromosomes retrieved by MA2 were decoded in the uncertain environment (using Split with Monte Carlo simulation) they used more vehicles than the best solutions found by MPMA. Indeed, the CNV of the best MPMA solutions is just over half of the CNV of the best MA2 chromosomes. Thus, neglecting the stochasticity and using deterministic solutions can conduce to a huge increment of the number of vehicles. This is specially important for type 1 instances in which due to tight time windows, deterministic solutions perform very poorly.

The eects of considering multiple populations

To assess the eects of using multiple populations in the MPMA we compare the performance on 100 customer instances of the single MA1 and the MA1 combined with MA2 and MA3, to the full MPMA.

Table 4.8 presents the results for the dierent congurations. It is shown that using the MPMA with the three populations improves the solutions in terms of the average number of vehicles by around 6% when compared to the MA1 used alone. Moreover, the best solutions of the full MPMA reduces by more than 4% the same metric of MA1, MA1+MA2 and MA1+MA3. In general, congurations with more than one population achieves to nd solutions with less vehicles than the single MA1. Besides this, the best costs in terms of distance are achieved by the MPMA with three populations, however it presents the higher recourse costs. Despite the increase of nearly 15% on the average total time of the full MPMA in relation to MA1, the former conguration has the best overall results.

MPMA + Log-normal approximation comparisons

To further test the proposed MPMA + log-normal approximation, we compared it to the works of Miranda and Conceição [34] and Nguyen et al. [39]. These are two of the closest problems related to the problem at hand. [34] consider the SVRPTW with stochastic travel and service times and time windows. Also, late services are allowed but the service must start at or after the opening time window (e i ∀i ∈ V c ), furthermore the authors set the customers service level at 80%. Their Iterated 4.5. NUMERICAL RESULTS [34]. Service levels α, β, γ were set to 80%, 0%, and 0% respectively. No recourse action is considered for the problem so the cost of the route (equation 4.7) only takes into account both distance and xed costs of the vehicle. Travel and service times were modeled as normal variables same with the same parameters used in [34]. The MAs stopping condition is set to 25 seconds running time and late services are allowed as in [34]. Table 4.9 presents the results of the two methods, which are gathered over 10 runs. For each instance the method, the average number of vehicles (Avg. # Veh.), the average distance (Avg. DC), the average minimum service level (Avg.

Min SL), i.e. min P ( AT rj ) ≤ l i | r j = i, i ∈ V c , the average service level for all customers (Avg. SL), and the average and maximum service level error (Avg. SLE -Max. SLE) are reported. Service level errors are dened as the absolute dierence between service level estimated by simulation minus the service level calculated with the approximation.

Results show that MPMA uses near to one percent fewer vehicles than ILS. This decrease does not negatively impact the operational costs (distance) which are improved upon by almost 23%.

Concerning the service levels, the two methods appear to give similar results to the extent that the average service level of MPMA is inferior to ILS by only 0.08%. Moreover, it should be noted that MPMA ends up using Monte Carlo simulation enabling the method to validate the chance constraints.

In fact, MPMA guarantees a service level of α in all tested instances rather than the ILS method which guarantees this constraint in only ve out of seven instances. In terms of service level errors, the metric tends to be higher in MPMA although solutions with a higher number of vehicles and distance cost have a signicant reduction in service level errors. It cannot be concluded which approximation performs better using the service level errors as far as this metric is gathered from the best solution found by each algorithm at each run, and they are also structurally dierent (see operational costs).

Computational times were not scaled up since operating systems, language programming and hardware characteristics were dierent, however, it can safely be said that ILS is faster than MPMA using nearly 40% less time.

Nguyen et al. [39] work on the SVRPTW with stochastic travel times and time windows, using a Satiscing Measure Approach (SMA) to guide its solution process. The authors use a tabu search (TS) with their SMA and tested it on Solomon [41] instances. Travel times are assumed to follow an ambiguous distribution, however, for simulation purposes, Gamma distribution is used with a mean equivalent to the deterministic distance between two nodes and standard deviation derived from a coecient of variation of 0.5. Hierarchical objectives are used by the authors, in descending order of importance: maximize the total number of customers served, maximize the overall satiscing measure, minimize the total number of vehicles used, and minimize total distance traveled. The following changes were made to our MPMA by aiming to compare it with the results in Nguyen et al. [39]. Service levels α, β, γ were set to 50%, 0%, and 0% respectively since these constraints are not considered in [39], therefore we let MPMA be driven as an SPR. Travel times are assumed to be Gamma distributed with the same parameters used in [39]. The α value of 50% serves to discard routes with high probability of failures and thus to accelerate the MPMA local search. Furthermore, as a result of the associated travel times variability is high, the MA3 is set to the 60 th percentile that

is t ij = φ ij -1 (60%) ∀ [i, j] ∈ E where φ ij -1 ∀ [i, j
] ∈ E represents the inverse probability functions of the travel times. The MAs stopping condition is set to 40 seconds running time and late services were allowed as in [39]. Table 4.10 presents the summary of results gathered for MPMA (out of 10 runs), and TS (out of 1 run). TS runs only once since it does not have any randomized components.

For each family of instances the method, the average number of vehicles (Avg. # Veh.), the average distance or expected travel time if travel times are set to their mean (Avg. DC), the average expected tardiness (Avg. ET), and the average running time (Time) are reported.

Conclusions

This chapter presented a Multi-population Memetic Algorithm (MPMA) to solve vehicle routing problems with time windows (SVRPTW) and stochastic travel and service times. The approach incorporates the fact that failures can take place at the same time that a service level must be satised.

The estimation of the arrival times is based on the assumption that they can be approximated by a log-normal distribution.

The MPMA exploits the characteristics of dierent populations guaranteeing high service levels for customers service. The MPMA was tested by comparing it to two recent variants of the SVRPTW where late services are permitted. The MPMA found on average better solutions for the two comparisons albeit MPMA has higher computational times.

The results conrms that our formulation and resolution method are exible enough to deal with some variants of the SVRPTW and that the MPMA + log-normal approximation is a valid and eective method to solve SVRPs. Undergoing work is concentrating on the use of the proposed approximation within exact methods and a multi-objective version of the algorithm where service levels are considered as objectives rather than constraints.

Contributions

Preliminary results of this chapter were presented at MIM2016 conference: Chapter 5

Wind farms maintenance

Let us be careful! If in our climates, the industry can avoid to use the direct solar heat, necessarily, a day will come when, for the lack of fuel, it will be forced to return to the work of other natural agents. That the deposits of coal and petroleum will still supply for a long time their caloric power, we do not doubt it. But these deposits will undoubtedly exhaust themselves...(about coal energy) One can not refrain from concluding that it is prudent and wise not to fall asleep with respect to this in a misleading security Augustine Mouchot

Introduction

Wind energy is one of the most important sources of renewable energies. According to the Global Wind Energy Council (GWEC) [17], the year 2016 nished with an installed capacity of nearly 487 GW, representing a 12.6% growth when compared to 2015. In 2016 wind power also achieved to account for almost a quarter of the worldwide renewable energies capacity (including hydropower) [65]. Moreover, the GWEC projections show that wind energy installed capacity will increase by almost 70% in the following ve years reaching nearly 830 gigawatts (GW). Despite the rapid growth of oshore installations [17] onshore still represents nearly the 97.2% of the total wind worldwide capacity [16] in 2015.

In this context, Operation and Maintenance (O&M) activities represent an important cost [77,71,23,64] in wind farm projects representing as much as 25% to 30% of the total energy production costs (Ding et al. [23]). In the oshore context Shaee [72] and Raknes et al. [64], state that O&M costs range from 25-33% of the total life cycle cost. According to Scheu et al. [71], El-Thalji and Liyanage [27], and Racknes et al. [64], O&M costs include transportation costs, technician salaries and cost of repair actions and spare parts, as well as loss of revenue caused by production stops. Therefore, optimizing the O&M activities is a mandatory task to improve the competitiveness of prices and attractiveness of the wind energy projects.

This chapter is divided in two parts. The rst one starts by making a review about maintenance activities in wind farms with a focus on operational decision level of analysis (section 5.2). Then, section 5.4 introduces a multi-objective maintenance scheduling problem in the context of onshore wind farms. The multicriteria approach is motivated by the presence of dierent stakeholders with conicting objectives in this kind of applications. The second part explores the strategic level of maintenance activities in wind farms in section 5.5. Then, in section 5.6 the eects of the operational scheduling problem within a long horizon and its relation with the maintenance strategy are studied.

The chapter ends in section 5.7 with concluding remarks on the found results and gives some research clues.

REVIEW ON OPERATIONAL MAINTENANCE ACTIVITIES

Review on operational maintenance activities

Wind turbines are the systems that allow to transform the kinetic energy in the air into electric energy [50,10]. Although wind turbines can be used as a stand alone system to provide the energy (or part of it) for particular houses, farms etc, most of the capacity is installed in wind farms to produce energy in an industrial scale [50]. Wind farms are the collection of several turbines working on grids to produce energy.

Maintenance activities are performed to keep systems working, and to x and prevent possible failures [73]. The fact that wind farm projects are devised for long periods of time [60] (20 to 40 years) makes the maintenance and repair activities a central issue to keep the turbines working. Maintenance actions can include inspection, changes of consumables (greasing, lubrication, oil lters), oil sampling, re-tightening of the bolts, repair, overhaul, and replacement of parts [5]. According to Manwell [50] the scheduled and unscheduled maintenance reduce the availability1 of wind turbines. Furthermore, maintenance activities have an inuence in the failure rates of wind turbines, since poorly maintained turbines can have higher failure rates [15]. Moreover, both maintenance activities and failures imply losses in production due to downtime [71,27] and generate important costs (spare parts, technicians, etc.).

The interaction of multiple actors involved in wind farms has eects on how maintenance activities are carried. Markard and Petersen [51] recognize ve stakeholders on wind farms projects, namely: turbine manufacturer, project developer, investor, operator, and load management and power distributor. The exact conguration of the stakeholders roles depends on each project. Still, dierent actors can play multiple roles in the project. For instance, an investor can also manage the operation of the wind farm. Situations in which roles are held by dierent entities might derive in cases where, one actor pursuing its objectives harms other actors targets.

The challenge raised by improving the maintenance of wind farms has attracted a large number of academic communities with dierent perspectives on how to achieve this goal. Some studies have focused on the materials used for the components (particularly the blades to avoid failures) of the wind turbines [9,43], wind farm topology [44], fault diagnostics techniques [26,29,37], prognostics [81,45,13,85,82,12,91], scheduling of activities and resources [8,40,79,19,35,31,30], or maintenance strategy selection [81,14,11,62,68,25]. These studies aim to improve the maintenance activities according to criteria in dierent moments of the wind farms life-cycle. Therefore, a powerful way to classify the dierent scientic contributions in this eld is to consider which decision-making level is involved.

As in production management, maintenance problems can be categorized in dierent decision levels according to the time horizon involved. Shaee [72] uses a three-echelon classication of decisionmaking with strategic, tactical and operational echelons as its levels. Strategic decisions aect in the long-term the maintenance activities. In this echelon, one can nd aspects such as the location [67,89],

layout [2,32,88], or the type of maintenance strategy [42,23,11,76]. Tactical level is related to midterm decisions which can vary from year to year, e.g. spare parts inventory management [36,82].

Operational echelon deals with short-term maintenance decisions. Within this, one can include the scheduling of resources to perform the required maintenance tasks [8,63,79,19,35,31,64].

The attention given by the Operations Research community to the Operational level (short-term maintenance decisions) focused mainly on the maintenance scheduling problem. As pointed out by Kovàcs et al. [40], maintenance scheduling plays a major role in the turbines availability and costs.

The rst study presented in this chapter is dedicated to such types of operational problems and starts with a review of studies dealing with this decision level.

OPERATIONAL MAINTENANCE LEVEL

Table 5.1: Summary of wind farms maintenance scheduling works, objectives, and types of models Based on the considerations presented in sections 5.2, in this section a model and a resolution approach are presented to deal with the multi-objective maintenance scheduling problem for onshore wind farms. It is assumed that O&M operator and investors are the two parts involved in the maintenance scheduling. O&M operator is paid to perform the maintenance activities, and it is assumed that he seeks to minimize its costs. Meanwhile, the investors expect to maximize the energy production. The problem is forthwith formally dened.

Problem Denition

A set of J turbines indexed by j are considered. Turbines might be distributed among dierent wind farms but it is assumed that they are close enough to be reached in despicable times. Furthermore, turbines might require more than one type of maintenance. Each task i ∈ I represents a maintenance activity. It is assumed that the set of tasks is specied prior to solving the problem. Also, each task is associated with a turbine j ∈ J, and I j dene the subset of tasks to be performed in turbine j.

To execute each task i a number of χ is technicians with the skill s are necessary. Actually, each task requires one or dierent skills held by technicians to perform the maintenance, e.g. mechanical, electrical, electro-mechanical, etc. S denes the whole set of skills. Besides, an execution time β i characterizes every task i. Time β i corresponds to the time elapsed from stopping the turbine until it is restarted and veried again. It is assumes that all maintenance activities in I requires that the associated turbine is stopped during the task execution. Additionally, every started task is performed until nished, and all tasks must be performed during the planning horizon.

Furthermore, every task is associated to a time window [e i , l i ], where the maintenance should take place. The opening time window is a hard constraint, so the task must start at e i or latter. The enforcement of this condition accounts for the necessity of spare parts or any special equipment (e.g. cranes) to perform the task. Other spare parts or consumables are assumed to be available at e i . The time window closure is a soft constraint; therefore, the task should preferably be nished at most at l i . In the case it nishes later, a penalization cost α is considered and is proportional to the delay.

Usually, periodic tasks can be performed all along the planning horizon, while corrective maintenance can be constrained for special equipment or special considerations based on their impacts.

Time is discretized in equally length periods t ( [19,31,35]). The set of all time periods is dened by T . Every time period stands for the same amount of time, e.g. one minute, two hours, etc. During each time period t, a turbine j generates an amount of energy or utility per production denoted as Θ jt . Moreover, the additional representation of time as done by Froger et al. [31] is also used. In this one, a set D denes the set of days in the planning horizon. Every day d ∈ D has a number of workable t periods called τ d . Furthermore, the time intervals within a day are divided as normal working and extra working periods. The subset of time periods representing extra working periods is dened as T e . At the end of each working day in the horizon plan, turbines continue to produce an amount of energy or utility per production Υ jd ∀ j ∈ J until the next day. Besides, for safety reasons, maintenance task can be only carried out during certain periods of time when weather allows it, e.g.

wind speeds are below a given value. To address this fact, the parameter ρ it takes value one if the safety conditions are met and zero otherwise. If a task i is being performed during time t -1 and ρ it = 0, the technicians must stop until the safety conditions are met again.

A limited set of technicians is dened by P . Technicians travel to the turbines to perform their assigned maintenance tasks. All technicians perceive a xed salary w p , which is linearly dependent on the amount of skills they possess. Furthermore, a technician receives an extra wage ew p for each extra working hour. It is assumed that a worker with a higher salary will enjoy a higher extra payment, i.e.

w p > w p then ew p > ew p ∀p, p ∈ P . Besides, the parameter π ps takes value one if the technician p has the skill s and zero otherwise. Let us also dene l s as the minimum wage among the workers p ∈ P | π ps = 1 ∀s ∈ S. That is, the less costly wage for a technician with the skill s. When a technician is assigned to perform a task, he must be present at the turbine until the task is complete. Moreover, technicians can only work at one task in each time period.

It is considered that O&M operator costs comprise the spare parts, the special equipment to perform tasks and the technician's wages. Moreover, three main components are examined. The rst one is the dierence of wages perceived by technicians assigned to tasks that could be performed by less costly ones. The second one, is the extra hour wages paid to technicians, and the third one, is a penalization due to maintenance tasks outside the time window. Spare parts and special equipment costs are disregarded since they are assumed to be unavoidable by the O&M operator. Hence, the objective is to minimize the costs of the O&M operator while maximizing the amount/prot of energy production.

Mathematical model

The problem described in section 5.4.1 is now formalized as an Integer Linear Problem (ILP). The following variables are thus dened:

x it :    1 if task i is scheduled to begin at time t 0 otherwise z it :    1 if task i is scheduled at time t 0 otherwise e it :    1 if task i is nished at time t 0 otherwise y pi :    1 if resource p is assigned to task i 0 otherwise v pti :    1 if resource p is scheduled at time t to task i 0 otherwise u id :    1 if task i is scheduled in day d 0 otherwise γ jt :    1 if turbine j is able to produce energy at period t 0 otherwise η jd :    1 if turbine j is able to produce at the end of day d 0 otherwise ζ i :
The number of days task i is delayed with respect to b i Using the abovementioned variables and the characteristics described in section 5.4.1, the ILP is dened as follows:

max Z 1 = j∈J t∈T Θ jt γ jt + j∈J d∈D Υ jd η jd (5.1) min Z 2 = i∈I p∈P s∈S κ i y pi (w p -l s ) + p∈P t∈Te i∈I ew p v pti + i∈I αζ i t∈T x it = 1 ∀ i ∈ I (5.2) 5.4. A MULTI-OBJECTIVE APPROACH TO THE MAINTENANCE SCHEDULING PROBLEM t∈T e it = 1 ∀ i ∈ I (5.3) z it = x it + x it-1 -e it-1 ∀ i ∈ I, t ∈ T (5.4) x i0 + e i0 + z i0 = 0 ∀i ∈ I (5.5) t∈T π ps y pi ≥ χ is ∀i ∈ I, s ∈ S (5.6) t∈T z it ρ it =β i ∀i ∈ I (5.7) dτ d t=1+(d-1)τ d z it ≤ τ d u id ∀ i ∈ I, d ∈ D (5.8) u id = 0 ∀ i ∈ I, d ∈ D | d < e i (5.9) ζ i ≥ du id -l i ∀ i ∈ I, d ∈ D (5.10) z it + y pi ≤ 1 + v pti ∀ p ∈ P. t ∈ T, i ∈ I (5.11) v pti ≤ z it ∀ p ∈ P. t ∈ T, i ∈ I (5.12) t∈T v pit ≤ M y pi ∀ p ∈ P, i ∈ I (5.13) γ jt ≤ (1 -z it ) ∀j ∈ J, t ∈ T, i ∈ I j (5.14) η jd ≤ 2 -(u id + u id+1 ) ∀ j ∈ J, d ∈ D, i ∈ I j (5.15)
x it , z it , e it ∈ {0, 1} ∀i ∈ I, t ∈ T (5.16)

y pi ∈ {0, 1} ∀ p ∈ P, i ∈ I (5.17) v pti ∈ {0.1} ∀ p ∈ P, t ∈ T, i ∈ I (5.18) u id ∈ {0, 1} ∀ i ∈ I. d ∈ D (5.19) γ jt ∈ {0, 1} ∀ j ∈ J. t ∈ T (5.20) η jd ∈ {0, 1} ∀ j ∈ J. d ∈ D (5.21) ζ i ∈ Z + ∪ {0} ∀ i ∈ I (5.22)
The objectives in equation (5.1) state that the amount/utility of energy produced is to be maximized while cost must be minimized 2 . This last is composed by three parts, the total dierence between the salaries of technicians assigned to a task and the minimum wage of a technician oering the same skill. The second part accounts for the extra time periods wages, and the third one is the penalization for not nishing the task within the time window. Constraints (5.2) and (5.3) ensure that all tasks are started and nished within the time horizon. The consistance between the start, the execution, and the end of a task is guaranteed by constraints (5.4), which ensure that tasks are executed until nished when they have been started. Constraint (5.5) state the initial conditions by considering that no task starts, nishes, or is assigned during time period zero. The number of technicians with the required skills to perform each task is imposed in constraints (5.6). These equations have a sense of greater or equal for some special cases, e.g. consider a task i which requires two skills s 1 , s 2 with a number of technicians of two and one respectively. If the equation is set to equality two technicians p 1 , p 2 both with skills s 1 , s 2 cannot be simultaneously assigned to the task since it will violate the equation for skill s 2 . Still, this scenario where p 1 and p 2 are assigned to perform task i is valid and therefore the constraint sense is kept. Constraints (5.7) guarantee that the amount of time periods spent on CHAPTER 5. WIND FARMS MAINTENANCE a task will be sucient to nish it, the parameter ρ it ensures that the weather allows to execute the task during these periods. The coupling between tasks executed during a day and the scheduled time periods is considered in constraint (5.8). Time windows are dened in constraints (5.9) and the (5.10), the rst ones safeguard that tasks must start after their associated opening time window, and the second ones, keep track of the number of days a task is delayed with respect to its time window closure. Equations (5.11) to (5.13) couple the technicians assignment to a task with its execution times. The parameter M in constraint (5.13) is used as a big M value. Moreover, constraints (5.14) and (5.15) determine if a turbine can produce energy at normal time periods, and at the end of the days respectively. Finally, equations (5.16) to (5.22) x the decision variables.

The Epsilon constraint approach

Multi-Objective Optimization Problems (MOOP) are dened as problems where at least two, often conicting objectives are to be optimized at the same time. Multi-objective optimization methods dier from single objective ones in which only one solution is found to minimize/maximize the criterion at hand. In fact, in MOOP more than one solution can be found, especially when the method used deals with Pareto-optimization. Thus, solving MOOP is a process of nding the ensemble of solutions called Pareto ecient or non-dominated. These solutions are those for which no objective can be improved without worsening at least another objective. The whole set of non-dominated solutions is called Pareto Front [48].

Several approaches are discussed in the literature to solve MOOP, e.g. weighted global criterion [52], goal programming [80], epsilon constraints [55], etc. To address the model the last method is selected. This approach is employed to construct an ensemble of Pareto Ecient Solutions (PES).

In general terms, the epsilon constraint works by iteratively solving single objective problems. Consider the following multi-objective problem max Z (x) = (z 1 (x) , z 2 (x) , . . . , z n (x)) , x ∈ Ω. Each z i (x) ∀ i = 1 . . . n represents an objective, and the condition x ∈ Ω stands for the feasible set of solutions. Epsilon constraint method solves a group of problems max z j (x) | z i (x) ≥ ∆ i ∀ i = 1..n ∧ i = j, x ∈ Ω by changing the values of ∆ i . Each solution found is ecient and kept devising (partially) the Pareto Front. Although simpler approaches can be used such as the Weighting Method in which the dierent objectives are reduced to a simple objective using weights for each component, the epsilon constraint is preferred for the following reasons. First, epsilon constraint can nd non-supported solutions, i.e. solutions not in the convex envelope of the Pareto Front. Second, there is no need to scale the objectives. Third, epsilon constraint iterations can be coded so new ecient points are found at each iteration, this can avoid iterations which need considerable amounts of time for solving the ILPs.

Algorithm 5 Epsilon Constraints adapted method for maximizing

Z = (Z 1 , -Z 2 ) Require: numSteps 1: Solutions ← { } 2: {z1 * , z2 . } ← max Z1 3: {z1 * , z2 * } ← max -Z2| (Z1 = z1 * ) 4: Solutions ← Solutions ∪ {z1 * , z2 * } 5: {z1 . , z2 * * } ← max -Z2 6: stepSize ← (z2 * -z 2 * * ) numSteps 7: i ← 1 8: while (z2 * -i • stepSize) > z2 * * do 9: {z1 * , z2 . } ← max Z1| (Z2 ≤ z2 * -i • stepSize) 10: {z1 * , z2 * * * } ← max -Z2| (Z1 = z1 * )

11:

Solutions ← Solutions ∪ {z1 * , z2 * * * } 12:

i ← i + 1
13: end while 14: Sort and check(Solutions) Algorithm 5 shows the general steps of the used implementation of the epsilon constraint method.

A MULTI-OBJECTIVE APPROACH TO THE MAINTENANCE SCHEDULING PROBLEM

The procedures max Z 1 and max -Z 2 solve the ILP described in section 5.4.2 maximizing the produced energy and minimizing the costs respectively. Moreover, the procedures retrieve both the values of Z 1 and Z 2 given the solution of the ILP. Line 1 starts by initializing an empty array in which the solutions will be kept. In line 2, the quantity/utility of the energy production is maximized, while line 3 minimizes the costs with the additional constraint that energy production must match the one found in line 2. The solution is then saved in the proper array of ecient solutions. It shall be noticed that it is necessary to solve both problems (energy and costs), to retrieve an ecient point. Preliminary tests show that solving only the problem maxZ 1 gives suboptimal solutions in terms of costs. The algorithm continues solving the problem max -Z 2 in line 5, this is done to dene the interval within the objective Z 2 is comprised, i.e. [z 2 * , z 2 * * ]. In line 6, a step size is calculated to use it for the epsilon constraint method. This value depends on the minimum and maximum value attained by Z 2 as well as a parameter called numSteps. This last allows to control the number of iterations performed during the loop. Higher values for the number of steps permit to better determine the Pareto Front at the expense of higher computational times. Nevertheless, the value for this parameter does not guarantee that a dierent ecient point will be found for each possible step value, it works as an upper bound.

Between lines 8 and 13 the main loop iteratively solves the ILPs for energy and costs using the proper constraints, while saving the ecient solutions. Finally, in line 14 the ensemble of solutions is sorted and solutions are checked for dominance. If solutions are optimal, they are ecient, nevertheless, if optimality is not guarantee, solutions are only potentially ecient.

Results

To test the model, the instances proposed by Froger et al. [31] are used. This testbed is composed by 160 instances grouped in families (5 instances per family) described as a_b_c_d_e where a, b, c, d, and e refer to the number of time periods in the planning horizon, the number of periods per day, the number of skills considered, the number of tasks and technician-to-work ratio, respectively.

Accordingly each core characteristics can have dierent values such as: time horizon lengths (10, 20 or 40), time periods per day (2 or 4), number of tasks (20, 40 or 80), number of skills (1 or 3), and the technician-to-work-ratio (A and B). Among these, only the subset of 40 instances with at most 20 tasks and 20-time periods are used. Moreover, for instances with two-time periods per day additionally period to stand for extra hours is considered. This extra time period has an implicit duration of four hours. The same procedure is performed for instances with four-time periods, however for these last, two extra time periods are added per day (each one representing 2 hours). To assign an energy utility to these time periods, the following procedure is employed. The original Υ jd for each day and turbine is taken and divide by 16 (number of hours between workable days). Then, extra periods use this coecient multiplied by the number of hours they stand for, to determine the amount/utility of energy production. It is assumed that weather conditions are safe to perform maintenance tasks during extra time periods, i.e. ρ it = 1. Additionally, since original instances consider multiple modes or ways to perform each task, a single mode is randomly selected. This picked mode includes the information on the number of technician per skill required to perform the task, and its duration (number of periods).

Besides, it is assumed that transport times between every pair of turbines are negligible. Tasks time windows are not considered in Froger et al. (2017) instances, thus e i = 0, l i = |D| ∀i ∈ I.

To constraint the computational eort required to solve the problems the amount of time that epsilon constraint method expends on the ILPs is limited. For problems in lines 2, 3 and, 5 in Algorithm 1, the time limit is set to 2000 seconds. ILPs in lines 9 and 10 are limited to 500 seconds.

The rst problems are led to run for more time since they are used to create the interval in which epsilon constraints will be dened. Moreover, the parameter associated to the number of steps is set to 50. This value guarantees a good trade-o between the approximation of the Pareto Front and the running times. Both time limits and the number of steps were selected after several preliminary tests.

Since times are constrained, optimal solutions are not guaranteed, therefore, the GAP metric for each single objective problem is saved.

All runs are conducted on a Dell Latitude E6420 personal computer with Intel Core TM i7-2760QM @2.4 GHz, running Windows 7 Professional 64 bits. The algorithms were coded on Java and compiled with JavaSE-1.8_45, with maximum allocated memory of 1 Gb. To solve linear problems, the Java interface with Gurobi 7 (2017) optimizer is used. each day and only one skill is considered. Moreover, the instance contains 20 tasks with a regular technician-to-work-ratio, that is, technicians can perform all the tasks during the planning horizon.

16 non-dominated points (solutions) are found, although only two are proven optimal, thus ecient.

Moreover, 12 points present a small gap in the energy component. Still, this gap is on average of 0.03%. The other two points present an average gap of 0.04% in the cost component. The graphic shows an overall concave behavior, displaying a smaller energy utility/production when costs are small and bigger production at the expense of higher costs. However, it shall be noticed that not all points rely on the convex hull, e.g. point nine, counting from the left-below part of the graph, showing hence the conicting nature of the objectives. All the tested instances present a similar behavior in terms of shape. Furthermore, it shall be noticed that for gure 5.1, an increase of 115% in the minimum cost can have an impact of almost 4% in the energy utility/production. Table 5.2 summarizes the results found for testbed instances. For each family of instances is reported, the average time (Avg. Time) per solved instance, the average number of non-dominated solutions found (Avg. Solutions), the average number of solutions proven to be ecient (Avg. E.

Solutions), the average gap for non-optimal solutions in terms of costs (Avg. Gap C), the average gap for non-optimal solutions for the energy objective (Avg. Gap E), the average number of solutions in which costs are not optimal (Avg. CNOP), and the number of solutions in which energy is not optimal (Avg. ENOP).

Table 5.2 shows that on average, instances take around 2.37 hours to be solved. A signicant 5.5. MAIN CONTRIBUTIONS ON STRATEGIC DECISION LEVEL increase in computational time is seen when the number of periods is incremented. The reason for this behavior is the increase in the size of the ILP models. In terms of the number of solutions, a very limited number is found for the whole set of instances, averaging 11.35 points. A reasonable explication to this limited number of potentially non-dominated solutions is that solution space of the problems is highly constrained. Therefore, the number of solutions, and more important, ecient solutions are limited.

Among the solutions found by the epsilon-constraint method, nearly one third are proven to be ecient. As little as this number might be seen, this fact can be explained by the complexity in solving the single objective ILP. One can see that the cost objective is by much the one with the higher gaps, averaging 13.9%. Furthermore, the gaps on cost component show a considerable rise when more skills are considered in the problems. This is especially important in instances with higher number of time periods, for instance family 20_4_3_20_B reaches a maximum of 34.6%.

Energy objective contrary to costs objective, show an excellent performance. Through the 40 instances, the average gap is only 0.6%. Type B instances (regular technician-to-work-ratio) consistently outperform type A (tight technician-to-work-ratio) instances in this metric. Despite the better performance in terms of energy gaps, results show that the number of solutions not proven ecient is mostly due to energy objective function. The average number of solutions for which energy presents a positive gap doubles the same metric for costs objective. Therefore, it is safe to say that single objective Z 2 is solved to optimality more consistently than Z 1 . However, when dealing with sub-optimal solutions, the average costs (Z 2 ) gap is over 20 times bigger than the energy (Z 1 ) average gap.

The results reect the importance of operational decisions with the wind farm performance indicators. The fact that schedules can change in as much as 4% the amount of energy produced in the current context is a matter of most importance. Actually, in closer to real life decisions, the time required to perform the activities, or the activities themselves might be stochastic. In such situations the changes in the production can be more pronounced. The same eect may be produced when bigger instances are solved.

Toward maintenance strategy selection

Evaluating which maintenance strategy is better adjusted is one of the most important decisions concerning O&M for a project [42,23,72,76]. Having no strategy on how to perform maintenance tasks usually derives in a very bad performance. Actually, Van Bussel and Schöntag [84] showed that a maintenance strategy is required to enhance the availability of wind farms. Additionally, the maintenance strategy policy denes scheduling parameters, in consequence the selection of this strategy will aect the solution of the maintenance scheduling at operational level. To evaluate the strategies, this part analyzes the eect of dierent maintenance strategies on the operational cost and energy production. Before presenting the developed model, a quick revision on strategic maintenance problems is presented in the following section.

Main contributions on strategic decision level

Strategies for maintenance activities rely essentially on two main decisions: when to take maintenance actions and what kind of maintenance (preventive,corrective,opportunistic, etc.) to perform. In the particular context of wind farms, Shaee [72] uses a general classication based on two major classes:

failure-based (reactive response), and proactive maintenance.

Failure-based strategy uses the idea of performing maintenance tasks only when failures take place, deriving in corrective actions. This actions can be performed immediately or delayed to a future time [5,86]. Nevertheless, there is no standardized way for dening a failure in the wind energy. For example, Caroll et al. [15] dene a failure as a visit to a turbine, outside the scheduled operation, in CHAPTER 5. WIND FARMS MAINTENANCE which material is consumed. Meanwhile, Wilkinson et al. [87] assume failures as events that require manual intervention to restart the machine generating a downtime of at least one hour. Consequently, it is important to dene what failures are in advance so maintenance strategies can be eciently selected.

Alternatively, Shaee [72] denes proactive maintenance as any maintenance that is performed before a failure occurs, aiming to avoid possible unexpected failures which usually derive in more costly actions [41,53]. Within this classication, two categories can be found: preventive and predictive maintenance. Preventive maintenance is used to refer to time-based maintenance schemes, whereas predictive maintenance is equivalent to condition-based maintenance schemes (decision is made according to the observed health state of the system).

Predictive maintenance usually requires condition-monitoring data. The collected information used to prevent failures before they become more costly to x. Condition-maintenance-systems are analyzed in Nilsson and Bertling [59] showing that an increase on the availability of around 0.43% can compensate the costs of such systems. The eects of false alarm of condition monitoring systems are studied by May and McMillan [56] showing that even with a decrease in condition monitoring system reliability, the availability of a wind farm is hardly changed. A recent review on condition-based maintenance optimization for stochastically deteriorating systems in a general context can be found in [3].

Prognostics can also be used in maintenance decision making. Examples of maintenance strategies based on the estimation of the Remaining Useful Life (RUL) of the asset can be found in [12,46].

Zhao et al. [91] analyze Supervisory control and Data Acquisition Systems (SCADA) information from turbines to detect normal and anomaly data of turbines generators during their runtime. This allows to make good estimations of the RUL, and provides sucient time before the appearance of failures so that a schedule (or a plan) can be implemented before they take place. A recent review on prognostics techniques for wind turbines can be found at Leite et al. [21].

When dealing with multi-unit systems, another type of maintenance is the one called opportunistic.

This type of maintenance aims to take advantage of breakdowns or stops of the system, due to a failure or a proactive maintenance, to perform additional proactive maintenance tasks on other components [81,22,12,1,47]. In the oshore context, Tian et al. [81] comment the possibility to perform additional maintenance tasks on nearby turbines when one was selected for a maintenance task. The main ideas is that grouping maintenance activities can lead to a cost reduction by sharing maintenance setup costs and downtime duration.

Maintenance strategy selection

Qualitative and quantitative approaches can be used to select the more appropriate strategy according to Andrawus et al. [6]. The former relies on subjective opinions [7,70] and may not be adequate to select the best possible strategy. Quantitative approaches aim to pick the best strategy based on mathematical models and are by far more studied than the qualitative ones.

In quantitative models, strategies are usually compared together. To evaluate which strategy is the best, two criteria are often considered: reliability (and availability, see [49]) maximization and cost minimization [23]. Other metrics such as environmental ones [53] has been used to compare dierent policies.

Strategies such as failure-based has been studied by Ding [24] given that corrective maintenance are performed after more than one failure occur. The author shows that this can conduce to an improvement in terms of cost compared to executing maintenance actions immediately after the rst failure. Nevertheless, corrective actions have shown to have considerable negative eects [42,41] on wind farms costs and availability. Alsyouf and El-Thalji [4] concludes that proactive maintenance is more suitable for wind farm maintenance problems. Condition-based and opportunistic maintenance are proposed in Van Bussel et al. [83] adapted to oshore context. Thus, it is safely to conclude that 5.5.2 Complex models for maintenance strategy selection Dierent authors have employed advanced simulation techniques to evaluate maintenance strategies under the interaction of wind farms elements. Normally, wind farms are simulated as an ensemble of turbines dened as single or multi-components systems. Moreover, weather models are used to recreate the wind speeds and wave heights (oshore) faced by the turbines. Weather data has two purposes, rst it enables to deduct the amount of the energy produced (wind speeds), and second it allows to check if maintenance safety conditions are met to conduct the tasks. The amount of energy produced depends on the rating of the turbine, speeds, air density, and many components, and is usually included as a function or Power Model. Furthermore, the inherent characteristics of each context creates dierences in the considered elements of the models. For example, maintenance operations are not constrained by transportation resources when the authors consider onshore wind farms. Conversely, in oshore this is an important characteristic since dierent types of vessels, corresponding to scarce resources are needed to perform the maintenance.

Table 5.3 presents a summary of the works using simulations approaches to evaluate maintenance strategies for wind farms maintenance. These references are chosen dependently in if they consider limited resources, incorporate weather models, or other operational or tactical levels. The selection criteria is used as this complex models allow to further research the interaction between strategy selection and lower levels decision problems. For each work the following nomenclature is used to classify the main aspects of the works: Onshore (ONS), oshore (OFF), number of turbines (NT), turbines power generation model (PG), weather model (WM), number of components per turbine (NC), failures modeling (FM), corrective maintenance (CM), preventive maintenance (PM), predictive maintenance (PAM), limited technicians (LT), limited resources for transportation (LRT), the presence of energy metrics to evaluate the strategy (EO), the evaluation of strategies using cost metrics (CO), and the consideration of inventories (IC). The abbreviation NS stands for non specied.

How the turbines are modeled is an important aspect of the models dealing with how the failures appear (and are solved). Either a single or a multi-component system are used as seen in column NC, table 5.3. The former considers only one of the components of the turbine (or aggregate them all in a single one). Meanwhile, multi-component systems permit to consider the individual degradation of each selected subsystem, allowing to consider strategies of opportunistic maintenance. Singlecomponent systems have been addressed by Carlos et al. [14] who propose to model the turbine using a reliability model in a onshore context. A multi-objective problem is designed by the authors using maintenance costs and the amount of energy produced as the two selected criteria. The failures process is considered to be Weibull distributed. Moreover, the authors consider that maintenance activities are imperfect and use a Proportional Age Set-back as in [54] while Wind velocity is modeled using a Weibull distribution. Furthermore, Carlos et al. [14] consider that the maintenance time are also random variables following an uniform distribution. Although the authors do not justify this choice, it can be explained by the deterioration of the system, or the technicians eciency. Using a simulation of one turbine, the authors derive an optimum time interval to perform maintenance activities. Similarly, Sahnoun et al. [66] design turbines as a single component for which degradation is modeled with ten status scale in an oshore context. Turbines status changes according to a function that includes, an exponential variable representing the components time to failure, and the eect of weather on turbines degradation. This is very important since among the studied works, this is the only one that considers the degradation as a function of the weather. Three strategies are compared using a multi-agent systems simulation: periodic, condition-based maintenance, and a hybrid version combining both. Also, maintenance task are scheduled based on rules such as the turbine with the maximum degradation, or if a time window is available to perform the task. Simulating a wind farm with 80 turbines, the authors show the hybrid strategy is the most eective considering the costs and the produced energy. Although Byon et al. [11] also model turbines as one component, this last stands for a subsystem of the turbine. This contrast with the approach of [14,66] where the single component aimed to model the whole turbine. 5.3). These usually include the gearbox and the generator, and other components such as the blades, shafts, or electrical ones.

These components are usually selected since they represent most of the downtime experienced by a turbine [15], or have the larger impacts in terms of costs. Andrawus et al. [6] include turbines subsystems such as the main shaft, main bearing, gearbox and generator. The last two include also some sub-components such as the bearings. The lifetime of subsystems and components are modeled using Weibull distributions. The parameters are derived from maximum likelihood estimators from reported data failures. A simulation is carried out for 26 wind turbines composing a wind farm.

Maintenance crews are limited and inventories policies are also addressed. Maintenance strategies are selected for each subsystem based on the parameters of the Weibull distribution, but no comparison among dierent strategies is given. Meanwhile, the benet of installing condition monitoring systems for the onshore context is tackled by McMillan and Ault [58,57]. Similarly to [6], the authors consider turbines as a four-component model which includes: the generator, gearbox, blades, and electronic related parts. This study uses the items presenting the most important failures, but also focuses on monitored components. It is assumed in [58] that monitoring reveals the true component state.

Moreover, the turbine is modeled with a Markov Chain taking account the state of each component.

The authors conclude that for most onshore wind turbine components, condition-monitoring is costeective. Furthermore, according to the authors the same conclusions should hold in the oshore case.

An extension of the work in [11] is presented by Perez et al. [62] to consider multiple components in the onshore context. As in other related works [6,58,81], the number of components is limited to four, namely the gearbox, power generator, blades, and control system. Moreover, the authors include the inherent constraints for schedules due to limited maintenance teams and to lead times in spare parts. Condition-based-maintenance presents better performance when compared to periodic maintenance. Additionally, a condition-based-maintenance strategy which also includes opportunistic maintenance, presents the best results in terms of costs and number of failures. Santos et al. [68] consider a one turbine model. The turbine is simulated by a four-component that includes the gearbox, generator, pitch system, and rotor. The authors optimize the maintenance strategies through the use of generalized stochastic Petri nets coupled with Monte Carlo simulation. The authors optimize two 5.6. MAINTENANCE STRATEGIES: RELATION WITH OPERATIONAL PLANNING parameters: the preventive repair threshold and the age reduction ratio. The rst is a proportion of the mean time to failure of the component to perform a preventive task, while the second stands for how much of the age of the component is reduced with the maintenance task.

A multi-objective opportunistic maintenance optimization considering limited resources is tackled by Abdollahzadeh et al. [1]. Turbines are dened by multiple components (rotor, main bearing, gearbox, generator) with lifetime modeled by a Weibull distribution. Using the reliability of the components, the authors dene several thresholds to execute opportunistic maintenance. By mean of a Particle Swarm Optimization, the threshold values are optimized to maximize the energy produced while minimizing the maintenance costs. Several conclusions can be retained from this work. First, it is shown that opportunistic maintenance coupled with preventive actions outperforms the classical corrective strategies. Second, the number of maintenance teams has an important eect in the Pareto ecient solutions and the addition of maintenance teams can greatly increase the energy produced.

Third, maintenance strategies must be compared with basis on more than one objective, as far as these can be in conict and lead to dierent conclusions.

A variant of the work of Sahnoun et al. [66] is introduced by Dahane et al. [18] considering the impact of spare parts re-manufacturing. Turbines are modeled by two components, a single turbine representation (as in [66]) and the gearbox. Both elements use similar degradation functions that include the eects of weather or the accelerated degradation of re-manufactured components. The comparison is performed on a simulated 80 turbines wind farm. A concluding remark of the authors is that whichever the strategy, the average production of energy does not vary. The NOWIcob tool is introduced by Hofmann and Sperstad [34]. The tool uses a multi-component approach, although the authors do not mention the number of modeled components. The example presented shows the importance of transportation vehicles in the oshore context. The authors show that the use of a mothership can increment by more than 3% the availability when compared to the use of a platform. Dalgic et al. [20] also model turbines as a series of subsystems. Each subsystem is modeled by its reliability function derived from historical failures rates and expert judgment. Nevertheless, no information is given about the number of components. The authors also model other uncertain parameters such as the weather (wind and waves). Although the maintenance strategy is not the focus of the work, the authors conclude that remote-monitoring (and thus condition-based maintenance)

can lead to important improvements on the oshore wind farms performance metrics.

Overall, the results of [25,64,62,1] show that under the presence of limited resources, maintenance strategy outputs can signicantly vary. This follows also the conclusion of Van Horenbeek [85] who claims that assuming maintenance operations duration as negligible can conduce to bad decisions.

Therefore, one can conclude that although models are simplication of the real wind farms, several aspects must be considered to compare dierent maintenance strategies. Moreover, most of the studied works consider both the costs and energy related objectives, but only Abdollahzadeh et al. [1] explicitly construct a set of solutions for which costs and energy production change. It is also interesting to see how the scheduling of limited resources can be tackled. Nearly all of the cited works schedule the resources as fast as they can be assigned, ignoring the possible gains (and eects) of optimizing the resources. Only the work of Sahnoun et al. [66] consider a rule to prioritize some of the tasks. These two issues, the multiple objectives and the rules to schedule operational resources, are considered in section 5.6 to evaluate the impact of operational decisions on dierent metrics in a long term horizon.

Maintenance strategies: relation with operational planning

In this section, the maintenance scheduling problem presented in section 5.4, is embedded in a longhorizon simulation model to evaluate how costs and energy production objectives perform through the life cycle of a wind farm project. Nevertheless, since the model and solution approach of section 5.4

requires much time to be practically often solved in a long-horizon evaluation, two actions are taken.

First, the model is simplied to consider only one type of technicians and no extra-hours. Second, it incorporates straightforward rules to schedule the limited resources. Therefore, these rules are used to heuristically solve the scheduling problem instead of using an exact solution method.

Problem Description

Let T to be the planning period for a wind farm. During this time the O&M performs its activities while a set J of turbines produce energy. Each turbine j ∈ J produces P E jt units of energy that depends on the wind speed during a period t ∈ T . Two dierent strategies to perform the maintenance are considered: "'On failure mode"' and "Preventive mode" strategies. On failure mode is a reactive strategy that considers a maintenance tasks if a failure takes place. Besides, the "'Preventive mode"' considers maintenance tasks based on xed time periods. Furthermore, in this mode if a failure takes place then a corrective maintenance task is performed. Besides, whichever the strategy, twice a year the turbines are visited to execute some maintenance on the components not explicitly considered and to change consumables. The dierent types of maintenance are identied by the set U indexed by u.

Additionally, each turbine j is composed by a set of components K indexed by k. Moreover, turbines (and their components) are subject to deterioration process which derive in failures. The degradation process of one component is independent from the others. Turbines components are maintained by a limited set of technicians (P ) who perform the maintenance. Technicians travel to the turbines to perform their assigned maintenance tasks within their shifts.

To perform a maintenance task type u on component k, a number of χ ku technicians are required, with an associated cost c ku . Proactive maintenance is considered as an imperfect repair, that is, a preventive maintenance action does not restore the component to an as good as new condition.

Conversely, a replacement brings the component back to a state as good as new. Furthermore, an execution time β ku characterizes every task. Time β ku corresponds to the time elapsed from stopping the turbine until it is restarted and veried again. It is assumed in this section that β ku is a random variable with known probability distribution. When a task is started by a set of technicians it is not stopped until nished, therefore, the assigned technicians will not work on other tasks until they nish the started activity. Still, at the end of each shift, the technicians stop their work and continue on the next workable shift. Moreover, if the conditions are not safe (high wind speeds) the technicians wait until they become safe again.

Maintenance Strategies are evaluated through several metrics gathered from the whole planning horizon: the average availability (A) of the wind farm, the amount of total energy produced (E ), the total number of failures (N F), and the total costs (T C). The objective is to compare the different strategies, taking as decision variables how and when to perform the schedule of the resources (operational problem). Again, this schedule could be executed with the model presented in section 5.4. Nevertheless, this is computationally intractable due to the length of the planning horizon (and thus the number of operational problems to solve). For this reason, other characteristics such as the technicians skills are not considered in this part. Also, to overcome the computational problem simple priority rules to dene the schedules are devised.

Simulation model

The simulation model is formed by four big modules: the wind turbines, the weather model, the scheduler model and the resources module. A scheme of the modules is presented in gure 5.2.

Wind Turbines

This module represents the turbines, their degradation process, the time between failures, and the power generation. It also allows to keep information as availability and energy production. In this module a set of |J| identical turbines are considered, each of them composed of four components

Modeled components

A turbine is modeled by four key sub-components3 : the rotor, the main bearing, gearbox and generator. Such components have proven to be one of the more important cost inductors for maintenance operations. It is also assumed as reported by [81,22,6,68,1] that the component times between failures can be modeled using a two-parameters Weibull distribution. The parameters reported by Abdollahzadeh et al. [1] for three dierent types of turbines as shown in table 5.4 are used. Since imperfect preventive maintenance is considered the virtual age va jk of each component k in turbine j (∀k ∈ K, j ∈ J) [28]. As time pass the virtual age of each component continues to increase. It is only modied in two cases: if a failure takes place, in which case the component is replaced and the age is restarted (va jk = 0), and in the case of preventive maintenance. In the second case the new va jk is set to va jk = va jk (1 -q). That is, the virtual age is reduced as percentage of the old age. The term q (age reduction ratio) can be component dependent. For simplicity in this work it is assumed that q is the same for every component in every turbine. It is further presumed that q follows a continuous uniform distribution [14] q ∼ U (0.8, 0.95).

Power generation

Equation (5.23) shows the power generated by a turbine j based on the formulas of Karki and Patel [38] given a wind speed v t at time t in the turbine location l. day are known to the O&M operator in a hourly basis. Moreover, next ve days previsions are also known with certitude. At this time the operator also knows which tasks are needed to be performed (except for possible imminent failures).

The following rules are tested to design the tasks scheduling:

O1 : schedule the maintenance tasks as fast as possible when there are available resources O2 : schedule the maintenance tasks to start at times where the wind speeds are low and always in the morning

Resources

The resources considered are technicians working in a shift starting at 8:00h in the morning and nishing at 17:00h from Monday to Friday. Technicians perform the maintenance tasks if wind speeds are below 10 m/s. At higher speeds the technicians need to wait. All other resources as vehicles, spare parts, cranes or any other requirement are always available to perform the maintenance tasks.

Model implementation

The model described in section 5. The model is tested on a virtual instance with the following characteristics. A total of 100 turbines with rated power of 2MW are considered. Each turbine is modeled following the assumptions presented in in section 5.6.2. For the sake of simplicity, turbines are located in a square location of 9 square kilometers in a grid with intervals of 300 meters.

Additionally, consider C k1 and C k2 to be the preventive and corrective costs for a component k ∈ K. Table 5.5 presents this costs and the number of technicians (persons) and the expected time required to carry out the activities. It is assumed as in [1]

that β k1 ∼ N (µ β k1 , σ 2 β k1 ) and β k2 ∼ N (µ β k2 , σ 2 β k2
). The values of σ 2 β k1 and σ 2 β k2 are inferred from a coecient of variation of 0.2 and 0.3 respectively [1]. Table 5.5 values are gathered from [6,81,1,28] and are reported for Type 1 turbine from table 5.4. In addition, the number of required technicians for the recurrent visits (u = 3) is two, with a duration uniformly distributed between two and four hours 4 . It should be noticed that even if most of the data is gathered from the literature, it is not claimed that it to be representative for any particular wind farm.

Simulation model results

To gather the results 1000 simulations are run for each of the strategies and scheduling rules: CM-O1, CM-O2, PM-O1, PM-O2. The "`CM"' strategies are those based on the reactive strategy, meanwhile "`PM"' use a preventive approach. The "`O1"' and "`O2"' parts stand for the rules presented in the sum of the cost of all maintenance tasks performed during the planning period. It is noticeable that PM type strategies present lower costs than CM strategies. These results are expected since the costs of failures are much larger than those due to preventive activities (nearly four times larger, see table 5.5). Moreover, O2 shows higher costs than O1 specially for the PM strategy. Due to that O2 tries to schedule activities more widespread in time, to exploit low wind periods. This eect might be also amplied by selecting mornings to schedule the activities and thus, avoid between days loss of production. Therefore, preventive activities can be executed on larger horizon times (several weeks), augmenting the probability of failures. Figure 5.5, conrms that the costs are increased due to the total number of failures. One should also note that O2 strategies present more variability. This appears to be the result of the vulnerability of O2 subject to longer scheduling periods. It also oers a clue to derive new rules in which the uncertainties (appearance of failures) are considered so the plan of preventive maintenance follow an intelligent order. This analysis lead us to another characteristic to include in the model, the condition base maintenance strategies. Due to the imperfect maintenance tasks, each component of each turbine presents dierent degradation levels and therefore, dierent probabilities of stopping working. Such information can be integrated so new improved rules exploit it to make better decisions. Figure 5.6 presents the total energy production for all maintenance strategies at the end of the simulation. The O2 policy seems to oer the higher production. The dierence in the total energy production between PM-O1 and PM-O2 gives an improvement of nearly 1.7% for the latter. This is further conrmed with the fact that PM-O2 shows the highest availability (see gure 5.7). The explanation of this dierences is again, the way in which O2 plans the activities. The avoidance of scheduling tasks, for example on a Friday afternoon, imply that the turbine will be operative during the weekend, and the tasks will begin preferably on the Monday.

To further explore the dierences between the tested strategies and scheduling rules, a life cycle analysis is conducted. The same four metrics (costs, energy production, availability, and number of failures) are gathered at the end of each year during the simulation of each strategy. Around the sixth and seventh year of operation, the number of failures start to increase in strategies that do not use preventive maintenance, this puts in evidence the eect of preventive strategies in the long term. This fact coincides roughly with mean time to failure of the components (see table 5.4), time at which failures probability becomes more important. The reduction in the number of failures is reected one or two years on the costs and the variation between strategies becomes more prominent at the tenth and eleventh year. The lag between failures and costs can be explained by the preventive costs matching the corrective costs of CM strategies. Nevertheless, after one to two years PM strategies take the edge and become less costly. Other important aspect to remark is that after the tenth year the total costs seems to be more stable for the PM strategy while the CM strategy exhibit a steady increase.

Figure 5.8 shows the behavior of the wind farm availability. This metric is mainly aected by the scheduling rule, particularly for the preventive strategy PM-O1. Such strategy greatly diminishes the availability metric around the sixth year of its life cycle. This is due to the preventive behavior of that strategy since turbines are stopped to execute the maintenance. Nevertheless, after a minimum value of 95% in the eight year, the availability increases for PM-O1 since very few failures occur.

Meanwhile, strategies based on O2 rule show a more stable behavior, decreasing only by 1% during the whole horizon-length. Even though O2 allows a bigger availability this also comes with increasing costs for the PM strategy.

The presented results conrm the idea that maintenance strategies are importantly dependent on how operational activities are planed and executed. Although the results were gathered for a virtual example in which resources were highly constrained5 the conclusions might just hold for more realistic environments. Actually, the reduced number of technicians compensates with the lack of lead times result of unavailable spare parts or need for special cranes.

Conclusions

In this chapter Maintenance Planning problems in the wind farm context are presented. The rst part tackled the operational level problem, in which the resources demanded to perform maintenance tasks are assigned and scheduled. In the second part the integration between the operational and strategic level is shown and the eects of dierent maintenance strategies on the operational cost and energy production is studied.

To solve the maintenance scheduling problem, a bi-criteria mixed integer program is proposed and solved within an epsilon constraint method for the scheduling problem. The results show that O&M costs and the energy produced objectives are in conict an have signicant consequences on each other. Moreover, these results imply that the resources will be assigned depending on the objective function. Thus, the eects on the long-term of using schedules based on one type of objective need to be addressed.

Using a Monte Carlo simulation within an Discrete Events framework we compared two classical strategies widely used in the literature and in real systems: corrective and preventive strategies are compared. Moreover, two scheduling rules at the operational level are tested within each strategy.

The results shown that the preventive strategy impacts in a positive way the costs and number of failures comparing with the corrective strategy: the cost are reduced on nearly a quarter and failures by half. Meanwhile, O2 scheduling rule benets the energy production and availability over rule O1

showing increments of nearly 2% on both metrics at the expense of a 9% increment in costs.

Undergoing work is concentrated on integrating the complete model presented in section 5. as re-optimization can be assessed against the classical recourse policy used in our methods. Currently, our results can serve as the baseline for evaluating the tradeo between algorithms eciency and the use of more complex policies to face uncertainties.

Second, a VRP with stochastic travel and service times and hard time windows was introduced since the majority of related works considered soft time windows, or coupled soft and hard time windows.

This problem also considers continuous random variables with dierent probability distributions which is fairly uncommon in the related literature. A parallel Memetic Algorithm framework is used to tackle the problem. Although the base of each thread of the algorithm is based on the same principles, each one undertakes the uncertainties in a dierent way. Some of the populations are created to solve deterministic problems for which the uncertain parameters are replaced by the mean values or a percentile of the considered probability distributions. To the best of our knowledge, this is the rst time a SVRP is solved using such a solution approach. The results conrm a signicant improvement on the eciency due to the multi-populations parallel algorithm. Moreover, the proposed algorithm is compared to other works published recently, and the results proved that it is very competitive.

Besides this, it can easily be extended to other SVRPs which make the designed parallel framework a powerful and exible tool.

The frameworks of the methods designed to solve the previous problems allowed us to state that good solutions in a deterministic related problem can be used as departure points for good stochastic solutions. Nevertheless, they require to be further improved by integrating adequately the uncertainties. In our approach, this strategy led to improvements on the algorithms performances. In particular, it increased the probability of nding near to best known solutions in shorter running times.

In the wind farm turbine maintenance context, a multi-objective optimization problem dealing with scheduling of maintenance resources is studied. The two objectives taken into account in our problem are the costs and energy production. A bi-criteria mixed integer program is proposed and solved within an epsilon constraint method. To the best of our knowledge, this is the rst time that this kind of approaches is conducted in the wind farm operational decision level context. The computational results suggest that with our formulation small size instances can be solved by providing the exact Pareto-front. Nevertheless, they also show that augmentation of the costs, mainly driven by overtime, can signicantly aect the amount of produced energy. Therefore, we extended this work to evaluate how these objectives (costs and energy production) behave through the life cycle of a wind farm project under dierent maintenance strategies. To do so, as is common in the literature, we devised a Monte Carlo simulation approach based on discrete events simulations. The appearance of failures and the time required to carry out maintenance activities are considered as stochastic.

Unlike what is usually used in other works, we derived a policy to schedule maintenance tasks in short term horizon, considering the wind speed for the subsequent days as well as the energy production.

This allows to emulate good solutions to the operational scheduling problem. The results show an important eect in the life cycle maintenance costs and produced energy. Although other works have employed simulation to compare maintenance strategies, most of them rely on assigning the resources to execute the maintenance tasks as fast as they are available. Further research must be devoted to the introduction of more information to schedule the resources, such as the system status, components history of failures, etc.

Further work on both SVRPs and wind farm maintenance planning should focus on the independence assumption of the random variables. Concerning stochastic demands, the interdependent scenario can arise when demands of close geographical assets are aected by localized events. Moreover, one can expect that variables representing travel times are highly correlated due to factors such as trac jams. Similarly, in the wind farms context, the degradation processes of the turbines components are usually assumed to be independent. This assumption is rather questionable.

Other interesting topics of research are related to the evaluation of the distributions used to model uncertain parameters. For example, the widespread devices carrying Geographical Positions Systems (GPS) can be used to conduct several analysis of travel times in dierent areas. This will allow to evaluate the accuracy and pertinence of using some distributions. In this sense, more experimental studies are needed to conrm the importance of SVRP solution in real applications. In theoretical works the Value of the Stochastic Solution (VSS) has proven to be considerable. Our studies showed that the VSS increases with the size of the problems and also when they are more constrained.

We expect to see an increase in the number of metaheuristic approaches to deal with SVRPs.

The experience of deterministic combinatorial problems must be capitalized and benet to stochastic variants. This must be accompanied with new strategies to incorporate the parameters variability.

Currently, one of the main problems in such methods remains the time overload generated by probability calculation to asses either a recourse action or a probabilistic constraint. This is particularly important in local search procedures which are recognized as time-consuming components. Besides, CHAPTER 6. CONCLUSIONS branch-and-price methods can be used to increment the amount of exact methods on SVRPs. Approximations of the recourse costs based on sampling procedures distributions (simulation) can be used to solve the pricing problems. Additionally, more recourse actions are expected to appear but they will likely be more complex.

To enable comparisons among dierent recourses, standardized benchmarks are needed to be set so that fair competition should contribute to improve the results and highlight new properties on these problems. Finally, although it was not the focus of this dissertation we expect to see an increase in the number of articles dealing with robust optimization as an important methodology to handle uncertainties when the parameters probability distributions are unknown. The trends of the literature on SVRPs and Maintenance Planning make us believe that it will continue to increase in the following years. We hope that the contributions of this thesis and the research clues identied through it will help future developments. Enn, nous pensons que l'optimisation robuste connaitra un développement important pour complémenter et enrichir ce domaine de recherche quand les lois de probabilités sont inconnues. Nous espérons que les pistes de recherche proposées et que les travaux présentés dans cette thèse serviront pour de noubreux développements futurs.

2. 1 .

 1 DETERMINISTIC VRP the service starts at customer i ∈ V c by vehicle l. Model M2VRPTW minimizes (2.7) the total cost assuming a xed eet of vehicles. M 2V RP T W : min l∈L i,j∈V x ijl c ij (2.7) l∈L j∈V

Figure 3 .

 3 1 presents an example of the Split procedure on a VRPSD instance composed by ve customers (V c = {1, 2, 3, 4, 5}) and where vehicles have a limited capacity Q = 50. Moreover, it is assumed that all clients have a stochastic demand which follows a Poisson distribution with mean 20 units, that is q i ∼ P oisson (20) ∀i ∈ V c and costs (c ij ) are equal to the distances. Part (a) of gure 3.1 represents a permutation of clients that is going to be decoded. Part (b) gives the necessary distances information to decode the permutation. Part (c) stands for the Split auxiliary graph, note that each arc has an associated cost (calculated by means of equation (4.2)), furthermore, the dashed arcs are the ones used in the solution of the shortest path. Take for example arc (2, 1), it represents the route r = {0, 4, 1, 0} and has an expected cost of E [C r ] = 88.68, the cost can be further divided into deterministic cost 85 units, and expected recourse cost 3.68 = 88.68 -85.

Figure 3 .

 3 1 part (d)

Figure 3 . 1 :

 31 Figure 3.1: Split example for the VRPSD.

Figure 3 . 2 :

 32 Figure 3.2: Broken pairs distance example.

Figure 3 .

 3 Figure 3.2 shows two solutions for which the number of broken pairs is to be calculated. Solution 1 serves as base for the comparison while solution 2 is a candidate to enter into the population. The distance measure is two since the arcs 2 -3 and 4 -5 present in solution 2 are not present in solution 1.

  3.4. NUMERICAL RESULTS

Figure 3 . 3 :

 33 Figure 3.3: Avg. gap of MA+GRASP against the lling coecient

Figure 3 .

 3 4 presents the graph of both variables. The graph uses the information of all instances but Tai385 which also is marked as an extreme point when compared to the rest of the testbed5 . A clearer pattern of a positive relation between the variables is devised for this relation. Nevertheless, a further inspection on gure 3.5 shows that the number of vehicles is highly correlated with the number of nodes per instance. Moreover, this same gure shows a very clear pattern in the relation of the average gap against the number of nodes in the instance. The results are not surprising as far as this parameter has an important inuence on the speed of the local searches of the MA+GRASP and is of importance in the framework of the MA. Consequently, new methods relying on local search procedures should aim to minimize the complexity of the task.

Figure 3 . 6 :

 36 Figure 3.6: MTTT Plots for the dierent MA methods -5% and 1%

Figure 3 . 7 :

 37 Figure 3.7: MTTT Plots for the dierent MA methods -0.5% and 0%

  terms of quality and eciency. Moreover a new testbed with nearly 2.4 more customers on average per instance is derived. The results on this new benchmark conrm the pertinence of the MA+GRASP which can eciently solve instances with as much as 385 customers. The new instances, which are closer to real life size problems, and the presented results can be used for future comparison. Research currently underway includes the evaluation of new recourse actions, the extension of the problem to consider other stochastic parameters in order to solve problems closer to reality and the introduction of correlation among the demands. Finally, undergoing work also search the means to combine heuristic methods with exact algorithms to tackle large instances of dierent SVRPs.ContributionsPreliminary results of this chapter were presented at CIE45 conference: Gutierrez, A., Dieulle, L., Labadie, N., Velasco, N. (2015) A Memetic Algorithm for the Vehicle Routing Problem with Stochastic Demands In Proceedings of the 45th International Conference on Computers & Industrial Engineering CIE45 Metz, France, 2830 October, 2015. An article version of this chapter has been published in the Computers & Operations Research journal. Please cite it as follows: Gutierrez, A., Dieulle, L., Labadie, N., Velasco, N. (2018). A Hybrid metaheuristic algorithm for the vehicle routing problem with stochastic demands. Computers & Operations Research, 99, 135-147. https://doi.org/10.1016/j.cor.2018.06.012. The Vehicle Routing Problem (VRP) is one of the most studied combinatorial problems; its importance is inherent in the number of real applications where it arises, as well as its signicant theoretical contributions to combinatorial optimization eld. First proposed by Dantzig and Ramser [10], many variants of the problem has been considered, e.g. Vehicle Routing Problem with Time Windows (VRPTW), Distance Constrained Vehicle Routing Problem (DVRP), Orienteering Problem, amongothers. Each of the proposed variants aims to address the problem in a more realistic way to nd more reliable solutions in real world applications. Nevertheless, one of the most important drawbacks in classical models is that available information can be limited because many of the input parameters are uncertain, such as travel times, service times, customers demands, or customer presence. Neglecting the lack of complete information can conduce to, for example, suboptimal solutions, so implying higher operational costs or customer dissatisfaction. Stochastic Vehicle Routing Problems (SVRP) arise when parameters are modeled as random variables, making models closer to reality but harder to solve than their deterministic counterpart. SVRP are usually solved by means of stochastic programming. Two approaches are often used to model and solve a stochastic optimization problem: Chance Constrained Programming (CCP) and Stochastic Programming with Recourse (SPR). CCP aims to solve the problem by bounding the probability of constraints violations to a threshold. SPR uses recourse which are actions to recover the feasibility of the solution when failures occur. The expected costs related to these actions are taken into account in the objective function. Both approaches (CCP and SPR) rely on a two-stage approach: in the rst stage a priori solution is created and then at the second stage the parameters are revealed. However, since CCP does not consider the cost associated with failures, the quality of the solutions might be inferior to those provided by SPR models.

4. 1 .

 1 INTRODUCTIONvehicle routing problem with hard time windows and stochastic travel and service times (SVRPTW)

Figure 4 .c rj rj+1 + k j=1 PFigure 4 . 1 :

 4j=141 Figure 4.1 presents a graphical example of both arrival times and initial service times. The example presents a technician arriving at customer i (a), performing its service if he arrives before l i (b), or continuing its route, and then going to node j (e). Furthermore, gure 4.1 part (a) presents vertical lines representing e i and l i at time 500 and 540. The parts (c) and (d) of the gure show the density function of service time for customer i and the travel time from i to j, respectively. The arrival time at node j presented in part (e) also shows the eects of time windows. Since time windows are hard, the use of convolution properties are precluded, making it harder to properly model the arrival times. Thus, AT rj is usually approximated by a random variable AT rj ∀j ∈ V with a known distribution allowing tractable computations. The route cost equation can then be rewritten in terms of the approximated arrival times.

3 .

 3 The MA starts by creating a population (line 1 of Algorithm 3), which is a set of encoded solutions to the problem at hand. The MA iterates until a stopping criteria is met (line 12 Algorithm 3). At each iteration a new chromosome (C) is processed. The chromosome C can suer a mutation process (line 4 Algorithm 3) or it can be decoded to perform a local search (lines 5 and 6 of Algorithm 3).Then at line 11 of Algorithm 3 the population is updated considering C. Each of the components of the MA and its procedures is now explained in detail.

  an encoded solution to the tackled problem. The MA chromosomes that we use are coded as permutations of the n clients. The decoded version of a chromosome is composed of a set of routes that satises the problem constraints. Henceforth we use the term chromosome to make reference to the permutation representation while the term individual is used to dene its associated decoded solution.

r∈K ζ r ≤ θ ( 4 . 12 )

 412 To overcome this, we rewrite constraint (4.3) as done by Errico et al.[12]. Note that settingζ r = -ln (U r ) and θ = -ln (γ), constraint (4.3) is equivalent to equation (4.12). By considering TIMES Chromosome cost: 75.54+40+60.75 = 176.29 Probability no failures: 0.97x1x0.98 = 95.06%

Figure 4 . 3 :

 43 Figure 4.3: Split example for the VRP with stochastic travel and service times

Figure 4 . 4 :

 44 Figure 4.4: Example of OX crossover.

  i.e. Or-opt movements are limited to movements involving, at most, three customers, and CROSS exchanges use sequences to at most two clients. Also, LS is performed to an individual with probability p ls . Mutations are performed by selecting a random position among the chromosome, then the selected customer and the following are relocated to another random position. Mutations are executed with a probability of p m . CHAPTER 4. VEHICLE ROUTING PROBLEM WITH STOCHASTIC TRAVEL AND SERVICE TIMES 4.4.6 MPMA framework

CHAPTER 4 .

 4 VEHICLE ROUTING PROBLEM WITH STOCHASTIC TRAVEL AND SERVICE TIMES

  to 80 tasks, several modes, skills and locations solved.Besnard et al. [8] Algorithm -Simulation ApproximatedCosts can be reduced delaying the schedule of preventive tasks Stålhane et al.[79] Optimal/Approximated Simulation is used to address a dynamic version of the problem Irawan et al.[35] A MULTI-OBJECTIVE APPROACH TO THE MAINTENANCE SCHEDULING PROBLEM 5.4 A multi-objective approach to the maintenance scheduling problem

Figure 5 . 1 :

 51 Figure 5.1: Approximate Pareto Front for Froger et al. [31] instance 10_2_1_20_B_5

Figure 5 .

 5 Figure 5.1 presents an example of the solutions found. The instance originally contains 10 period times but this value was increased to consider extra time periods. Three periods are dened for

Figure 5 . 2 :

 52 Figure 5.2: Scheme of the simulation modules.

6 . 1

 61 was implemented using a Discrete Events Simulation program designed for this purpose. All runs are conducted on a Dell Latitude E6420 personal computer with Intel Core TM i7-2760QM @2.4 GHz, running Windows 7 Professional 64 bits. The horizon length is set to a 15 years horizon length assuming the wind farms working 365 days of the year and 24 hours per day. A total of 1000 replications are run for each combination of strategy and rule of decision from the scheduler.

5. 6 .

 6 Figure 5.4: Total maintenance costs by strategy/rule for the wind farm simulation.

Figure 5 .

 5 Figure 5.4 presents the box-plot for the total costs measured in monetary units (MU) calculated as

Figure 5 . 8 presents

 58 the values averaged from the 1000 runs across 15 years, and for each one the cumulative value is reported.

Figure 5 . 5 :

 55 Figure 5.5: Total number of failures by strategy/rule for the wind farm simulation.

Figure 5 . 6 :Figure 5 . 7 :

 5657 Figure 5.6: Total produced energy under dierent strategy/rules for the wind farm simulation.

  0

Figure 5 . 8 :

 58 Figure 5.8: Temporal analysis for wind farm metrics under dierent strategies.

1 )(A. 2 )

 12 Pour modéliser le VRPSD nous avons utilisé la programmation stochastique avec recours. Le recours utilisé est le même que celui proposé par Bertsimas[3]. Celui-ci prévoit que lorsqu'un véhicule atteint sa capacité Q, il retourne au dépôt pour déchargement. Ensuite, le véhicule reprend sa tournée depuis le client où la rupture de charge est survenue pour compléter la demande qui n'a pas pu être satisfaite, puis le véhicule continue avec sa tournée. L'objectif du VRPSD avec recours est de minimiser le côut total des tournées y compris ceux relatifs aux recours. L'évaluation du coût d'une tournée est donnée par l'équation A.1. Elle inclut les coûts déterministes et la valeur moyenne du recours (ERC).Pour cette dernière, la demande cumulée jusqu'au client r i est dénie par D ri = i j=1 q rj . Du fait de l'hypothèse d'independence des variables aléatoires relatives aux demandes, l'espérance et la variance de la demande cumulée se calcule parE [D ri ] = i j=1 E [q r l ], V ar [D ri ] = i j=1 V ar [q r l ].La valeur moyenne du recours est donc calculée par l'équation A.2. Il faut remarquer que la otte est supposée illimitée, et que le nombre de véhicules n'a pas d'impact sur la fonction objectif. Ainsi, pour éviter des solutions avec des tournées avec des échecs fréquents, la contrainte E [D r k ] < Q ∀r est aussi rajoutée.ERC ri = 2 • c 0ri • ∞ u=1 P D ri-1 ≤ uQ -P (D ri ≤ uQ)La méthode de résolution proposée pour résoudre le VRPSD est une hybridation entre un algorithme mémétique (MA)[21] et une méthode de type Greedy Randomized Adaptive Search Procedure (GRASP)[10]. L'objectif de la méthode baptisé MA+GRASP est de trouver des solutions de qualité en peu de temps même pour des grandes instances. Pour échapper aux optima locaux, le MA utilise une méthode de re-démarrage basé sur la méthode GRASP. Le but de cette procédure est d'accroître l'exploration de l'espace de recherche tout en controlant son impact négatif sur les temps d'exécution. L'algorithme 6 présente un aperçu global du fonctionnement de la méthode. Elle utilise une population P op à taille xe dans laquelle chaque individu représente une solution du problème. L'algorithme exécute les lignes 4 et 19 jusqu'à ce que un critère d'arrêt soit satisfait. À chaque itération, MA+GRASP crée une nouvelle solution en croisant deux solutions présentes dans la population. Ainsi, A.2. UNE MÉTHODE HYBRIDE POUR LES VRP AVEC DEMANDES STOCHASTIQUES de procéder est que, même si les solutions déterministes n'ont pas la meilleure performance dans le contexte stochastique, elles peuvent servir de base pour être améliorées rapidement. Alors, une fois que des bonnes solutions déterministes ont été trouvées, elles sont améliorées avec la recherche locale qui prend en compte les coûts du recours. Il faut remarquer qu'une partie de la population est complétée par des solutions aléatoires pendant le redémarrage. Cette partie a comme but d'augmenter la diversité et donc l'exploration de l'espace de recherche. Notre approche MA+GRASP permet d'obtenir de très bons résultats. Une comparaison avec les travaux de Mendoza et al. [17] basé sur un GRASP, et ceux de Goodson et al. [14] ( présentant un récuit simulé) conrme la supériorité de notre méthode. MA+GRASP arrive à trouver toutes les meilleures solution connues (BKS en anglais). De plus, elle présente un écart moyen de juste 0.004%. Cela veut dire que l'algorithme retrouve des solutions très proches du BKS à chaque exécution et pour chaque instance (10 exécutions par instance). De plus, MA+GRASP utilise moins de dix secondes pour trouver ces valeurs, ce qui permet de montrer que MA+GRASP est très ecace pour résoudre le VRPSD. Pour tester la méthode sur des instances plus diciles, on a proposé un nouvel ensemble d'instances basé sur les ensembles A et P d'Augerat qui n'ont pas été considérés par Christiansen and Lysgaard [4] (22 instances), ainsi que des chiers test de Cristodes, Mingozi et Toth CMT (4 instances), et de Rochat et Taillard (13 instances) 1 . Le nouvel ensemble présente en moyenne 2.4 plus de clients que l'ensemble original de Christiansen and Lysgaard. Les caractéristiques et les BKS du nouvel ensemble sont reportées dans le tableau A.2. On constate comme illustré dans le tableau A.2 que le coût des solutions BKS est principalement composé du coût déterministe (Det) sur les instances testées. Cela explique pourquoi les solutions qui considèrent juste le coût derministent sont de bonnes solutions aussi pour le problème stochastique.Cependant, si seuls les meilleurs solutions du problèmes déterministe (BDS) sont évaluées, une forte augmentation intervient dans les coûts du recours lorsque ce dernier est considéré. On peut également clairement noter que pour des instances de plus grande, les économies générées par les solutions stochastiques sont plus importantes.Finalement, pour analyser l'impact de la méthode GRASP sur les perfomances de la méthode hybride MA+GRASP, un ensemble de tests numériques a été eectué. Deux versions modiées de MA+GRASP ont été utilisées pour résoudre le VRPSD pour le nouvel ensemble d'instances. La première considère le redémarrage mais, ne fait pas appel au GRASP (MA+RANDOM). Dans ce dernier algorithme, tous les individus crées pour remplacer l'ancienne population sont générés aléatoirement.La deuxième variante est une version sans redémarrage (NR-MA). Une série de calculs eectuée sur les trois algorithmes a été réalisée, le but étant d'analyser le temps que chacun prend pour atteindre les valeurs cibles pour chaque instance. Cet ensemble de données a servis pour construire des graphiques de probabilités cumulés comme proposé par Reyes and Ribeiro[22] (graphiques MTT). Les gures A.2 et A.3 montrent les graphiques pour diérentes valeurs cibles en fonction de l'écart à la meilleure solution connue. Il est évident que la méthode MA+GRASP a l'avantage quand il s'agit de trouver des solutions proches du BKS. En eet, MA+GRASP a plus de chance de trouver tous les BKS avec une même quantité de temps en comparaison avec le MA+RANDOM ou NR-MA. Donc, ceci conrme que le redémarrage avec la méthode GRASP, qui utilise de "`bonnes"' solutions déterministes comme base pour le contexte incertain, génère un eet positif sur les résultats trouvés par la méthode.On peut conclure par les résultats obtenus que MA+GRASP a une performance supérieure à celle des autres méthodes de la littérature. En comparaison avec les versions sans redémarrage ou avec redémarrage complètement aléatoire, le MA+GRASP retrouve des solutions de meilleure qualité en moins de temps. Dans la continuité de ce travail, on peut s'intéresser à l'utilisation de nouveaux recours, à la considération d'autres paramètres stochastiques (comme le temps), et l'extensions du problème pour considérer des demandes corrélées.

Figure A. 2 :Figure A. 3 :

 23 Figure A.2 : grapheique MTT pour les versions de MA -5% and 1%

Algorithm 9 i ← i + 1 13: end while 14 :

 9114 Epsilon Constraints adapté pour la maximisation de Z = (Z 1 , -Z 2 )Require: numSteps1: Solutions ← { } 2: {z1 * , z2 . } ← max Z1 3: {z1 * , z2 * } ← max -Z2| (Z1 = z1 * ) 4: Solutions ← Solutions ∪ {z1 * , z2 * } 5: {z1 . , z2 * * } ← max -Z2 6: stepSize ← (z2 * -z 2 * * ) numSteps 7: i ← 1 8: while (z2 * -i • stepSize) > z2 * * do 9: {z1 * , z2 . } ← max Z1| (Z2 ≤ z2 * -i • stepSize)10: {z1 * , z2 * * * } ← max -Z2| (Z1 = z1 * ) 11: Solutions ← Solutions ∪ {z1 * , z2 * * * } 12: Sort and check(Solutions) Un modèle bi-objective est présenté pour résoudre le problème. Les deux objectifs sont dénis par Z 1 et Z 2 et cherchent à maximiser la quantité d'énergie produite tout en minimisant l'ensemble des coûts respectivement. Il existe diérentes façons de résoudre les problèmes avec plusierus objectifs. Dans cette étude une adaptation de la méthode epsilon constraints expliquée dans l'algorithme 9 est utilisée. La méthode itère pendant un nombre prédéni de pas. À chaque itération (lignes 8 à 13) l'algorithme cherche à optimiser chaque objectif (séparément) en utilisant des contraintes pour les valeurs de l'autre objectif. Ceci permet de trouver des solutions Pareto optimales ou des solutions pour lesquelles aucun objectif ne peut être amélioré sans endommager un autre. Les lignes 2 à 6 permettent de limiter la recherche entre les bornes dénies par la quantité maximale d'énergie produite et le coût minimal. L'approche Epsilon Contraintes est testée sur les instances proposés par Froger et al. [11]. Certaines modications ont été introduites pour tenir compte des heures supplémentaires. La gure A.5 présente un exemple de solutions trouvées pour trois périodes (créneaux) par jour et une seule compétence. De plus, l'instance contient 20 tâches et considère une période de cinq jours. Un total de 16 solutions non dominées sont trouvés mais seulement 2 sont prouvées optimales. Cela signie que l'exécution des problèmes linéaires pour optimiser z 1 et z 2 s'est arrêtée à cause de la limite imposée sur le temps d'exécution avant de garantir l'optimalité. Toutefois, les écarts moyen sont très faibles et c'est possible

Figure A. 6 :

 6 Figure A.6 : Coûts par stratégie et par règle d'assignation.

Figure A. 7 :

 7 Figure A.7 : Énergie produite par stratégie et par règle d'assignation.

  2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS

			Solution approaches
			for VRPs
	Exact			Approximate
	methods			methods
	Branch and	Heuristics	Metaheuristics	Matheuristics
	Price
	Branch and	Constructive	Tabu Search
	Cut Branch and Price and Cut	Clarke & Wright Nearest Neighborhood	Simulated Annealing Iterated Local Search
		Two-phase	Ant Colony
				Optimization
		Cluster-first, route second	Route-first, cluster second	Genetic Algorithms
		Figure 2.1: Solution approaches scheme.

Table 2 .

 2 1: VRPs taxonomy based on Pillac et al. [139] -Information evolution and quality

			2.2. VRP UNDER UNCERTAINTIES
			Information Quality
			Deterministic input	Uncertain inputs
	Information evolution	Input known beforehand Input changes over time	Static and deterministic Static and uncertain Dynamic and deterministic Dynamic and uncertain

  Constraint(2.18) is equivalent to constraint(2.4), but the VRPSD version uses the expected value of demands. M4VRPSD extends the M1CVRP but constrains the probability of a route failing at least once to a level β. In this model, it is assumed that V β (S) is the minimum integer value which satises that the subset S ⊂ V has a probability of failure not exceeding β. M4VRPSD minimizes the total travel costs (2.20) while the routes probability of failure constraint is handled by(2.21). M3VRPSD and M4VRPSD present the SPRs and CCPs models. The rst one considers the expected cost of recourse (2.17), while the second one bounds failures probability (2.21).

3) , (2.5) -(2.6) (2.19) 2.2. VRP UNDER UNCERTAINTIES The objective function (2.17) minimizes the total travel costs plus the expected recourse costs.

Table 2 .

 2 After the VRPSD introduction, it has been widely studied. Table

	2: Summary of VRPSD literature

Table 2 .

 2

	3: Summary of VRPST literature

  Ehmke et al. [48] studied the VRPST with stochastic travel times. The opening time windows are considered hard while the time windows closure time windows are soft. No penalization is considered for late services. The CCP formulation is solved by means of a TS. The approach used from Ehmke et al. to estimate arrival and initial service times is quite similar to the one used by Chang et al.

	CHAPTER 2. LITERATURE REVIEW -VEHICLE ROUTING PROBLEMS
	[88]). Two approximate methods are proposed to solve the problem, a TS, and an Adaptive Large
	Neighborhood Search (ALNS).
	More recently, [36]. Thus, arrival times are assumed to be normally distributed and initial service times are modeled
	by truncated (at e i ) normal random variables. A statistical study shows that even if travel times follow
	shifted exponential or shifted Gamma distributions, assuming normality for the arrival times allows a
	reliable estimate of metrics such as the mean, standard deviation and the 95th percentile for the start-
	service times. This approach has the advantage of being easily incorporated into classical methods
	for the VRPTW. Stochastic travel and service times are considered by Miranda and Conceição
	variable, the mean and variance of the arrivals times are calculated in a closed form if customers
	have no service times. If service times are considered, a rst order approximation is used (see also

  . A CCP model is deducted to guarantee a service level on the vehicles capacity and customers time windows. The CCP model is solved by means of a GA and limited results are shown for an instance with more than 200 customers. Zheng and Liu [179] propose a VRPTW with uncertainties in the travel times. A CCP is used to model the problem, considering a service level on customers time windows. Using Monte Carlo simulation embedded within a GA, the authors solve a problem with 18 customers given

triangular fuzzy variables. Xu et al.

[START_REF] Xu | Vehicle routing optimization with soft time windows in a fuzzy random environment[END_REF]

consider the VRPTW with soft time windows. Uncertainty is considered for the tolerable starting time, which is modeled through a random fuzzy variable. That is, services can start earlier than e i or later than l i , but within endurable earliness and lateness times. Therefore, uncertainty aects the time windows. The proposed model aims to minimize the operational costs and maximize the average satisfaction level of customers. This satisfaction level depends on the service starting time, whether it is outside the time window, during endurable times or within the time window. A particle Swarm Optimization (PSO) is used to solve the problem and results are presented for a case study application with 18 customers.

Table 3 .

 3 1: VRPSD Christiansen and Lysgaard[7] Testbed comparison

			Method
	Metric	MA + GRASP	GRASP + HC	SA
	Avg. Gap	< 0.01%	0.02%	0.35%
	Max. Gap	0.14%	0.19%	1.89%
	Avg. Best Gap	0.00%	0.00%	0.04%
	NBKS	40/40	40/40	33/40
	Max. CPU (s)	10.13	102.43	603.80
	Min. CPU (s)	0.69	1.69	9.00
	Avg. CPU (s)	7.39	36.09	268.66
	per route.

  3.4. NUMERICAL RESULTSare better in comparison to all the published methods in the literature. Although times are not scaled due to dierences in programming languages, operating systems, compiler versions and characteristics of the computers, MA+GRASP seems to oer the best performance in terms of execution time. The time metric is specially important in its Max CPU version since the tested instances were small (at most 60 nodes). The MA+GRASP among the three methods oers the best results quality with reasonable running times, suggesting it as an ecient method to deal with real life size problems.Detailed results for MA+GRASP and GRASP+HC are given in table 3.2.To conrm the ability of MA+GRASP to eciently solve the VRPSD, we propose a new testbed of 39 instances composed by the instances of Augerat test sets A and P, which are not considered in[7] (a total of 22 instances), instances from the Cristodes, Mingozi and Toth CMT test set (4 instances) and the Rochat and Taillard instances(13 instances). This new benchmark contains from 22 to 385 customers and from 3 to 46 vehicles necessary to satisfy customers demands. The average number of customers per instance is almost 97. It shall be noticed that some of these have been already

																				3.4. NUMERICAL RESULTS
	3.4.2 New proposed testbed
	Table 3.2: Christiansen and Lysgaard [7] testbed results Table 3.2: Christiansen and Lysgaard [7] testbed results: Continued	MA+GRASP GRASP + HC GRASP + HC MA+GRASP	Instance BKS Avg. Cost Best Cost Avg. Time (s) Avg. Gap Avg. Cost Best Cost Avg. Time (s) Avg. Gap Avg. Cost Best Cost Avg. Time (s) Avg. Gap Instance BKS Avg. Cost Best Cost Avg. Time (s) Avg Gap	A-n32-k5 853.60* 853.60 853.60 5.32 0.00% 853.60 853.60 14.79 0.00% 4.76 0.00% 233.05 233.05 1.86 0.00% P-n20-k2 233.05* 233.05 631.58	A-n33-k5 704.20* 704.20 704.20 4.93 0.00% 704.20 704.20 13.15 0.00% 6.04 0.00% 218.96 218.96 2.53 0.00% P-n21-k2 218.96* 218.96 218.96	A-n33-k6 793.90* 793.90 793.90 4.77 0.00% 793.90 793.90 13.34 0.00% 7.10 0.00% 231.26 231.26 2.83 0.00% P-n22-k2 231.26* 231.26 231.26	A-n34-k5 826.87* 826.87 826.87 5.37 0.00% 826.27 826.87 14.48 0.00% 4.72 0.00% 681.06 681.06 1.18 0.00% P-n22-k8 681.06* 681.06 681.06	A-n36-k5 858.71* 858.71 858.71 7.43 0.00% 858.71 858.71 20.26 0.00% 5.45 0.00% 619.53 619.52 1.21 0.00% P-n23-k8 619.52* 619.53 619.52	A-n37-k5 708.34* 708.34 708.34 7.47 0.0% 708.34 708.34 23.23 0.00% 26.49 0.00% 472.50 472.50 8.72 0.00% P-n40-k5 472.50* 472.50 472.50	A-n37-k6 1030.73* 1030.73 1030.73 6.81 0.00% 1030.86 1030.73 20.14 0.01% 36.25 0.06% 533.83 533.52 10.03 0.00% P-n45-k5 533.52* 533.52 533.52	A-n38-k5 775.13* 775.13 775.13 8.08 0.00% 775.13 775.13 20.11 0.00% 40.40 0.00% 758.76 758.76 10.09 0.03% P-n50-k10 758.76* 759.04 758.76	A-n39-k5 869.18* 869.18 869.18 8.26 0.00% 869.18 869.18 27.89 0.00% 44.17 0.00% 582.37 582.37 10.09 0.00% P-n50-k7 582.37* 582.37 582.37	A-n39-k6 876.60* 876.60 876.60 8.15 0.00% 876.60 876.60 25.33 0.00% 39.7 0.02% 669.33 669.23 10.05 0.01% P-n50-k8 669.23* 669.33 669.23	A-n44-k6 1025.48* 1025.48 1025.48 10.10 0.00% 1025.92 1025.48 33.93 0.04% 52.72 0.00% 809.70 809.70 10.08 0.00% P-n51-k10 809.70* 809.70 809.70	A-n45-k6 1026.73* 1026.73 1026.73 10.02 0.00% 1026.81 1026.73 31.93 0.01% 56.26 0.00% 742.41 742.41 10.08 0.00% P-n55-k10 742.41* 742.41 742.41	A-n45-k7 1264.83* 1264.83 1264.83 10.06 0.00% 1267.05 1264.83 38.47 0.18% 72.10 0.00% 1068.05 1068.05 9.38 0.00% P-n55-k15 1068.05* 1068.05 1068.05	A-n46-k7 1002.22* 1002.22 1002.22 10.10 0.00% 1002.22 1002.22 46.23 0.00% 64.54 0.03% 588.76 588.56 10.11 0.00% P-n55-k7 588.56* 588.56 588.56	A-n48-k7 1187.14* 1187.14 1187.14 10.14 0.00% 1187.32 1187.14 55.05 0.02% 73.36 0.00% 803.60 803.60 10.06 0.01% P-n60-k10 803.60* 803.73 803.60	A-n53-k7 1124.27* 1124.27 1124.27 10.04 0.00% 1124.27 1124.27 80.22 0.00% 85.52 0.00% 1085.49 1085.49 10.11 0.14% P-n60-k15 1085.49* 1087.02 1085.49	A-n54-k7 1287.07* 1287.07 1287.07 10.09 0.00% 1287.41 1287.07 86.17 0.03% BKS column: Marked with * when optimal proven	A-n55-k9 1179.11* 1179.11 1179.11 10.10 0.00% 1179.11 1179.11 66.16 0.00%	A-n60-k9 1529.82 1529.88 1529.82 10.09 0.00% 1529.82 1529.82 102.43 0.02%	E-n22-k4 411.57* 411.57 411.57 1.67 0.00% 411.57 411.57 1.24 0.00%	E-n33-k4 850.27* 850.27 850.27 5.67 0.00% 851.87 850.27 24.66 0.19%	E-n51-k5 552.26 552.26 552.26 10.07 0.00% 552.26 552.26 56.75 0.00%	P-n16-k8 512.82* 512.82 512.82 0.69 0.00% 512.82 512.82 1.69 0.00%	P-n19-k2 224.06* 224.06 224.06 1.72 0.00% 224.06 224.06 3.51 0.00%

Table 3 .

 3 |V | < 100 the time limit is set to ten seconds, it is increased to 80 seconds for instances with |V | <= 200 and 160 seconds for instances with |V | > 200. Demands are assumed to be Poisson distributed, with expected values equal to the deterministic demand values. Moreover, travel costs are calculated as the Euclidean distance between two nodes rounded to the nearest integer.

3 presents a summary of the new testbed instances. For each instance is reported: the number of nodes (|V |), the minimum number of vehicles (Min veh) needed to satisfy the customers demands, the vehicles capacity (Q), the lling coecient (FC) which stands for i∈V

E[qi]

Table 3 .

 3

	3: Summary of new testbed BKS Instance |V | Min veh Q FC Veh Total Det Rec Total Rec VSS BDS A-n61-k9 61 9 100 0.98 10 1144.23 1084 60.23 1215.37 181.37 5.85% A-n62-k8 62 8 100 0.92 9 1430.81 1375 55.81 1533.07 245.07 6.67% A-n63-k10 63 10 100 0.93 11 1459.49 1412 47.49 1581.17 267.17 7.70% A-n63-k9 63 9 100 0.97 10 1847.69 1734 113.69 1991.75 375.75 7.23% A-n64-k9 64 9 100 0.94 10 1569.85 1534 35.85 1699.68 298.68 7.64% A-n65-k9 65 9 100 0.97 10 1313.30 1266 47.30 1421.52 247.52 7.61% A-n69-k9 69 9 100 0.94 10 1259.35 1214 45.35 1339.55 180.55 5.99% A-n80-k10 80 10 100 0.94 11 1987.17 1918 69.17 2109.56 346.56 5.80% E-n23-k3 23 3 4500 0.75 3 569.72 569 0.72 569.72 0.72 0.00% E-n30-k3 30 3 4500 0.94 4 504.55 503 1.55 569.92 32.92 11.47% E-n76-k10 76 10 140 0.97 11 885.11 861 24.11 911.70 81.70 2.92% E-n76-k14 76 14 100 0.97 16 1118.90 1086 32.90 1187.18 166.18 5.75% E-n76-k7 76 7 220 0.89 7 698.95 692 6.95 723.96 41.96 3.46% E-n76-k8 76 8 180 0.95 8 771.23 744 27.23 791.88 56.88 2.61% E-n101-k14 101 14 112 0.93 15 1164.149 1116 48.15 1233.84 166.84 5.65% E-n101-k8 101 8 200 0.91 8 839.47 824 15.47 878.36 63.36 4.43% P-n55-k8 55 7 160 0.93 7 607.71 581 26.71 631.82 43.82 3.82% P-n65-k10 65 10 130 0.94 10 854.06 802 52.06 861.52 69.52 0.87% P-n70-k10 70 10 135 0.97 11 882.01 851 31.01 929.82 102.82 5.14% P-n76-k4 76 4 350 0.97 4 609.54 593 16.54 609.62 16.62 0.01% P-n76-k5 76 5 280 0.97 5 648.11 628 20.11 649.10 22.10 0.15% P-n101-k4 101 4 400 0.91 4 686.81 684 2.81 694.47 13.47 1.10% CMT12 101 10 200 0.91 10 982.80 827 155.81 984.31 164.31 0.15% CMT11 121 7 200 0.98 8 1201.15 1187 14.15 1226.62 189.62 2.08% CMT4 151 12 200 0.93 12 1072.19 1036 36.19 1130.31 114.31 5.02% CMT5 200 16 200 1.00 18 1378.85 1355 23.85 1486.96 209.96 7.27% Tai75a 76 10 1445 0.95 11 1653.85 1649 4.85 1718.39 102.39 3.76% Tai75b 76 9 1679 0.99 10 1353.09 1343 10.09 1438.62 99.62 5.95% Tai75c 76 9 1122 0.94 9 1349.87 1332 17.87 1396.98 110.98 3.37% Tai75d 76 9 1699 0.93 9 1392.86 1391 1.86 1444.39 87.39 3.57% Tai100a 101 11 1409 0.98 12 2107.71 2089 18.71 2388.31 350.31 11.75% Tai100b 101 11 1842 0.96 12 1985.33 1986 29.33 2118.93 181.93 6.31% Tai100c 101 11 2043 0.93 11 1421.66 1413 8.66 1526.73 126.73 6.88% Tai100d 101 11 1297 0.95 12 1602.53 1599 3.53 1748.94 174.94 8.37% Tai150a 151 15 1544 0.94 15 3211.55 3188 23.55 3574.95 530.95 10.10% Tai150b 151 14 1918 0.95 15 2792.39 2789 3.39 3113.52 397.52 10.31% Tai150c 151 14 2021 0.99 15 2406.13 2380 26.13 2546.29 201.29 5.50% Tai150d 151 14 1874 0.97 15 2718.39 2692 26.39 3044.91 406.91 10.64% Tai385 386 46 65 1.00 55 29364.03 28007 1357.03 31360.45 7008.45 6.22%

(Avg. Cost). Additionally, in table

  as well as NR-MA reach 22 of the BKS while MA+RANDOM nds 21 BKS. In addition, the Best gap metric averages only 0.16% with a maximum of 1.21% in instance CMT5 for the MA+GRASP.

		3.4. NUMERICAL RESULTS
	2,00%	2,00%
	1,80%	1,80%
	MA+RANDOM and NR-MA perform similarly with an average best gap of 0.19%. Computational times present a similar behavior, the MA+GRASP uses 40.92 seconds on average to solve each instance, MA+RANDOM uses slightly less time with an average of 40.43 seconds, while the NR-MA reaches the lowest value with 39.81 seconds. Although computational times cannot be directly compared to those of table 4.1, obtained results suggest that our method has competitive times even for larger instances. Table 3.4: VRPSD Proposed testbed results -Costs BKS MA+GRASP MA+RANDOM NR-MA Instance Best Cost Avg. Cost Best Cost Avg Cost Best Cost Avg Cost 607.71 607.71 607.71 607.71 607.71 607.71 P-n65-k10 854.06* 858.30 859.57 859.17 859.71 858.30 859.43 P-n70-k10 882.01* 882.01 883.20 882.01 883.83 882.01 884.37 P-n76-k4 609.54 611.59 616.56 614.04 616.44 614.04 616.48 P-n76-k5 648.11 652.16 652.71 652.16 653.80 652.16 653.02 P-n101-k4 686.81 686.81 687.78 687.67 688.57 686.81 687.91 CMT12 982.80* 982.80 982.80 982.80 982.80 982.80 982.80 CMT11 1201.15 1202.26 1208.66 1202.75 1207.03 1202.98 1208.14 Table 3.5: VRPSD Proposed testbed results -Time, Gaps MA+GRASP MA+RANDOM NR-MA Instance Time(s) Avg. Gap Best Gap Time(s) Avg. Gap Best Gap Time(s) Avg. Gap Best Gap A-n61-k9 10.12 0.04% 0.00% 10.03 0.15% 0.00% 10.02 0.19% 0.00% A-n62-k8 10.05 0.04% 0.00% 10.02 0.05% 0.00% 10.02 0.02% 0.00% A-n63-k10 10.17 0.00% 0.00% 10.02 0.01% 0.00% 10.01 0.02% 0.00% A-n63-k9 10.02 0.22% 0.05% 10.01 0.24% 0.24% 10.01 0.25% 0.19% A-n64-k9 10.36 0.01% 0.00% 10.02 0.01% 0.00% 10.02 0.01% 0.00% A-n65-k9 10.35 0.00% 0.00% 10.02 0.00% 0.00% 10.01 0.00% 0.00% A-n69-k9 10.42 0.04% 0.00% 10.03 0.05% 0.00% 10.02 0.00% 0.00% A-n80-k10 10.17 0.30% 0.00% 10.04 0.41% 0.00% 10.05 0.36% 0.05% E-n23-k3 2.84 0.00% 0.00% 2.21 0.00% 0.00% 2.23 0.00% 0.00% E-n30-k3 5.90 0.00% 0.00% 4.49 0.00% 0.00% 4.18 0.00% 0.00% E-n76-k10 10.15 0.19% 0.00% 10.03 0.44% 0.00% 10.02 0.29% 0.00% E-n76-k14 10.12 0.13% 0.00% 10.02 0.19% 0.00% 10.01 0.06% 0.00% E-n76-k7 10.51 0.05% 0.00% 10.06 0.00% 0.00% 10.02 0.02% 0.00% E-n76-k8 10.16 0.42% 0.09% 10.03 0.31% 0.09% 10.02 0.43% 0.09% E-n101-k8 80.18 0.11% 0.00% 80.18 0.05% 0.00% 80.07 0.08% 0.00% E-n101-k14 80.22 0.30% 0.00% 76.74 0.56% 0.00% 64.30 0.54% 0.00% P-n55-k8 10.14 0.00% 0.00% 10.01 0.00% 0.00% 10.01 0.00% 0.00% P-n65-k10 10.20 0.64% 0.50% 10.01 0.66% 0.60% 10.01 0.63% 0.50% P-n70-k10 10.29 0.13% 0.00% 10.02 0.21% 0.00% 10.01 0.27% 0.00% P-n76-k4 10.22 1.15% 0.34% 10.04 1.13% 0.74% 10.05 1.14% 0.74% P-n76-k5 10.20 0.71% 0.62% 10.07 0.88% 0.62% 10.03 0.76% 0.62% P-n101-k4 80.59 0.14% 0.00% 80.16 0.26% 0.12% 80.19 0.16% 0.00% CMT12 80.69 0.00% 0.00% 79.38 0.00% 0.00% 71.74 0.00% 0.00% CMT11 80.16 0.63% 0.09% 80.22 0.49% 0.13% 80.25 0.58% 0.15% MA+GRASP MA+RANDOM CMT4 80.65 CMT5 80.53 Tai75a 10.24 Tai75b 10.04 Tai75c 10.18 Tai75d 10.19 Tai100a 80.24 Tai100b 80.92 Tai100c 80.41 Tai100d 80.62 Tai150a 80.19 Tai150b 80.93 Tai150c 81.15 Tai150d 81.09 Tai385 164.64 Average 40.93 Max 164.64 Min 2.84 Instance Time(s) Avg. Gap Best Gap Time(s) Avg. Gap Best Gap 1.18% 0.43% 80.39 1.24% 0.98% 1.86% 1.21% 80.54 1.57% 0.90% 0.03% 0.00% 10.02 0.03% 0.00% 0.20% 0.00% 10.02 0.21% 0.00% 0.35% 0.34% 10.03 0.35% 0.34% 0.02% 0.00% 10.03 0.02% 0.00% 0.24% 0.08% 79.48 0.65% 0.27% 0.29% 0.17% 77.57 0.28% 0.17% 0.08% 0.06% 78.71 0.12% 0.08% 0.04% 0.00% 80.08 0.01% 0.00% 1.42% 0.98% 80.27 1.42% 0.43% 0.64% 0.07% 80.16 0.53% 0.00% 0.41% 0.17% 80.16 0.73% 0.06% 1.36% 0.13% 80.31 1.47% 0.48% 1.48% 0.79% 165.52 1.58% 0.91% 0.380% 0.16% 40.44 0.418% 0.18% 1.86% 1.21% 16552 1.57% 0.98% 0.00% 0.00% 2.21 0.00% 0.00% Table 3.5: VRPSD Proposed testbed results Time, Gaps: Continued NR-MA 80.24 1.49% 0.97% 80.46 1.53% 0.73% 10.03 0.00% 0.00% 10.25 0.32% 0.00% 10.02 0.56% 0.34% 10.04 0.03% 0.02% 78.00 0.61% 0.16% 75.76 0.24% 0.17% 79.99 0.10% 0.00% 80.06 0.00% 0.00% 80.19 1.36% 0.88% 80.19 0.41% 0.00% 80.25 0.73% 0.06% 80.20 2.19% 0.73% 163.91 1.58% 0.74% 39.82 0.434% 0.18% 163.91 2.19% 0.97% 2.23 0.00% 0.00% 0,00% 0,20% 0,40% 0,60% 0,80% 1,00% 1,20% 1,40% 1,60% 0,88 0,90 0,92 0,94 0,96 0,98 1,00 Avg. 1,60% Gap -MA+GRASP 0,00% 0,20% 0,40% 0,60% 0,80% 1,00% 1,20% 1,40% 0,65 0,70 0,75 0,80 0,85 0,90 0,95 1,00 Avg. Gap -MA+GRASP Time(s) Avg. Gap Best Gap Filling capacity Filling capacity* -BKS number of vehicles

Table 3 .

 3

		6: Two-way test values for Friedman method
			Metric
			Average Cost	Best Cost
		|R 1 -R 2 |	15.0	8.0
		Table 3.7: Friedman test ranks
		MA+GRASP	MA+RANDOM
	Instance	Avg. Cost rank	Best Cost rank	Avg. Cost rank	Best Cost rank

  tackles the estimation of arrival times using phase-type (PH) distribution for the Distance-Constrained Capacitated VRP with Stochastic Travel and Service Times. Although no time windows are considered, the authors conclude that PH distributions can be used to handle them.[34] focus on the VRP with time windows and stochastic travel and service times modeled as a CCP and solved using an Iterative Local Search procedure mixed with a discrete approximation on arrival times. Travel and service times are considered as normally distributed. Results were presented for seven Solomon[41] instances with 100 customers.

The purpose of this study is:

(1)

to tackle the SVRP with hard time windows, instead of soft time windows as is frequently executed in the literature;

(2)

to propose a recursive approach to estimate mean and variance of arrival times, including the eect of late arrivals at previous customers; (3) given the fact that arrival times probability distribution are in general unknown, and in most cases expensive to calculate, a log-normal approximation is selected among other probability distributions; and (4) to propose a Multi-Population Memetic Algorithm (MPMA) exploiting dierent characteristics in each

Table 4 .

 4

	Instance Family	1: Basic experiment information per family of instances Num. routes Total nodes Avg. customers per route % Fitted log-normal % Fitted Gamma % Fitted Normal
	C1 C2 R1 R2 RC1 RC2	18 16 24 22 16 16	195 460 219 363 143 241	9.38 27.75 8.13 15.5 7.94 14.06	35.03% 23.42% 46.67% 40.47% 28.35% 30.67%	36.16% 21.17% 40.00% 34.31% 18.11% 25.33%	15.82% 10.81% 18.97% 23.17% 9.45% 16.44%
	Table 4.2: Arrival times average absolute gaps between simulated values and log-normal approximation
	Notable percentiles Family Mean SD 90% 95% 99% Mean Min. Max. 1 st to 99 th percentiles Instance
	C1 C2 R1 R2 RC1 RC2	0.10% 0.06% 0.07% 0.03% 0.08% 0.03%	1.41% 2.92% 0.76% 0.73% 1.57% 0.77%	0.21% 0.31% 0.70% 0.19% 0.36% 0.78% 0.13% 0.28% 0.73% 0.06% 0.12% 0.29% 0.25% 0.58% 1.45% 0.08% 0.17% 0.41%	0.26% 0.15% 0.97% 0.21% 0.10% 0.78% 0.20% 0.11% 0.73% 0.08% 0.04% 0.29% 0.44% 0.17% 1.74% 0.11% 0.04% 0.42%

only for the best three distributions. A total of 112 routes were used in the analysis accounting for a total of 1621 evaluated nodes. The number of routes diers from family to family since the number of instances in each family is dierent (except for C2, RC1, and RC2 with 8 instances). Instances of type 2 have a larger number of clients per route as expected because they have larger time windows. Table

4

.1 presents the percentage of nodes that tted log-normal, Gamma and normal distributions

Table 4

 4

	CHAPTER 4. VEHICLE ROUTING PROBLEM WITH STOCHASTIC TRAVEL AND SERVICE
	TIMES
			6,00%
		95 percentile Mean absolute error	1,00% 2,00% 3,00% 4,00% 5,00%					30,00% 25,00% 20,00% 15,00% 10,00% 5,00% 0,00%
			0,00%
			0,00%	5,00%	10,00%	15,00%	20,00%	25,00%	30,00%
					% of noise Travel Times
	Figure 4.2: Average absolute error of 95% percentile with service and travel time noises.
	the fact that larger time windows prevent the truncation of the arrival times.
	Table 4.3: Arrival times mean and standard deviation absolute gaps between simulated values and
	three approximations Instance Family Mean Log-Normal SD		Normal Mean	SD	Gamma Mean	SD
	C1 C2 R1 R2 RC1 RC2	0.10% 0.06% 0.07% 0.03% 0.08% 0.03%	1.41% 2.92% 0.76% 0.73% 1.57% 0.77%		0.22% 0.15% 0.11% 0.03% 0.12% 0.03%	2.96% 4.79% 1.99% 0.90% 3.11% 1.21%	0.22% 0.15% 0.11% 0.03% 0.11% 0.03%	3.09% 4.51% 1.80% 0.84% 2.40% 1.06%

.2 also reects the accuracy of the estimations, with a mean gap across the families and from the 1 st to 99 th percentile inferior to half a percent. Again, RC1 shows the highest gaps in terms of mean, minimum and maximum average gap when compared with the rest of the families. Type 2 instances (C2, R2, RC2) look to have lower gaps than their counterparts (C1, R1, RC1) which might be explained by

  Algorithm 4 the MAs are created. Then, each MA starts to solve its problem (in parallel) in lines three to six of Algorithm 4. While the MAs are working they cooperate by letting chromosomes be copied from one MA to another (line 8 of Algorithm 4). The aim of communication is to transfer valuable information from one MA to another and to enhance diversication.Once all the MAs nish, the MPMA continue its procedure. Since the MAs solutions are created through the approximation of arrival times (see section 4.4.7), they may not be feasible in the uncertain environment. Therefore, to guarantee a feasible solution we proceed as follows. The best ve chromosomes of each MA are extracted (lines 10 to 12 of Algorithm 4) and decoded into individuals using Split (lines 14 to 17 of Algorithm 4). Monte Carlo simulation is used to check the feasibility of the routes, the probability of having no failures, and to estimate their costs.

	Algorithm 4 MPMA
	9: end while
	10: for all f = 1...F do
	11:	Best ← Best ∪ Get best Chromosomes(M A f )
	12: end for
	13: P ool ← {}
	14: for all C ∈ Best do
	15:	Split(C, α, β, γ, numSimulations) // Using Monte Carlo Simulation
	16:	P ool ← P ool ∪ Simulated routes during Split
	17: end for
	18: return SolveSetPartitioning(P ool, γ)

Require: α, β, γ, τ, numSimulations, runT ime, F

1: Create(M A f ) ∀f = 1...

F 2: initialT ime ← CurrentTime 3: In parallel : // The MAs start and run in parallel 4: for all f = 1...F do 5: run(M A f (α, β, γ, τ, runT ime)) // See algorithm 3 6: end for 7: while CurrentTime -initialT ime ≤ runT ime do 8: Communicate M A f ∀f = 1...F

Table 4 .

 4 4: Best solutions found by MPMA for C type instances the average number of vehicles (Avg. # Veh.), the percentage of times the best solution is found (% Times BSF), and the average from 10 runs. The average time stands for the time expended by MPMA in one run of one instance. It comprises from the start time of the algorithm until the best solution is retrieved, i.e. it does not take into account the total computational time generated by parallel threads. To establish if two solutions are equal, we simulate the best solution found at each run of the MPMA, which allows us to construct a condence interval of 95% on the total costs 2 .

	Best Solution -50 Customers Instance ERC # Veh. DC RC ERC # Veh. DC RC Best Solution -100 Customers C101 9741.22 9.00 642.90 98.32 21519.83 20.00 1497.54 22.29 C102 8635.27 8.00 590.81 44.46 18391.69 17.00 1330.67 61.02 C103 7534.46 7.00 514.58 19.88 15227.97 14.00 1193.63 34.33 C104 5382.54 5.00 369.45 13.10 11986.00 11.00 961.25 24.75 C105 7555.61 7.00 522.88 32.73 16292.06 15.00 1199.72 92.35 C106 8609.68 8.00 564.34 45.34 16259.28 15.00 1212.13 47.15 C107 6506.86 6.00 484.18 22.68 14171.11 13.00 1095.82 75.29 C108 6467.17 6.00 440.31 26.86 13069.71 12.00 956.64 113.07 C109 5369.22 5.00 364.11 5.11 10846.15 10.00 836.91 9.24 C201 3432.76 3.00 417.42 15.34 5821.75 5.00 717.08 104.67 C202 2462.64 2.00 433.43 29.20 5755.96 5.00 718.78 37.18 C203 2514.33 2.00 501.01 13.32 4715.41 4.00 679.56 35.85 C204 2384.27 2.00 372.21 12.06 4667.72 4.00 654.73 12.99 C205 2464.49 2.00 460.92 3.57 4688.49 4.00 671.74 16.75 C206 2460.33 2.00 456.51 3.82 4666.89 4.00 656.51 10.38 C207 2449.09 2.00 448.27 0.82 4655.02 4.00 649.24 5.78 C208 2366.15 2.00 364.61 1.54 4632.91 4.00 625.22 7.69
	Total	86336.10 78.00 7947.94 388.16 177367.9 161.00 15657.17 710.78
	Table 4.5: Best solutions found by MPMA for R type instances
	Best Solution -50 Customers Instance ERC # Veh. DC RC ERC # Veh. DC RC Best Solution -100 Customers R101 16232.47 15.00 1191.51 40.96 28960.67 27.00 1934.28 26.38 R102 14082.91 13.00 1049.13 33.78 25744.05 24.00 1725.49 18.56 R103 10890.01 10.00 867.52 22.49 19409.84 18.00 1371.98 37.86 R104 7734.53 7.00 722.66 11.88 14154.45 13.00 1128.20 26.25 R105 12007.66 11.00 986.34 21.32 19602.65 18.00 1527.25 75.40 R106 9898.11 9.00 867.16 30.95 17420.11 16.00 1369.40 50.71 R107 8795.05 8.00 775.87 19.18 14258.45 13.00 1194.59 63.86 R108 6659.90 6.00 631.58 28.32 12068.34 11.00 1040.52 27.82 R109 9863.25 9.00 849.35 13.90 15308.71 14.00 1270.43 38.28 R110 8779.34 8.00 762.32 17.02 14239.67 13.00 1193.42 46.25 R111 8751.49 8.00 745.91 5.58 14184.60 13.00 1151.82 32.78 R112 7678.27 7.00 672.68 5.58 13035.89 12.00 1029.60 6.29 R201 3898.66 3.00 892.30 6.36 6333.61 5.00 1306.96 26.65 R202 3745.80 3.00 739.81 5.99 6090.63 5.00 1079.22 11.41 R203 2726.72 2.00 703.38 23.34 4927.94 4.00 916.71 11.23 R204 2512.02 2.00 511.28 0.74 3792.57 3.00 792.53 0.03 R205 2774.60 2.00 768.96 5.64 6011.41 5.00 1011.05 0.36 R206 3642.05 3.00 642.02 0.04 3981.86 3.00 981.50 0.35 R207 2587.16 2.00 587.16 0.00 3835.66 3.00 835.66 0.00 R208 2493.69 2.00 493.69 0.00 3715.32 3.00 710.28 5.04 R209 2694.99 2.00 693.19 1.80 4926.19 4.00 918.85 7.34 R210 2697.69 2.00 692.05 5.63 4993.53 4.00 992.32 1.21 R211 2561.32 2.00 561.28 0.03 4778.41 4.00 778.20 0.20
	Total	153707.6 136.00 17407.14 300.54 261774.5 235.00 26260.27 514.25
	Table 4.6: Best solutions found by MPMA for RC type instances
	Best Solution -50 Customers Instance ERC # Veh. DC RC ERC # Veh. DC RC Best Solution -100 Customers RC101 12133.54 11.00 1091.58 41.96 20920.01 19.00 1819.27 100.73 RC102 9947.87 9.00 901.19 46.68 18714.96 17.00 1660.41 54.55 RC103 8883.27 8.00 848.98 34.29 15475.45 14.00 1427.97 47.48 RC104 6797.88 6.00 728.56 69.32 13300.61 12.00 1260.09 40.51 RC105 10989.11 10.00 969.13 19.98 18783.43 17.00 1704.47 78.96 RC106 8935.65 8.00 904.74 30.91 16627.30 15.00 1597.60 29.70 RC107 7827.67 7.00 803.86 23.80 15433.93 14.00 1387.96 45.97 RC108 6738.07 6.00 710.82 27.25 14321.21 13.00 1309.87 11.34 RC201 4825.24 4.00 820.78 4.46 7402.04 6.00 1387.02 15.02 RC202 3836.37 3.00 830.57 5.80 6184.45 5.00 1181.73 2.72 RC203 3640.20 3.00 633.32 6.88 5098.22 4.00 1096.44 1.78 RC204 2498.75 2.00 490.95 7.80 4801.08 4.00 800.55 0.53 RC205 4684.67 4.00 676.94 7.73 7245.54 6.00 1229.47 16.07 RC206 2914.94 2.00 856.71 58.23 5229.15 4.00 1195.54 33.60 RC207 3617.36 3.00 614.62 2.74 5071.24 4.00 1067.12 4.11 RC208 2531.56 2.00 525.87 5.68 4879.98 4.00 866.94 13.04
	Total	100802.1 88.00 12408.62 393.52 179488.5 158.00 20992.47 496.11
		Table 4.7: Average performance of MPMA
	Instance Family	50 Customer instances Avg. # Veh. % Times BSF Avg. Time(s)	100 Customer instances Avg. # Veh. % Times BSF Avg. Time(s)
	C1 C2 R1 R2 RC1 RC2	7.04 2.38 9.36 2.63 8.31 3.24	52.3% 57.5% 65.0% 26.4% 61.25% 32.50%	37.61 61.00 35.14 57.94 36.18 51.49	14.33 4.28 16.48 4.51 15.54 5.29	57.7% 70.0% 29.2% 14.54% 38.8% 11.25%	132.14 188.82 128.05 186.06 127.65 170.59

proposed 56 les divided among six families, C1, R1, RC1, C2, R2, RC2 for the VRPTW. Families C1, R1, and RC1 have tight time windows, so routes in these instances tend to have less number of customers. On the other hand C2, R2, and RC2 have larger time windows therefore more clients are visited by each vehicle. Another classication can be done based on the position of customers in the space. Families C1 and C2 have clustered sets of customers, while they are randomly positioned in R1 and R2. RC1 and RC2 have a mix of clustered and randomly located 4.5. NUMERICAL RESULTS

Table 4 .

 4

	8: Comparison of single MA1 to MPMA -100 customer instances MA1 MA1 + MA2 MA1 + MA3 MA1 + MA2 + MA3 # Vehicles Best Solutions Metric 579 575 576 554 # Vehicles Average 616.7 612.7 608.4 577.3 Average Distance Cost 1156.56 1158.88 1138.24 1123.87 Average Recourse Cost 3.86 3.90 3.88 31.16 Avg. Time on pure MA (s) 100 100 100 100 Avg. Time on simulation+set partitioning (s) 35.28 65.06 61.85 54.70 Avg. Total Time (s) 135.28 165.06 161.85 154.80
	Table 4.9: MPMA comparison to Miranda and Conceição [34] ILS Instance Method Avg. # Veh. Avg. DC Avg. Min SL Avg. SL Avg. SLE Max. SLE Time (s)
	R105 R109 C101 C106 RC101 RC106 RC107 Avg.	ILS MPMA 17.30 1470.89 80.55% 97.30% 0.39% 17.67 1615.75 81.66% 99.31% 0.18% ILS 15.00 1488.81 79.53% 95.45% 0.17% MPMA 15.00 1270.27 80.96% 98.00% 0.28% ILS 17.00 2284.77 83.82% 95.52% 0.22% MPMA 17.00 1372.27 80.42% 95.42% 0.64% ILS 14.67 1722.32 81.33% 98.65% 0.18% MPMA 14.00 1141.17 80.57% 96.01% 0.72% ILS 19.67 2012.76 78.32% 100.22% 0.16% MPMA 19.60 1841.32 80.98% 97.53% 0.46% ILS 14.67 1584.63 82.26% 95.57% 0.16% MPMA 15.30 1534.95 81.07% 97.87% 0.40% ILS 14.00 1569.10 83.03% 96.23% 0.16% MPMA 13.70 1368.53 81.46% 98.44% 0.23% ILS 16.10 1754.02 81.42% 97.28% 0.18% MPMA 15.99 1428.49 80.86% 97.20% 0.44%	1.81% 18.40 2.81% 31.82 0.89% 6.36 3.44% 32.21 0.87% 29.49 3.56% 31.77 0.70% 5.56 4.78% 32.25 1.42% 15.18 4.98% 31.77 0.82% 24.43 4.64% 32.08 0.86% 34.99 3.01% 32.40 1.05% 19.20 3.89% 32.04

Local Search (ILS) is used to solve some of the Solomon

[41]

instances with 100 customers. Service times are assumed to be normally distributed with mean equals to the deterministic service time and standard deviation derived from a coecient of variation generated by a uniform law U [0.1; 0.6]. Travel times are also assumed to be normally distributed with mean equals to the euclidean distance between the nodes and standard deviation derived from a coecient of variation generated by a uniform law U [0.1; 0.6]. Some modications were made to our MPMA to compare it with the ILS proposed by

Miranda and Conceição

Table 4 .

 4 10: MPMA comparison to Nguyen et al. [39] TS Instance Method Avg. # Veh. Avg. DC Avg. ET Time (s)

	C1	TS (AD1) TS (AD2) TS (AD3) MPMA	10.33 10.33 10.22 9.89	944.55 929.85 918.61 908.73	2.82 3.43 4.19 5.55	9.20 7.94 8.62 51.41
	C2	TS (AD1) TS (AD2) TS (AD3) MPMA	3.25 3.13 3.00 3.00	654.41 648.60 638.40 606.30	0.00 0.04 0.02 0.08	14.34 9.97 8.52 51.60
	R1	TS (AD1) TS (AD2) TS (AD3) MPMA	15.67 15.83 15.58 15.51	1521.12 1495.02 1447.70 1409.47	30.77 24.02 26.02 19.25	16.10 14.96 18.95 51.76
	R2	TS (AD1) TS (AD2) TS (AD3) MPMA	3.45 3.55 3.45 3.99	1134.96 1050.33 1086.00 991.07	5.33 1.48 1.77 1.96	51.50 42.38 43.90 53.56
	RC1	TS (AD1) TS (AD2) TS (AD3) MPMA	15.25 15.38 15.13 15.86	1736.57 1713.98 1700.99 1630.91	49.79 34.81 41.07 22.55	11.24 10.78 9.46 51.71
	RC2 Avg.	TS (AD1) TS (AD2) TS (AD3) MPMA TS (AD1) TS (AD2) TS (AD3) MPMA	4.00 3.88 4.00 4.56 8.91 8.95 8.82 9.04	1396.98 1296.20 1234.10 1182.94 1241.83 1198.80 1181.67 1131.34	4.95 2.64 3.17 3.86 15.91 11.35 12.92 9.19	33.49 25.82 23.41 52.67 23.48 19.46 19.98 52.16

  Although MPMA needs more time than TS, it nds better results on average, with low expected tardiness and barely incrementing the number of vehicles used. Times are not scaled up since dierences in the operating system, language programming, and hardware, yet, TS procedures work on average 2.5 times faster than MPMA.

Results show that MPMA is very competitive. If compared to TS (AD1) the MPMA increments the number of vehicles by only 1.5%, additionally MPMA can reduce distance costs by almost 10% and more importantly, reduces the expected tardiness by more than 40%. The MPMA presents the best results among the compared methods in both distance (total travel time) and expected 4.6. CONCLUSIONS tardiness.

Table 4 .

 4 11: 50 Customer instances best and average solutions found by MPMA per instance

	Gutierrez, A., Dieulle, L., Labadie, N., Velasco, N. (2016)
	A multi population memetic algorithm for the vehicle routing problem with time windows and
	stochastic travel and service times
	In 8th IFAC Conference on Manufacturing Modelling, Management and Control MIM 2016 Troyes,
	France, 2830 June 2016
	Extensions of the problem at hand were also presented at CLAIO2016 conference:
	Gutierrez, A., Dieulle, L., Labadie, N., Velasco, N. (2016)
	An approximate column generation for the vehicle routing problem with hard time windows and
	stochastic travel and service times. In XVIII CLAIO, the Latin-Iberoamerican Conference on Opera-
	tions Research Santiago de Chile, Chile, 26 October 2016.
	An article version of this chapter has been published in the Computers & Industrial Engineering
	journal. Please cite it as follows:
	Gutierrez, A., Dieulle, L., Labadie, N., Velasco, N. (2018). A multi-population algorithm to solve

the VRP with stochastic service and travel times. Computers & Industrial Engineering, 125, 144-156. https://doi.org/10.1016/j.cie.2018.07.042. MPMA Tables with average results for Solomon instances

Table 4 .

 4 12: 100 Customer instances best solutions and average found by MPMA per instance

	Best Solution

Table 5 .

 5 2: Epsilon Constraints summary results for Froger et al. [31] instances Family Avg. Time (h) Avg. Solutions Avg. E. Solutions Avg. Gap C Avg. Gap E Avg. CNOP Avg. ENOP

	10_2_1_20_A	1.85	12.2	4.4	4.7%	0.7%	2	7.2
	10_2_1_20_B	1.22	10.2	4	0%	0.3%	0	6
	10_2_3_20_A	1.96	10.8	5	5.7%	0.6%	2	4.6
	10_2_3_20_B	2.39	13.4	6.2	9.5%	0.3%	1.2	6
	20_4_1_20_A	2.81	12.8	4.2	16.2%	0.5%	4.6	7.6
	20_4_1_20_B	3.41	15.2	3	12.1%	0.1%	5.2	10.8
	20_4_3_20_A	2.65	8.4	1	28.4%	1.6%	5.0	5.1
	20_4_3_20_B	2.70	7.8	1	34.6%	1.4%	5.4	4.6
	Total	2.37	11.35	3.6	13.9%	0.6%	3.2	6.5

Table 5 .

 5

	3: Summary of literature with complex models for maintenance strategy selection	ONS OFF NT PG WM NC FM CM PM PAM LT LRT EO CO IC
		Reference	Andrawus et al. [6]

  Byon et al. picked the gearbox since the turbines most critical failures are related to this component. The gearbox degradation is modeled through a Hidden Markov Chain with a transition matrix of one week periodicity. This type of model allow the authors to integrate the incomplete information received through sensors output. That is, the real state of the components are hidden and the sensor outputs can be generated by multiple of these states. Using a Discrete Events System Specication (based on the work of Perez et al. [61]) they simulate 100 turbines letting to conclude that condition-based maintenance enables more wind power generation as it reduces the number of failures. Multi-component systems have been preferred to model the turbines. Half of the reviewed works use four components to represent the turbines (see column NC, table

Table 5 .

 5 Furthermore, v 0 represents the cut-in CHAPTER 5. WIND FARMS MAINTENANCE Figure 5.3: An example of a turbine modeled a multi-component system. 5: Costs, technicians, and time requirements per maintenance component based on Abdollahzadeh et al. [1]C k1 C k2 χ k1 χ k2 µ β k1 µ β k2

			Turbine
	C1	C2	C3	C4
	SC1.1	SC2.1	SC2.2	SC 4.1
	SC1.2
				SC 4.1.1
	SC1.3	SC 2.2.1	SC 2.2.2
				SC 4.1.1.1

  ConclusionsThis thesis is dedicated to the Vehicle Routing Problems (VRP) under uncertainties and Maintenance Planning on wind farms. Although most of the works on these elds are devoted to deterministic problems, recently published literature has shown an increasing interest on both research subjects. This growing attention is due to the fact that in real world applications, data are seldom known with precision or with certitude when the decisions have to be taken. Furthermore, the eects of neglecting that data are imperfect when designing maintenance schedule or vehicles routes can have Therefore, we proposed a new data set including 39 instances based on a benchmark originally designed for the deterministic VRP. The average number of customers in this new set doubles the previous classical set, and comprises some instances with more than 200 customers. New results were reported for 29 instances while the remaining ten have been already optimally solved by other authors. The new benchmark and best solutions achieved by our algorithm enable in the future comparison on middle to big size instances for which results are available. In addition, new recourse actions and strategies such

4 within the simulation approach. Moreover, it is intended to incorporate the uncertainties in the scheduling problem, for example due to weather predictions errors. Besides, new rules to schedule resources and new maintenance strategies including opportunistic and condition based need to be tested. Finally,

Chapter 6

important impacts in terms of costs, systems unavailability, or customers' dissatisfaction. Within this context, this thesis makes the following contributions: rst it presents ecient solution approaches based on memetic algorithms with strategies specially tailored for two Stochastic VRPs variants: namely the VRP with stochastic demands and the VRP with random travel and service times in which maintenance tasks on distributed assets must begin within hard time windows. New results are exhibited to these two stochastic VRPs including new best solutions. Moreover, a natural extension of such problem to maintenance planning for wind farms is devised for a deterministic case, and then explored within a simulation approach.

First, a hybridized Memetic Algorithm with a restarting procedure based on a Greedy Randomized Adaptive Procedure Search (GRASP) was proposed to tackle the VRP with stochastic demands. The method is proven to be ecient in a classical available benchmark showing better results than current state-of-the-art methods. During this thesis, the literature review showed a lack of detailed results for bigger instances involving larger number of customers, as can be encountered in real life applications.

  Les activités de transport jouent un rôle très important dans l'économique. En eet, la logistique et les activités de transport ont généré 8% du Produit Intérieur Brut aux États-Unis pendant l'année 2015. Le transport rien qu'à lui même peut représenter jusqu'à 60% des coûts logistiques, ce qui le rend un sujet d'étude important. Sur le plan académique, la communauté de la Recherche Opérationnelle a largement contribué à la résolution de problèmes soulevés dans le domaine du transport, et en particulier les problèmes de tournées de véhicules (VRP en anglais pour Vehicle Routing Problem). et la réoptimisation. L'optimisation à priori est liée au cas statique dans lequel des décisions doivent être prises ici et maintenant bien avant la révélation des réalisations des paramètres stochastiques. C'est le paradigme préféré quand les incertitudes peuvent avoir des conséquences importantes sur la solution, et donc il est préférable de les anticiper. Parmi ces paradigmes à priori, nous citons les problèmes avec recours (SPR en anglais pour Stochastic Programming with INTRODUCTION AUX VRP le chapitre 2 une introduction générale aux problèmes VRPs et une revue des VRPs stochastiques sont présentées. Le chapitre 3 est consacré au VRP avec demandes stochastiques (VRPSD). Pour ce dernier une métaheuristique hybride composée d'un algorithme mémétique et d'une procédure de redémarrage, par une méthode de type Greedy Randomized Adaptive Search Procedure (GRASP), est proposée. Les résultats de cette approche hybride montrent son ecacité en comparaison avec d'autres méthodes publiées dans la littérature. De plus un nouvel ensemble d'instances de grande taille, inspiré d'un benchmark dédié au cas déterministe, est proposé pour servir de base à des comparaisons futures. Le chapitre 4 se concentre sur le VRP avec temps de trajets et de services stochastiques et fenêtres de temps dures. Le modèle proposé intégre l'impact de la violation des fenêtres de temps sous forme de recours, mais impose de garantir des niveaux de services. En outre, pour estimer les temps d'arrivée chez les clients, une approximation par une loi Log-normale est proposée et démontrée ecace par des tests statistiques. Pour résoudre le problème, un algorithme méméthique parrallèle et à populations multiples a été développé. Cette méthode a permis d'obtenir de très bons résultats en comparaison avec ceux disponibles dans la littérature. Le chapitre 5 présente une revue de la littérature consacrée à la planication des activités de maintenance pour un parc d'éoliennes abordé du point de vue décision opérationnelle. Dans ce problème nous considérons deux critères de décision : l'opérateur du parc veut minimiser les coûts de maintenance tandis que l'investisseur du parc veut produire la plus grande quantité d'énergie possible. Un modèle mathématique bi-objectif est utilisé pour modéliser le problème sous forme d'un programme linéaire à variables mixtes. Ce dernier est ensuite résolu par une méthode de type epsilon-contraintes. et le fait que le domaine a été l'origine du développement de diérentes méthodes, exactes et approchées, pour la résolution de problèmes combinatoires.Dans sa version de base, le VRP avec contraintes de capacité (CVRP) a comme objectif de construire un ensemble des tournées de coût minimal qui respectent les contraintes de capacité des véhicules. Le problème est déni sur un graphe complet non orienté G = (V, E). L'ensemble de n÷uds est noté V = {0, 1, . . . , i, . . . , n} et l'ensemble d'arêtes est E = {[i, j] ∀i, j ∈ V | i < j}. Le n÷ud 0 est associé à un sommet particulier appelé dépôt, et le reste des n÷uds V c = V \ {0} représentent les clients. De plus, un ensemble de véhicules ayant la même capacité Q sont disponibles au dépôt. Par ailleurs, chaque arête dans E est associée à un coût non négatif c ij , et chaque client dans V c est associé à une demande q i . La solution du problème est un ensemble de tournées visitant une et une seule avec fenêtres de temps (VRPTW en anglais). Le VRPTW généralise le CVRP par l'ajout des durées de trajet et de service, ainsi que par la présence de fenêtres de temps [a i , b i ] sur le début de service pour tout n÷ud ∀i ∈ V . Comme le CVRP, le VRPTW peut être dénit sur le graphe G. Le VRPTW rajoute un temps de trajet t ij ∀ i, j ∈ V pour chaque arête et un temps de service t i ∀ V c pour chaque client. Les fenêtres de temps sont classées en deux types : dures et souples[8]. La version avec fenêtres de temps dures considère le cas où les services chez les clients doivent impérativement commencer à l'intérieur de la fenêtre de temps. Par conséquent, si un véhicule arrive chez le client avant l'ouverture de la fenêtre de temps, il doit attendre jusqu'à ce moment-là. Dans le cas où le véhicule arrive après la fermeture de la fenêtre de temps, aucun service ne peut être eectué. Dans le VRPTW avec fenêtres souples, les services en dehors des fenêtres sont autorisés mais une pénalité proportionnelle à l'écart entre la date de début de service et la borne de la fenêtre de temps est souvent considérée pour ces évènements. La fenêtre de temps pour le dépôt pose une contrainte sur la date de départ et de retour à ce dernier. Par conséquent, il n'existe pas d'algorithmes de complexité polynomiale qui peut les résoudre pour toute taille de problème. Les diérentes approches utilisées pour obtenir des solutions peuvent être classiés en méthodes exactes et méthodes approchées comme le montre la gure A.1. Les méthodes exactes permettent de trouver la solution optimale, cependant leurs temps d'exécution sont très élevés à partir d'une taille donnée. Ceci est d'ailleurs le principal inconvénient de INTRODUCTION AUX VRP tableau A.1 présente une classication des VRPs, telle que suggérée par Pillac et al. [20] selon la qualité et l'évolution de l'information.

	Annexe A
	Résumé en français
			Solution approaches
			for VRPs
	Exact			Approximate
	methods			methods
	Branch and	Heuristics	Metaheuristics	Matheuristics
	Price
	Branch and	Constructive	Tabu Search
	Cut Branch and Price and Cut	Clarke & Wright Nearest Neighborhood	Simulated Annealing Iterated Local Search
		Two-phase	Ant Colony
				Optimization
		Cluster-first, route second	Route-first, cluster second	Genetic Algorithms
	Figure A.1 : Classication des méthodes de résolution pour les VRP.
	Recourse) et les problèmes avec contraintes probabilistes (CCP en anglais pour Chance Constraint ce type de méthodes. Actuellement, le CVRP est résolu pour des instances allant jusqu'à 200 n÷uds,
	Programming). alors que ce chire diminue à 100 n÷uds pour le VRPTW. Les méthodes approchées essaient de
	Les modèles type SPR utilisent des actions appelées recours qui permettent de réagir face aux trouver un compromis entre la qualité de la solution et le temps d'exécution. Néanmoins, la plupart
	situations de violation de contraintes (échecs) suite à la révélation des paramètres incertains. Les de ces méthodes n'orent aucune garantie d'optimalité de la solution obtenue. Parmi les méthodes
	modèles CCP quant à eux visent à limiter la probabilité de violation des contraintes. Les deux modèles approchées, les métaheuristiques sont très utilisées car elles donnent des solutions souvent assez proches
	de l'optimum.

Ceci a donné lieu à un grand nombre de publications depuis l'introduction du VRP en 1954. Toutefois, la plupart des travaux continuent à considérer que les informations, et donc les paramètres des problèmes sont connus à l'avance. Cette supposition est rarement vraie dans la réalité, en eet il existe plusieurs facteurs qui peuvent remettre en cause cette supposition sur la certitude des paramètres.

Dans le contexte du transport urbain par exemple, les temps de trajets peuvent être aectés par les embouteillages. Les temps nécessaires pour servir les clients peuvent aussi dépendre de la complexité des services à eectuer, etc.

Les problèmes de type VRP avec incertitudes ont sucité un intérêt grandisant ces dernières années.

Gendreau et al.

[13]

ont proposé récemment une des revues les plus complètes sur le sujet, démontrant une activité de recherche croissante. Dans cet article, les auteurs se sont concentré sur les diérents modèles de programmation stochastique dédiés aux problèmes de tournées avec incertitudes. Pour analyser et classier les travaux, deux caractéristiques principales sont en général considérées : le type de paradigme de résolution considéré et les paramètres entâchés par les incertitudes (demandes, présence de clients, temps...etc.). Les paradigmes de résolution se divisent en deux catégories : l'approche d'optimisation à priori présentent cependant l'inconvénient d'êtres très lourd à résoudre en termes du temps nécessaire pour le calcul des coûts relatifs aux actions de recours, et des probabilités de satisfaction des contraintes.

Ainsi, il est crucial de développer des méthodes de résolution qui permettent de gérer les incertitudes d'une façon ecace et ceci aussi bien pour les modèles SPR que CCP. Cette thèse utilise la programmation stochastique et le paradigme à priori pour étudier deux VRP stochastiques et un problème de planication de la maintenance de parcs d'éoliennes. Dans A.1. Dans la seconde partie de ce chapitre, le problème précédent est étendu sur un horizon de planication long tout en considérant les stratégies de maintenance. Ce dernier problème est étudiée et abordée par une méthode basée sur la simulation. Les resultats obtenus montrent l'impact du choix des règles de priorité utilisées pour l'ordonnancement des tâches de maintenance sur les coûts et la quantité d'énergie produite. La thèse se cloture avec le chapitre 6 en orant des pistes pour des recherches futures sur l'ensemble des problèmes étudiés. A.1 Introduction aux VRP Dans le chapitre 2 une introduction générale aux problèmes de tournées des véhicules (VRP en anglais) est présenté ainsi qu'une étude plus détaillée des problèmes de type VRP avec incertitudes. Les VRPs ont été largement étudiés depuis leur introduction par Dantzig et Ramser [7]. Ceci peut être expliqué par deux raisons : l'importance du transport dans les activités humaines (distribution des produits et services), fois chaque client. Chaque tournée est une séquence ordonnée de n÷uds r = {r 0 = 0, r 1 = 0, . . . , r j , . . . , r k , r k+1 = 0} qui démarre et nit au dépôt. Ainsi, le coût d'une tournée particulière r est calculé par C r = k j=0 c rj ,rj+1 . Beaucoup de travaux sont consacrés au CVRP, cependant l'existence de plusieurs cas particuliers dans la réalité a donné naissance à de très nombreuses variantes. Le lecteur peut se référer au livre de Toth et Vigo [25] pour une revue de la littérature sur les variantes du CVRP. Une de ces variante est le VRP Pour résoudre le VRPTW (et en général les VRP) beaucoup de méthodes ont été proposées. En eet, résoudre les VRPs n'est pas une tâche facile car cette catégorie de problèmes fait partie des problèmes NP-Diciles.

La littérature dédiée aux VRPs continue de s'accroitre mais une grande partie des travaux suppose que les paramètres des problèmes sont connus à l'avance ou d'une façon déterministe

[20]

. Cependant, dans la réalité certains paramètres ne peuvent pas être connus avec certitude à cause des conditions météorologiques, les accidents, de la présence ou non de la clientèle, etc. Les travaux publiés dans la littérature montrent qu'ignorer les incertitudes conduit à des solutions infaisables et couteuses. Le A.1.

Table A .

 A 1 : Taxonomie des VRPs basée sur l'article de Pillac et al.[20] Qualité de l'information Dans le cas statique et déterministe, les paramètres sont considérés comme connus dès la planication de la solution. Dans le cas déterministe et dynamique les paramètres (ou une partie d'entre eux) sont complètement inconnus et sont seulement révélés à des moments spéciques. Le cas avec incertitudes partage une caractéristique avec le cas dynamique, étant donné que les vraies valeurs des paramètres ne deviennent connues qu'à des moments précis. Dans le cas statique et incertain, on dispose d' informations exploitables sur l'incertitude des paramètres (leurs lois de probabilité, les intervalles dont lesquels ils prennent leurs valeurs, etc.). Ces informations sont donc utilisées pour résoudre le problème. Finalement dans le cas dynamique et incertain, les paramètres (ou une partie d'entre eux) sont inconnus, mais comme dans le cas statique et incertain, il existe des informations relatives Enn, la logique oue permet aussi de représenter les incertitudes en utilisant des variables oues. Certains travaux combinent deux approches en utilisant par example des lois de probabilité pour modéliser les paramètres du problème et des nombres ous pour modéliser leur espérance ou leur variance[15].Les méthodes du type à priori sont basées sur des solutions statiques et peuvent être divisées en deux catégories : les problèmes avec recours (SPR) et les problèmes avec contraintes probabilistes (CCP). Les premiers utilisent des actions appelées "`recours"' qui permettent de rétablir la faisabilité de la solution quand des échecs se produisent.Un exemple pour le VRPSD est de revenir au dépôt quand la demande d'un client dépasse la capacité disponible du véhicule. Après le passage au dépôt, le véhicule reprend la tournée depuis le client où l'échec s'est produit. En outre, le véhicule complète la demande restante avant de servir les clients suivants planiés dans la tournée. Les coûts associés à ce type d'actions sont rajoutés dans la fonction-objectif du problème. D'autre part, les problèmes avec des contraintes probabilistes cherchent à limiter la probabilité des échecs à un seuil. Les CCP sont recommandés quand la dénition du recours est trop dicile ou quand un niveau de service doit être garanti. Une façon commune de classier les VRP stochastiques (SVRP) consiste à considérer les paramètres entâchés par les incertitudes. Selon Gendreau et al. [13] trois catégories peuvent être considérées : les VRPs avec demandes stochastiques (VRPSD pour VRP with Stochastic Demands), VRPs avec incertitudes sur la présence des clients (VRPSC), et les VRPs avec temps stochastiques (VRPST). Le problème avec demandes stochastiques a été le plus étudié et les modèles avec recours ont été privilégiés par apport aux modèles avec contraintes probabilistes. Le recours classique considère que lorsque la capacité d'un véhicule est épuisée, le véhicule fait un retour au dépôt pour s'approvisionner puis revient chez le client où l'échec s'est produit. D'autres recours existent pour le VRP avec la même politique de réapprovisionnement mais des visites au dépôt pouvant être eectuées avant que la capacité du véhicule ne soit atteinte. De cette façon, des économies en temps et en distance parcourue peuvent être réalisées. Ils existent d'autres recours plus complexes toutefois les recours simples ont été favorisés. Pour résoudre le VRPSD des méthodes exactes ont été proposés pour les cas où des recours simples sont utilisés. Cependant c'est les mpethodes apporchées les plus utilisés. Ces méthodes sont souvent testés sur un ensemble standard d'instances, toutefois la taille de ces dernières demeurent petites.Les problèmes avec clientes stochastiques sont ceux dans lesquels la présence des clients est incertaine. Quand le véhicule arrive chez le client, ce dernier peut-être présent ou non avec une certaine probabilité. Les VRPCS sont les moins étudiés parmi les VRP stochastiques. En eet, les VRP avec temps (trajet ou service) stochastiques (VRPST) ont reçu plus d'attention que les VRPCS mais demeurent moins étudiés que les VRPSD. La plupart des travaux consacrés à cette catégorie de problèmes considèrent des fenêtres de temps sur le service. Ces fenêtres de temps peuvent être souples ou dures , les premières étant largement favorisées. Ceci peut être expliqué par l'eet des fenêtres de temps dures sur les temps d'arrivée chez les clients. Généralement, les temps de trajets sont représentés par des lois qui ont des propriétés de convolution, mais les fenêtres dures empêchent l'utilisation de ces propriétés. Des problèmes qui prennent en compte les deux types des fenêtres sont les plus étudiés.Ces travaux considèrent que la date au plus tôt de service (début de fenêtre de temps) doit absolument être respectée mais autorise le service après la date de fermeture des fenêtres. La complexité des VRPST explique le faible nombre de publications utilisant des méthodes exactes. L'utilisation de ces dernières est limitée aux cas dans lesquels les variables aléatoires sont discrètes, ou quand l'espace de scénarios réduit, où au cas de variables aléatoires additives (avec des convolutions possibles à calculer).De même que pour les autres VRP stochastiques, les méthodes approchées restent les plus privilégiées pour résoudre le VRP avec temps stochastiques.La complexité des VRP avec incertitudes a limité le développement des recherches sur ce sujet ainsi que la taille des instances résolues. Toutefois, comme pour le cas déterministe, il existe un réel besoin de méthodes puissantes pour résoudre les problèmes dans des conditions réelles, autant en termes de taille, de nature de paramètres aléatoires etde contraintes complexes. Donc, dans cette thèse on propose d'étudier les VRP de nature stochastique et de développer des approches adaptées capables de résoudres des instances de grande taille. Ce travail est l'un des rares à considérer des fenêtres de A.2. UNE MÉTHODE HYBRIDE POUR LES VRP AVEC DEMANDES STOCHASTIQUES temps dures sur le service tout ayant des temps de trajet et service stochastiques. A.2 Une méthode hybride pour les VRP avec demandes stochastiques Le chapitre 3 présente une méthode de résolution approchée pour le VRP avec demandes stochastiques (VRPSD). Un nouvel ensemble d'instances est proposé pour permettre des comparaisons futures. La méthode de résolution développée pour ce problème est un algorithme mémétique (MA) hybridé avec une méthode gloutonne du type Greedy Randomized Adaptive Search Procedure (GRASP). L'algorithme proposé montre une meilleure performance que les méthodes trouvées dans la littérature. En eet, pour les instances "`classiques"' de Christiansen et Lysgaard [4] notre méthode, trouve toutes les meilleures solutions connues dans un temps de calcul réduit. Le VRPSD est une généralisation du CVRP dans lequel les demandes des clients sont représentées par des variables aléatoires. Le problème étudié dans ce chapitre considère une version du VRPSD avec les caractéristiques suivantes. La demande q i de chaque client i est modélisée par une variable aléatoire qui suit une loi de probabilité ψ avec espérance E [q i ] > 0 et variance V ar [q i ] > 0. Il est supposé que la loi ψ est connue et que les demandes sont indépendantes entre elles. De plus, il est considéré que la loi de probabilité de la somme des variables ψ est aussi une variable aléatoire de même nature ψ.

			Données déterministes	Données incertaines
	Évolution de l'information	Données connues à l'avance Données dynamiques	Statiques et déterministes Statiques et incertaines Dynamiques et déterministes Dynamiques et incertaines
	Parmi les trois approches, la programmation stochastique est la plus répandue pour résoudre
	les VRPs avec incertitudes. Dans cette approche, les paramètres sont modélisés par des variables
	aléatoires de lois connues. De ce fait, les contraintes des problèmes peuvent ne plus être respectées
	par les solutions. Par exemple, dans le VRP avec demandes stochastiques (VRPSD), la réalisation
	de la demande d'un client peut dépasser la capacité restante du véhicule. Si une contrainte n'est pas
	respectée par la réalisation des parametres du problème un "`échec"' est dire de se produire. Pour les
	considérer, deux types des modèles sont utilisés [12] : la réoptimisation et les approches à priori avec
	recours.

à ces derniers. Une autre diérence fondamentale entre les problèmes statiques et dynamiques est la façon dont laquelle les solutions sont calculées. Dans le cas statique, une solution reste non modiable quelques soient les vraies valeurs des paramètres. Dans le cas dynamique en revanche, la solution peut être constamment modiée pour s'adapter aux informations qui arrivent au fur et à mesure.

Trois approches ont été principalement utilisées pour modéliser les incertitudes des paramètres des VRPs, à savoir : la programmation stochastique, l'optimisation par intervalles, et la logique oue. L'optimisation par intervalles modélise les paramètres incertains par des intervalles de valeurs possibles. L'objectif poursuivi par l'optimisation par intervalles est de trouver des solutions qui sont faisables pour toutes les réalisations possibles des paramètres

[2]

. Des VRPs avec incertitudes sur les demandes et les temps apparaissent dans la littérature. Adulyasak et Jaillet

[1]

traitent la version avec temps de trajets incertains. Les auteurs proposent une comparaison entre l'optimisation stochastique et l'optimisation par intervalle. Les auteurs montrent que les solutions robustes surpassent largement les solutions issues de l'optimisation stochastique dans le cas où les paramètres ne sont pas modélisés par la bonne loi de probabilité. Les modèles qui utilisent la réoptimisation ne reposent pas sur une solution xe. En eet la solution est construite et modiée au fur et à mesure que les informations apparaissent. Dès révélation de ces dernières, les tournées peuvent-être réoptimisées à nouveau. Toutefois, ce type d'approche rend la coordination des véhicules assez ardue et la vitesse à laquelle les solutions nécessitent d'être fournies reste un problème.

E [C r ] =

Table A .

 A 2 : Synthèse des résultats sur le nouvel ensemble d'instances BKS BDS Instance |V | Min veh Q FC Veh Total Det Rec Total Rec VSS A-n61-k9 61 9 100 0.98 10 1144.23 1084 60.23 1215.37 181.37 5.85% A-n62-k8 62 8

Table A .

 A 3 : Comparaison des diérentes MPMA Instances avec 100 clients

	Metric	MA1 MA1 + MA2 MA1 + MA3 MA1 + MA2 + MA3
	# Véhicules Meilleures solutions # Véhicules moyenne Coût déterministe moyen Coût du recours moyen	579 616.7 1156.56 3.86	575 612.7 1158.88 3.90	576 608.4 1138.24 3.88	554 577.3 1123.87 31.16
	Temps moyen pendant MA (s) 100 Temps moyen simulation (s) 35.28 Temps total moyen (s) 135.28 165.06 100 65.06	100 61.85 161.85	100 54.70 154.80

Table A .

 A 4 : Comparaison du MPMA avec la métaheuristique ILS de Miranda et Conceição [18] 16.10 1754.02 81.42% 97.28% 0.18% 1.05% 19.20 MPMA 15.99 1428.49 80.86% 97.20% 0.44% 3.89% 32.04 classiques de Solomon [24] dédiées au VRPTW. Les deux articles considèrent des lois continues pour les temps de trajet et de service. Toutefois chacune des deux études a ses propres particularités. Miranda et Conceição [18] modélisent le SVRPTW comme un CCP avec un taux de service α égal à 80%. Ainsi, les temps de service et de trajets suivent des lois normales avec des espérance égales à leurs valeurs déterministes et la variance est déduite à partir des coecients de variations générées par une loi uniforme U [0.1; 0.6]. tableau A.4 montre les métriques suivantes : nombre de véhicules, la distance moyenne parcourue, les niveaux de services (NS) et les écarts entre le niveau de service (ENS) calculé par simulation et celui prédit par l'approximation de AT i . On peut remarquer que MPMA utilise moins de véhicules en moyenne qu'ILS avec une distance parcourue beaucoup plus réduite. En plus MPMA est le seul capable de garantir le niveau de service requis pour toutes les instances (α = 80%). D'une autre part, le travail de Nguyen et al. [19] utilise des temps de trajets avec une distribution de probabilité Gamma. L'espérance des variables est égale à la valeur déterministe et l'écart type est de 50% l'espérance. Le tableau A.5 résume les résultats pour le nombre moyen de véhicules, le coût de la solution (distance), et le retard moyen (RA) chez chaque client. Dans un objectif de comparaison, le MPMA est modié pour n'utiliser que le coût du recours sans tenir compte des contraintes probabilistes. Le retard minimal moyen est atteint par MPMA de même que le coût moyen minimale. Ceci est réalisé au détriment d'une petite augmentation sur le nombre de véhicules en comparaison avec les diérentes versions de la Recherche Tabu (TS) de Nguyen et al. [19]. Ces auteurs considèrent le retard moyen comme objectif principal. On peut conclure que MPMA donne les meilleurs résultats par rapport à ces deux travaux. .3. UN ALGORITHME PARALLÈLE POUR LES VRP AVEC TEMPS DE TRAJET ET TEMPS DE SERVICE STOCHASTIQUES Table A.5 : Comparaison de MPMA avec la TS de Nguyen et al. [19] Il reste dans la continuité de ce travail de tester d'autres types de recours, et l'extensions du problème pour considérer les temps de trajets corrélés. Aussi, l'approximation de AT i peut être exploitée dans des méthodes exactes de type Branch-and-Price. Finalement, on prévoit d'étudier des problèmes multi objectif dont le niveau de service est un objectif et pas une contrainte. Des résultats présentés dans ce chapitre ont été présentés à la conférence MIM2016 : Gutierrez, A., Dieulle, L., Labadie, N., Velasco, N. (2016) A multi population memetic algorithm for the vehicle routing problem with time windows and stochastic travel and service times In 8th IFAC Conference on Manufacturing Modelling, Management and Control MIM 2016 Troyes, France, 2830 June 2016 Une extension de ce travail pour résoudre le problème avec une méthodes exacte a été présentée à la conférence CLAIO2016 : Gutierrez, A., Dieulle, L., Labadie, N., Velasco, N. (2016) An approximate column generation for the vehicle routing problem with hard time windows and stochastic travel and service times. In XVIII CLAIO, the Latin-Iberoamerican Conference on Operations Research Santiago de Chile, Chile, 26 October 2016. Un article regroupant les résultats de ce chapitre a été accepté pour publication dans le journal Computers & Industrial Engineering. A.4. GESTION DES RESSOURCES POUR LA MAINTENANCE D'UN PARC D'ÉOLIENNES cherche à satisfaire son propre objectif, à savoir, minimiser les coûts d'une part et maximiser l'énergie produite de l'autre part. Le problème considère un ensemble d'éoliennes pour lesquelles des tâches de maintenance doivent être exécutées. Un nombre limité de techniciens ayant chacun leurs propres compétences est dédié à la réalisation de ces tâches. Pendant l'exécution des tâches de maintenance, les éoliennes sont arrêtées et ne produisent donc pas d'énergie. De plus, toute tâche commencée doit être achevée. Ainsi, à la n de chaque journée, si la tâche n'est pas encore nie, les techniciens la poursuivront le lendemain. Chaque tâche est associée à une fenêtre de temps dont l'ouverture peut correspondre par exemple à la date de disponibilité des matériaux nécessaires à l'éxécution de la tâche de maintenace. La fermeture de la fênetre est considérée comme une date à laquelle la tâche devrait nir. Si la fenêtre de temps n'est pas respectée, une pénalité est appliquée. Chaque tâche nécessite un certain nombre de compétences et de techniciens par compétence. Le temps est discrétisé, et pendant chaque créneau de temps les éoliennes qui ne sont pas en maintenance produisent de l'énergie. Pour chaque créneau, la vitesse du vent est supposée connue ce qui permet de savoir si les conditions permettent ou non d'eectuer la tâche prévue. Dans le cas où la vitesse est trop élevée, les techniciens doivent arrêter leur travail. En plus le salaire des techniciens dépend de leurs compétences. Trois coûts sont considérés dans ce chapitre, le coût lié à l'aectation de techniciens surqualiés, le coût lié aux heures supplémentaires, et la pénalité en cas de retard d'exécution de la tâche de maintenance.

	Instance Mé-thode R105 ILS MPMA 17.30 1470.89 80.55% 97.30% 0.39% # Veh. Moyen CD Moyen NS Min Moyen NS Moyen ENS Moyen 17.67 1615.75 81.66% 99.31% 0.18% R109 ILS 15.00 1488.81 79.53% 95.45% 0.17% MPMA 15.00 1270.27 80.96% 98.00% 0.28% C101 ILS 17.00 2284.77 83.82% 95.52% 0.22% MPMA 17.00 1372.27 80.42% 95.42% 0.64% C106 ILS 14.67 1722.32 81.33% 98.65% 0.18% MPMA 14.00 1141.17 80.57% 96.01% 0.72% RC101 ILS 19.67 2012.76 78.32% 100.22% 0.16% MPMA 19.60 1841.32 80.98% 97.53% 0.46% RC106 ILS 14.67 1584.63 82.26% 95.57% 0.16% MPMA 15.30 1534.95 81.07% 97.87% 0.40% RC107 ILS 14.00 1569.10 83.03% 96.23% 0.16% MPMA 13.70 1368.53 81.46% 98.44% 0.23% Avg. C1 TS (AD1) 10.33 944.55 2.82 TS (AD2) 10.33 929.85 3.43 TS (AD3) 10.22 918.61 4.19 MPMA 9.89 908.73 5.55 C2 TS (AD1) 3.25 654.41 0.00 TS (AD2) 3.13 648.60 0.04 TS (AD3) 3.00 638.40 0.02 MPMA 3.00 606.30 0.08 R1 TS (AD1) 15.67 1521.12 30.77 TS (AD2) 15.83 1495.02 24.02 TS (AD3) 15.58 1447.70 26.02 MPMA 15.51 1409.47 19.25 R2 TS (AD1) 3.45 1134.96 5.33 TS (AD2) 3.55 1050.33 1.48 TS (AD3) 3.45 1086.00 1.77 MPMA 3.99 991.07 1.96 RC1 TS (AD1) 15.25 1736.57 49.79 TS (AD2) 15.38 1713.98 34.81 TS (AD3) 15.13 1700.99 41.07 MPMA 15.86 1630.91 22.55 RC2 TS (AD1) 4.00 1396.98 4.95 TS (AD2) 3.88 1296.20 2.64 TS (AD3) 4.00 1234.10 3.17 MPMA 4.56 1182.94 3.86 Avg. TS (AD1) 8.91 1241.83 15.91 TS (AD2) 8.95 1198.80 11.35 TS (AD3) 8.82 1181.67 12.92 MPMA 9.04 1131.34 9.19 Le InstanceMéthode# Veh. MoyenDC MoyenRA MoyenTemps (s) Max. ENS Temps (s) 1.81% 18.40 2.81% 31.82 0.89% 6.36 3.44% 32.21 0.87% 29.49 3.56% 31.77 0.70% 5.56 4.78% 32.25 1.42% 15.18 4.98% 31.77 0.82% 24.43 4.64% 32.08 0.86% 34.99 3.01% 32.40 9.20 7.94 8.62 51.41 14.34 9.97 8.52 51.60 16.10 14.96 18.95 51.76 51.50 42.38 43.90 53.56 11.24 10.78 9.46 51.71 33.49 25.82 23.41 52.67 23.48 19.46 19.98 52.16 Gamma 2 .

ILS

Le modèle proposé dans le chapitre 4 ainsi que la méthode de solution MPMA ont montré leur exibilité. Les comparaisons avec d'autres méthodes de la littérature donnent l'avantage à MPMA, même si la méthode utilise des temps d'exécution légèrement supérieurs. Il est également important de noter que l'approximation de AT i avec des loi Log-normales permet d'avoir des résultats compétitifs pour les problèmes même si les temps de trajet (et de service) sont distribués par des loi normales où A

Nevertheless, the parameters can be also modeled by sets (Robust Optimization) or by Fuzzy Variables. More information is given in chapter 2

The authors use the date 1954 as the rst VRP record in the literature considering the work of Dantzig et al.[40] on the Traveling Salesman Problem (TSP), a particular case of the VRP

Optimal deterministic solutions retrieved from http://vrp.atd-lab.inf.puc-rio.br or https://www.coinor.org/SYMPHONY/branchandcut/VRP/data/index.htm for all instances except tai385 for which the deterministic BKS from the rst link is used.

 2 We use the same approach of Gauvin et al.[15], by evaluating the routes in the optimal deterministic solution to the right and the reverse, retaining the one with the minimum cost

Instance CMT12 appears in Gauvin et al.[15] as M-n101-k10

E-n23-k3 instance has a FC which is nearly

standard deviations from the mean FC5 Tai385 has a number of vehicles in the BKS which is over 5 standard deviations from the mean

Except for instances CMT4, Tai 150a, 150c, 150d, and 385 for which the BKS reported was found during this long run, a slightly bigger BKS was used during this tests

I thank professor Christian Prins for pointing this out during the dissertation questions and comments.

The number of trials was increased for this test to

30,000 to have a better accuracy on the condence intervals, and thus conclude if two solutions are dierent on average

Refers to the percentage of time that a wind turbine is available to produce power

Note that min Z 2 is equivalent to max -Z 2

Each of the sub-components is not further divided into sub-components

Based on Carlos et al. [?] and meetings with experts on the eld.

Only four technicians were considered.

Toutes les instances sont disponibles dans la page : http ://vrp.atd-lab.inf.puc-rio.br

Avec des coecients de variation de plus de

35%

Acknowledgments

While considering that Solomon [41] instances are designed for the VRPTW, we performed some modications to adapt them to deal with the stochastic nature of the SVRPTW. The service levels α, β, γ are set to 95%, 95%, and 90% respectively. To guarantee a feasible solution, the elementary routes (routes visiting only one customer) dened as r = {r 0 = 0, r 1 = i, r 2 = n + 1} ∀i ∈ V c must at least respect constraints (4.1) and (4.2). To guarantee the feasibility of these routes, each one was tested through one million trials by simulation. If a given route r is unfeasible because of constraint (4.1), the customer time window closure l i is set to l i = l i + 5. If infeasibility arises because of constraint (4.2) then, the depot time window closure is set to l n+1 = l n+1 + 5. This process is repeated until feasibility is reached for every route visiting a single customer.

General discussion

The MPMA presented in section 4.4 was tested on the modied instances of Solomong [41] introduced in section 4.5.1. We used the instances including the rst 50 customers and the complete instances.

MPMA was implemented on a Dell Latitude E6420 personal computer with Intel Core i7-2760QM @2.4 GHz, running Windows 7 Professional 64 bits. The algorithms were coded in Java and compiled with JavaSE-1.8_45, with maximum allocated memory of 1 Gb. Random variables and computation probabilities were generated by the library of Stochastic Simulation in Java ( [30]). Parameter values were selected after several preliminary tests. The vehicle cost M is xed to 1000, this value allows the MPMA to minimize the number of vehicles as the rst objective followed by the total distance plus recourse cost. MA populations size P opSize is set to 25 for the three MAs. When simulation is used during the MPMA it performs 1000 trials to evaluate the feasibility, estimate the cost and required probabilities for each route. The MPMA performs 10 runs for the whole set of instances (50 and 100 customers). is set to n/10 seconds, while τ value is 0.05 seconds. The p ls and p m are set to 0.2 for all MAs although MA1 starts without performing LS. The MAs stopping condition is set to 20 seconds for the 50 customer instances while this value is increased to 100 seconds for the instances with 100 customers. Set partitioning problem is solved by means of Gurobi 6.0 ( [20]). Tables 4.4 to 4.6 presents the best solutions found out of the 10 runs for C, R, and RC instances respectively. Meanwhile, average results can be found at tables 4.11 and 4.12. For each instance the expected route cost (ERC), the number of vehicles (# Veh.), the distance (DC), and the recourse cost (RC) are reported. Type 1 instances show higher recourse costs than type 2 instances. This is quietly pronounced in family R (100 customers) since the average recourse cost is increased from R2 to R1 by 605%; similarly families C and RC display the same behavior. This is completely coherent since type 1 instances have tighter time windows which imply that failures are more likely to occur, and consequently recourse costs are higher. Furthermore, the family C instance shows an increment of nearly 106% in the cumulative number of vehicles (CNV) required for 100 customer instances when compared to 50 customer instances. The rise of the total costs is near to 105%, while distance and recourse increases by 96.9% and 83.1% respectively. Family R on the other hand, shows the minimum increments (among C, R, and RC families), with an augmentation of nearly 72.8% in the CNV, 70.3% in terms of cost and 50.8% and 71.1% for distance and recourse metrics. Moreover, RC family raises the CNV by 79.5% and the total cost by 78.1%. The distance is increased by nearly 69.2%, and recourse presents the lowest increment of only 26.1%. Also, it should be noted that the best CNV for 50 customer instances is 302 while the same metric for 100 customer instances is 554. This dierence of 252 vehicles is translated into an increase of 83.4%. A similar behavior was observed for the deterministic best solutions found by MA2, which incremented the CNV by a factor of three quarters when passing from 50 customer to 100 customer instances. Thus, the increments in the CNV seem to be driven mostly by the inherent combinatorial nature of the problems whether the stochastic or deterministic versions are considered.

Operational maintenance level

Operational decision level makes reference to day-to-day decisions. Within this level, maintenance planning refers to the problems in which maintenance tasks are scheduled to be performed on a xed short length horizon. The detailed scheduling may include constrained resources, such as technicians, vessels, vehicles, etc.

In this context, it is usually assumed that the maintenance actions are known prior to their execution [8,40,31]. Additional information on the tasks such as their release dates and due dates [72] are also known. Since the horizon is limited to very few days, it is a common practice to assume that stochastic parameters such as wind speed or the waves height can be perfectly predicted. In cases where parameters can change, a roll over approach is used, re-optimizing each current day scheduling with the new predictions of stochastic parameters and the realizations over previous periods. Table 5.1 shows a summary of the principal works dealing with the maintenance scheduling problem. It focuses on the types of objectives used in the dierent literature works. The column costs makes reference to objectives associated with costs, while the energy one refers to objectives associated with the energy production. The following conventions are used for the table: Onshore (ONS), oshore (OFF), penalties (PEN), transportation (TRAN), technicians wage (TWAG), technicians extra wage (EWAG), production loss (PL), energy produced (EP), and revenue (REV).

From table 5.1 one can see that few onshore works have relied on energy objectives. In this context, Kovàcs et al. [40] assume that crews formed by two technicians are used to perform maintenance tasks.

A Mixed Integer Linear Problem (MILP) on a rolling horizon is designed to minimize the total loss of production due to downtime of the turbines, or when they are still working but in a degraded state.

Two types of degradations are considered, general and peak, the rst one diminishes the power output by a percentage in any operation condition, while the second one reduces the production during high speed winds but not on lower speed ones. Albeit no detailed computational results are given, the authors claim to solve instances with up to 50 maintenance tasks, 4 teams, and 7 wind farms. More recently Froger et al. [31] have also considered the onshore wind farm maintenance scheduling. The authors dealt with multiple technician skills, dierent types of execution modes for the tasks, as well as dierent farm locations. Moreover, the objective in the proposed model is to maximize electricity production over a short-time horizon. It must be noticed that maximizing the revenue and the energy produced are only equivalent when the energy prices are constant. Since this is the case in [31,30] both objectives are selected. To solve the problem, two formulations based on Integer Linear Programming (ILP) are proposed. A constraint programming large neighborhood approach is devised to solve the problem. The same authors worked a second paper dealing with this previous problem. They tackled it in this last study [30] with a branch-and-check approach. This method can consistently produce optimal or near optimal solution for instances with up to 80 tasks, several modes, skills, and farms locations.

CHAPTER 5. WIND FARMS MAINTENANCE

In the oshore context minimizing the transportation costs is omnipresent. This can be explained by the high costs associated to vessels and helicopters. Besnard et al. [8] propose a stochastic optimization model for wind farms maintenance planning. Stochasticity is incorporated through scenarios in which the wind, waves and production take dierent values. The model is based on a rolling horizon which considers opportunistic and corrective maintenance actions. The horizon planning is constituted of seven days and is discretized in steps of one day. Scenarios are used to characterize the expected hourly production power, the wind speed and wave height; and the objective is to minimize production losses, transportation costs, and extra hours penalties. A similar model is proposed by Kennedy et al. [39]. The authors consider only one vessel and one maintenance team.

Using a Genetic Algorithm, the authors optimize the order and time at which tasks must be carried out. The objective is to minimize the production losses, transportation and crew costs. Also, it is shown that a signicant saving (from 13% to 21%) can be achieved when the maintenance schedule is optimized instead of repairing the turbines as fast as weather conditions allow access to the turbines.

Problems involving more than one vessel are the most common type in oshore context. Two models based on arc-ow and path-ow formulations are presented in Stålhane et al. [79] for the routing and scheduling of vessels that perform maintenance tasks at oshore wind farms. The presented models consider as objective the minimization of transportation, downtime, and penalties costs.

Instances considering a workday and at most eight tasks and ve vessels are solved to optimality.

A similar problem is tackled in Dai et al. [19]. The authors use a four index MILP where vessels daily availability depends on the type of the vessel and the weather conditions. Numerical results are provided for instances with eight turbines and an horizon of three days, aiming to minimize the costs and production loss. A combined vessel routing and maintenance scheduling problem is tackled by Raknes et al. [64]. The model consider multiple wind farms which are managed by the same O&M enterprise. Dierent types of vessels are included, some of which can stay oshore for longer periods, and other need to return to the depot at the end of each shift. The authors propose a large Mixed Integer model to minimize the sum of production loss, transportation costs, and penalty costs related to tasks that are not carried out. Furthermore, a dynamic version of the problem is considered dealing with new task arrivals and updated weather forecasts. The results show the importance of evaluating the strategies by simulation rather than considering static models in dynamic contexts.

The maintenance routing and scheduling at oshore wind farms is tackled by Irawan et al. [35] thanks to a model that considers multiple vessels, periods (days), bases, and wind farms. Weather conditions are considered by dening maximum working hours for a vessel. A Dantzig-Wolfe decomposition is used, where for each vessel and each period, all the feasible scheduled routes (customers visits) are generated. The time horizon varies from three to seven days, with three types of technicians.

The objective function contains the cost of technicians, transports, and penalty costs when turbines are visited after a given deadline. Results are presented for literature instances with eight turbines and three periods. The authors also show that O&M with multiple bases and attending multiple wind farms can produce savings of around 12% when compared to solutions with single O&M base and wind farm pairs. Yet, it is interesting that Irawan et al. [35] nd that costs are mostly composed by crew costs. Certainly, that behavior could be expected to appear in the onshore case where technicians are the most costly resource (apart from spare parts).

From the explored literature almost all the works have used exact methods to solve the problems. This presents the drawback of ignoring the stochastic parameters to have tractable models. Consequently, the works assume that weather conditions, or the time to perform a task are known in advance. Moreover, no work has addressed the multi-objective case in which both minimization of the costs, and the maximization of energy related objectives are tackled at the same time. Nowadays, most of the proposed strategies are hybrids in which corrective maintenance is combined with some proactive maintenance, either preventive and predictive approaches [78]. For example, Lu et al. [47] studied a condition-based maintenance strategy coupled with opportunistic maintenance for oshore wind farms. In this work, an articial neural network (ANN) approach is used to predict the component life percentage based on monitoring information. Also, by using the information from the ANN, a conditional failure probability is determined. The authors use a two failure probability thresholds to decide when to perform a preventive or an opportunistic maintenance task. The authors conclude that the proposed strategy reduces the cost by nearly 30% when compared to a time-based (preventive) strategy in onshore and oshore contexts. A very similar method can be found in [81], the authors compare a condition-based maintenance against a periodic maintenance policy. The conditionbased maintenance use two threshold values to decide when the turbines needs to be maintained. The authors conclude that a condition-based strategy can reduce by more than 40% the maintenance cost for a ve turbines wind farm example. Shaee et al. [75] propose an opportunistic condition-based maintenance for wind turbine blades.

The authors optimize two decisions: the length of a crack at which a major repair must be performed, and the operational age of the blades that triggers a preventive maintenance. In major repairs and preventive actions, the other blades also receive a maintenance. The proposed method achieves to reduce by more than 20% the O&M costs when compared to a reactive strategy. A higher saving (30%) on O&M costs is generated by the proposed strategy against an individual strategy (per blade) without using opportunistic maintenance. A similar work is presented in Shaee et al. [74] where the intermediate speed shaft and the high speed shaft are considered. The same parameters (crack threshold and operational age) are optimized, although, each component has its own threshold. In Hameed and Vatn [33] the opportunistic strategy is extended to consider components of a turbine and among turbines. Still, opportunistic maintenance is preferred to be performed on the same turbine to avoid traveling costs from turbine to another. By using this strategy setup costs are shared among dierent maintenance activities, thus reducing the total costs. The approach to group the tasks works as follows. First, an optimal interval for individual components for individual activities is found.

Then, using a heuristics the groups are created. An example of three wind turbines each with eight components for a horizon length of four years is analyzed. The results show that approximately the group of activities have a periodicity of around three months. However, further experimentation is needed to evaluate the pertinence of this type of grouping. Especially, since the appearance of failures or bad weather times can completely change the planning.

Other studies on strategies selection involving opportunistic maintenance can be found in [22,69,92,90,28]. In [90] two level thresholds are used to decide on the basis of the reliability functions of component if minimal maintenance is performed or if a replacement is preferable. Moreover opportunistic maintenance is also considered as soon as a minimal maintenance or a replacement takes place.

In Erguido et al. [28] a dynamic opportunistic policy is evaluated. Dynamism is added through reliability thresholds to perform activities, changing in function of the weather conditions. The authors conclude that using dynamic thresholds surpass the use of static thresholds.

Although the previous works address the comparison of dierent strategies, they usually rely on many simplications. For example, considerations of the weather, the number of resources (technicians, vessels, etc), the time to perform the maintenance tasks, the necessary spare parts, among others are not fully addressed issues, or not considered at all. Therefore, other works have focused on incorporating such considerations so the complexity of wind farms operation and maintenance is better modeled. In this context, simulation has been the preferred approach to enable more tractable comparisons.

MAINTENANCE STRATEGIES: RELATION WITH OPERATIONAL PLANNING

wind speed, v 1 the rated wind speed, and v 2 the cut-out speed.

P E

Thus, when the weather module provides a speed to each turbine, this one can estimate by means of equation (5.23) the amount of energy produced.

Weather model -Wind speeds

Performing simulation of wind farms require to address a fundamental aspect, that is, a synthetic creation of wind data speeds. Several ways of dealing with this are reported in the literature, for example Carlos et al. [14] ws t -µ = Λ (ws t-1 -µ) + t (5.24)

Although the approach presented in [11] seems to oer a more robust way to generate correlated winds for many locations, the auto-regressive model of McMillan and Ault [57] is kept for its simplicity.

Moreover, it allows to consider the season patterns depending on the month, by varying the value of Λ, and the value of µ.

Schedule Modeler

Whichever the type of strategy used, at the operational decision level the resources utilization needs to be optimized. To the best of our knowledge, no other study has before considered to optimize the maintenance scheduling plan during a life cycle simulation of a wind farm. This is mainly due to the computationally demanding time to carry on such a task. Certainly, to embed the whole model and solution approach presented in section 5.4 in a long-term simulation model is unrealizable. This is why an approximate solution approach based on simple heuristics rules to assign the maintenance tasks is used.

In this part, failure-based maintenance, and periodic maintenance are compared. Strategy (CM) performs only maintenance tasks when a failure occurs. Meanwhile, (PM) policy also use corrective tasks when failures occur, but also, it performs periodic maintenance to avoid possible future failures.

To optimize the schedules a rolling basis to make operational decisions is devised. That is, each day the assignment of resources to maintenance tasks is redened. Nevertheless, tasks which have already begun are not rescheduled. Each day, at 7 a.m. it is assumed that the wind speeds for the 5.7. CONCLUSIONS the studied models will be extended to the oshore case to bring more understanding on the how to select the best strategy to investigate the maintenance tasks in dierent environments.

Contributions

The results of the multi-objective problem in this chapter were presented at IEOM conference: end if 20: end while cette nouvelle solution peut être modiée par des procédures de mutation et de recherche locale (lignes 6 et 8). Comme MA+GRASP est basé sur un algorithme génétique, des méthodes permettant de coder et de décoder une solution sont utilisés. Le décodage utilise la méthode Split de Prins [21], alors que le codage concatène toutes les tournées de la solution sans considérer le dépôt. Dans la ligne 17, et si le nombre d'itérations sans amélioration de la meilleure solution trouvée dépasse une valeur φ, La loi de probabilité de AT i est nécessaire pour la vérication des contraintes probabilistes et le calcul de la fonction-objectif. Dans ce travail AT i est aprochée par la variable AT ri de distribution Log-normale à l'issue d'une étude comparative de plusieurs lois. En utilisant les instances de Solomon [24], une série de tournées ont été construites avec l'heuristique de Clarke et Wright. Les temps de trajet ont été simulés en utilisant la loi Log-normale et les temps de service avec la loi Gamma.

Les espérances de ces variables sont égales aux valeurs déterministes des instances originalles, et les variances ont été calculées en utilisant des coecients de variation de 0.