A robust approach for solving a vehicle routing problem with time windows with uncertain service and travel times

The main purpose of this paper is to study the vehicle routing problem with hard time windows where the main challenges is to include both sources of uncertainties, namely the travel and the service time that can arise due to multiple causes. We propose a new approach for the robust problem based on the implementation of an adaptive large neighborhood search algorithm and the use of efficient mechanisms to derive the best robust solution that responds to all uncertainties with reduced running times. The computational experiments are performed and improve the objective function of a set of instances with different levels of the uncertainty polytope to obtain the best robust solutions that protect from the violation of time windows for different scenarios.


Introduction
Since the pioneer paper of Dantzig and Ramser (1959) on the truck dispatching problem appeared at the end of the fifties of the last century, work in the field of the vehicle routing problem (VRP) has increased exponentially.Using a method based on a linear programming formulation, the authors of this work produced by hand calculations a near optimal-solution with four routes of a fleet of gasoline delivery trucks between a bulk terminal and twelve service stations supplied by a terminal.Nowadays, vehicle routing problem is considered as one of the most outstanding research achievements in the story of operations research and particularly in practice.There are important advances and new challenges that have been raised during the last few years such as radio frequency identification, and parallel computing (e.g.Pillac et al., 2013;Montoya-Torres, 2015).The class of VRP problems involves minimizing a travel distance of vehicles, starting and ending from a depot, to serve some customers.Typically the solution has to obey several other constraints, such as the consideration of travel, service, and waiting times together with time-window restrictions.This variant is called in the literature vehicle routing problem with time windows (Bräysy & Gendreau, 2005;Kallehauge et al., 2005;Rincon-Garcia et al., 2015).
For instance, three types of solution approaches can be used to solve these types of problems.First, the exact methods assert that the optimal is found if the method is given sufficiently in time and space.We cannot expect to construct exact algorithms which solve NP-hard problems.Second, the heuristics are solution methods that can quickly achieve a feasible solution in a reasonable quality.A special class called metaheuristics provides a high solution quality (Labadie et al., 2016).The third class of solutions is also a special class of heuristic which provides a near-optimal solution and an error guarantee.An interesting topic on solving VRP consists in considering parameters affected by uncertainty, making the problem more realistic.
Different approaches have been proposed to deal with uncertainties in a VRP problem either in demand, travel time and/or service time.Among them, the stochastic approaches of vehicle routing problem SVRP have been treated in series of papers (Dror et al., 1993;Dror & Trudeau, 1986;Gendreau et al., 1996).The aim of the SVRP methodology is to find a near-best solution of the objective function responding to all possible data uncertainty.An alternative approach to handle the uncertain parameters is the robust optimization in which one can optimize against the worst scenario that can be generated from the source of uncertainty by using bi-objective function (Yousefi et al., 2017) and is immunized against this uncertainty (Sungur et al., 2008).In this context, the literature coats a large number of applications such as scheduling (Goren & Sabuncuoglu, 2008;Hazir et al., 2010), facility location (Minoux, 2010;Baron et al., 2011 ;Alumur et al., 2012 ;Gülpinar et al., 2013), inventory (Bienstock & Özbay, 2008;Ben-Tal et al.,2009), finance (Fabozzi et al., 2007;Gülpinar & Rustem, 2007;Pinar, 2007).In particular, the authors proposed a mathematical model for the robust optimization with uncertain demands (Moghaddam et al., 2012) and heterogeneous fleet (Noorizadegan et al., 2012) and routing with capacity (Sungur et al., 2008;Gounaris et al., 2013).For instance, this is equivalent to the deterministic case studied by Miller-Tucker-Zemlin formulation of the used VRP.We refer the reader to an excellent survey and tutorial of the robust vehicle routing proposed by Ordóñez (2010).We note that uncertainty in travel cost could be handled using the robust combinatorial optimization approach.Wu et al. (2017) proposed a linear model evaluated on a set of random instances for the vehicle routing problem with uncertain travel time to improve the robustness of the solution which enhance its quality compared with the worst case on a majority of scenarios.In the same spirit, Toklu et al. (2013) treated the VRP problem with capacity and uncertain travel costs based on a variant of the ant colony algorithm to generate sets of solutions of uncertainty levels and to analyze their effects on the problem.
The stochastic approach is also applicable for the vehicle routing problem with time window constraints (VRPTW).Errico et al. (2016) formulated the VRPTW with stochastic service times as a set partitioning problem and solve it by exact branch-cut-and-price algorithms.More precisely, they elaborated efficient algorithms by choosing label components, developing lower and upper bounds on partial route reduced the cost to be used in the column generation step.Unlike this approach, robust optimization seeks to get good solutions for the VRPTW problem by only considering nominal values and deviations possible uncertain data.Many works tackled the vehicle routing problem with time windows and uncertain travel times (Sungur et al., 2008).Agra et al. (2012) presented a general approach to the robust VRPTW problem with uncertain travel times.Travel times belong to a demand uncertainty polytope, which makes the problem more complex to solve than its deterministic equivalent.The advantage of the addition in complexity is that the model from Agra et al. (2012) is more usefule than the one from Sungur et al. (2008) and leads to less conservative robust solutions.Toklu et al. (2013) adapted their approach to solve the problem of VRPTW with uncertain travel times, whose objective was to minimize window time violation penalties by providing the decision makers a group of solutions found over several degrees of uncertainties considered.Agra et al. (2013) studied the VRPTW with uncertain travel times and proposed two robust formulations of the problem.The first extends the formulation of inequalities of resources and the second generalizes the formulation of inequalities of way.Their results show that the solution times are similar for both formulations while being significantly faster than the solutions times of a layered formulation recently proposed for the problem.This work is devoted to the robust VRPTW including both uncertainties in travel times and service times.It is worth mentioning that incorporating service time as a source of uncertainty in addition to the travel times is considered in this work, for the first time to our knowledge.Our contribution to all previous works, lies in the introduction of an effective way of modeling uncertain data, the choice of a mathematical model and the methods to evaluate that robust solution that can solve the whole problem, and also the selection of an adaptive large neighborhood search heuristic to solve each problem related to each scenario (with the use of Monte-Carlo Simulation to generate scenarios).
The remainder of this paper is structured in the following way.First of all, we define the problem, and we introduce its mathematical formulation.Next, we present our robust approach that is meant to solve the problem including a presentation of the ALNS preprocessing, destruction and insertion heuristics.Moreover, we evaluate the new approach using both the exact algorithm and the best-known heuristics.
A detailed computational and comparative study is presented in Section 4 in order to provide perfectly robust conclusions.Finally, some concluding remarks are discussed.

Problem statement
This section is devoted to the statement of the vehicle routing problem with time windows under travel time and service time uncertainty.First, we introduce the deterministic model of the VRP problem which consists of an optimization of the total distance traveled by all vehicles under four constraints.Next, the service at any customer starts within a given time interval and it is not allowed to arrive late.Furthermore, if the vehicle arrives too early at a customer, it must wait until the time window opens.Taking into consideration these two constraints on time windows we transform the VRP problem to its VRPTW variant.To complete our statment of the problem, we introduce the source of uncertainties namely travel times and service times which makes the problem harder to solve than its deterministic counterpart.We suggest a new formulation of the uncertainty that was inspired by the work of Bertsimas and Sim (2007) including only the travel time which belongs to a demand uncertainty polytope.
The complexity of this problem leads us to look for robust solutions and therefore to min-max the objective function, this is the last part of the state of the art of our problem.Now, in order to describe our problem, let us denote the set of nodes by N, using i and j to denote general nodes, the depot will be denoted by o.The set of arcs is denoted as A and contains pairs of nodes, (i, j).The set of vehicles is called V with elements k.Now we can assign to each edge (i, j) a cost  , and to each node i a time window [𝑎 ,𝑏 ].Then  are binary decision variables that take the value 1 if vehicle k uses the edge (i, j) and 0, otherwise.The deterministic model of the VRP can be stated as follows: Each customer must be visited once, which is ensured by the first constraint.The second constraint ensures that each tour passes through the depot.The constraint (3) is a flow conservation constraint.Finally, the last constraint guarantees that each tour ends at the depot.Since the service time  at any client i by vehicle k begins inside a given time interval [ , ], we require an additional constraint.

𝑎 ≤ 𝑃 ≤ 𝑏 ∀(𝑖 ∈ 𝑁) ∀(𝑘 ∈ 𝑉).
(5) The time windows considered here is hard, i.e. they cannot be violated, if the vehicle arrives earlier than required at a client i, it must hold up until the time window [ , ] opens and moreover it is not permitted to arrive late.
where M is a great value.To model the uncertainty in travel times and service times in the presence of time windows, a step-wise (layered) formulation is used.Based on the approach of Bertsimas and Sim (2007), we assume that the travel times and service times are uncertain and that they take their values respectively in the intervals [ ,  +  ] and [ ,  +  ], where  and  are the nominal values,  and  represent the maximum deviations.We also define the sets of uncertainties associated with these times by: and where Г and Λ are two degrees of uncertainties defined to control the number of travel times and service times uncertain.They vary respectively between 0 and || + ||, and 0 and ||.Thus, when Г=0 and Λ=0 the robust case coincides with the case deterministic and when Г=|| + || and Λ=|| is the worst case where all travel times and service times are supposedly uncertain and simultaneously reach their worst values.Robust optimization seeks to obtain good solutions for all the possible realizations of the uncertainties without needing to define the laws of probability and considering only the nominal values and the possible deviations of the uncertain data.We introduce uncertainties by replacing the function objective by: And the constraint (6) treating the time windows by this: where  and  are two indicator functions. takes the value 1 when i ∈  and  takes 1 when (i, j) ∈  3. New robust approach for the VRPTW Along these lines, we propose an adaptive large neighborhood search (ALNS) heuristic to integrate into our approach in order to deal with robust VRPTW.The proposed metaheuristic (ALNS) is an extension of the Large Neighborhood Search (LNS) heuristic, which was first introduced by Shaw (1998), ALNS is a metaheuristic proposed by Ropke and Pisinger (2006).It is a common technique used to enhance a locally optimal solution and can prevent getting stuck in premature convergence to local optima within tightly constrained search space.Given an initial solution obtained by a construction method, it is based on the idea of improving the initial solution by applying various destroy and repair operators to generate large neighborhoods through which the search space is explored (Palomo-Martínez Pamela et al., 2017).The ALNS has already been adapted to several transportation problems including vehicle routing (Ropke & Pisinger, 2006), arc routing (Angélica Salazar-Aguilar et al., 2012), inventory-routing (Coelho et al., 2012), and the reliable multiple allocation hub location problem (Chaharsooghi et al., 2017).The ALNS seems well-suited for the VRPTW, its power is manifested in the fact that each new solution is obtained by first removing a number of vertices, then re-inserting these vertices into the solution.ALNS was chosen because it outperforms other mono-objective algorithms applied to the same problem while keeping the simplicity and high performance that characterize local search algorithms.

Adaptive Large Neighborhood Search
We will now describe the ALNS that we have used in the present paper.We believe that ALNS can be applied to a large class of difficult optimization problems.In order to design an ALNS algorithm for a given optimization problem we need to: -Choose a number of fast construction operators which are able to construct a full solution.
-Select a number of destruction operators.It might be sufficiently important to choose the destruction operators that are expected to work well with the construction operator, but it is unnecessary.Here is the detailed algorithm:

Initial solution generation
In order to deal with the initial solution, we apply a greedy algorithm, which will be used in the reconstruction phase of the ALNS.This operator aims to insert the non-inserted nodes by testing the different possible configurations and then giving a feasible solution.It is not necessary that the completion time of the initial solution be minimal, as this solution will be further enhanced using the ALNS method.

Solution destruction
During the destruction phase, we put forward three different removal methods to maintain the diversity during the searching process and to define the neighborhood to explore at each iteration.Each removal method aims to remove a predefined number of nodes.The first operator is known as proximity operator.Its objective is to select close clients according to a spatio-temporal measure (Shaw, 1998), and then remove clients engendering the higher value of this measure.Using the same technique, the route portion operator comes to give more exibility than the proximity operator to change the routes.The principle consists in choosing a pivot client owned to a road and remove it as well as its adjacent clients.Then, we calculate the spatio-temporal measure, with the objective to select a second client belonging to another route but close to the initial pivot.The second pivot is removed from the solution as well as its adjacent clients and so on until all clients will be removed.The third operator is referred to as longest detour operator.The interest of this operator is to remove the customers that lead to the largest cost increases for servicing them.For more details, we refer the reader to (PrescottGagnon et al., 2009).The algorithms of the used destroy operators can be found in Annexe.

Solution reconstruction
Solomon's insertion heuristic (1987) presented a technique for choosing the new customer to be inserted into a route using two criteria  (, , ) and  (, , ) to select customer  for insertion between adjacent clients  and  in the current partial route.The primary criterion  calculate the best feasible insertion place in the current route for each unrouted client  as  (),  , () = min Customer  * is then inserted into the route between ( * ) and ( * ).The measurement of an insertion place  depends on factors: the increase in total distance of the current route after the insertion, and the delay of service start time of the customer following the new inserted customer.To be more precise,  (, , ) is calculated as:  (, , ) =   +  −  +   −  ,  +  = 1,  ≥ 0,  ≥ 0,  ≥ 0 where  +  is the new distance between two nodes  and  after the insertion,  is the previous service start time,  is the old distance between  and  and  is the new service start time of customer  after the insertion of customer .The criterion  (, , ) is calculated as following  (, , ) =  −  (, , ),  ≥ 0 where the parameter  is used to define how much the best insertion place for an unrouted customer depends on its distance from the depot and extra time required to visit the customer by the current vehicle.

Roulette wheel
For each iteration of the destruction phase, a roulette-wheel procedure is applied to select a method for generating the neighborhood (nodes to be removed).During the search process, the ALNS maintains a score φ which measures the best performance of an heuristic  in the past iterations.The roulette wheel selection consists in selecting an heuristic  with a probability ∑ .During M iterations, the score φ is reset and the probabilities of choosing an heuristic are recalculated (Pisinger and Ropke (2007)).

ALNS applied to the robust VRPTW
In this section, we apply the adaptive large neighborhood search (ALNS) to solve the VRPTW, assuming that the displacement and the service times are both objects of uncertainty.The robustness of this approach has been tested on several scenarios generated by the Monte Carlo tool of simulation.We will now describe how we have adapted the general ALNS to the robust VRPTW.Our goal is to provide, for each pair of degrees of robustness (, Г) considered, a robust solution that best protects from the violation of time windows, or a solution that minimizes the total delay compared to the dates at the latest where Г and  are two degrees of uncertainty defined to control the number of uncertain displacement times and the number of uncertain service times.They vary respectively between 0 and || + ||, and 0 and ||.Thus, when Г = 0 and  = 0, the robust case coincides with the deterministic case, and when Г = || + || and  = || is the worst case where all the displacement times and all service times are assumed uncertain and simultaneously reach their worst values.
For each pair of degrees of robustness (, Г) given, we first generate a set of possible realizations.Each realization  ,Г is defined by the assignment of Г displacement time ( service time) to their maximal values, and || − Г (resp.||− ) that remain at their nominal values.Then, on each achievement, we seek a robustly feasible solution or a solution with a minimum total delay.The reached solutions for the realizations considered are evaluated on all possible realizations and at the end of each iteration, we keep the solution that offers the worst minimum assessment.We note that the diversification of solutions is ensured by having solutions calculated independently from different realizations.Our algorithm is presented in detail in the subsections: (3.

Generation of realizations
An iteration of the robust approach algorithm starts with the generation of a set of realizations (by Monte-Carlo Simulation), each realization  ,Г represents a possible scenario in which the displacement times associated with a subset of arcs  ⊂  of cardinality Г take their maximum values  =  +  , and service times of a subset of vehicles  ⊂  of cardinality  take their maximum values  , =  +  .While the other arcs and the other nodes take respectively their nominal values  and  .

Research of solution
For each realization  ,Г , we apply the Adaptive large neighborhood search in order to obtain a feasible solution noted  satisfies each scenario that we have already generalized by Monte-Carlo.

Check of robustness
Even if the solution  is achievable on the realization  ,Г , it may violate window time constraints if it considers other realizations.Thus, to verify the robustness of this solution we apply Algorithm 2, where we seek to verify the robustness of the solution without the need to test it on all possible realizations.Indeed, a solution  is robustly feasible if its tours respect the windows of time at any node visited, where at most, Г displacement times and  service times are uncertain.
Formally, a tour  is robustly feasible, if and only if, on each path  ∀ℎ ∈ 2,3, … , , the maximum arrival date  does not violate time windows to the last node.The displacement times are determined in distinguishing between two cases: the case where the degree of uncertainty Г is greater than or equal to the number of arcs of the path  .In this case, we consider only the worst realization that can arise, where all the displacement times associated with the arcs of path  take their maximum values.On the other hand, in the case where Г is less than the number of arcs of the path  , we assign to the Г arcs of the set  Г, having the greatest deviations the maximum values, and to arcs that do not belong to this set the nominal values.The same procedure is applied to calculate the service times.

Evaluation on the worst case
In this step, we evaluate the solution  on the worst of possible cases, which corresponds to the realization where the displacement times associated to the Г arcs with the worst deviations and belonging to this solution, reach simultaneously their maximum values.First, we put in descending order all the arcs of the solution according to their maximum deviations.Then, we assign to the first Г arcs their maximum displacement time and the remaining arcs nominal displacement time.Finally, we carry out a summation of the times obtained to determine the worst-case evaluation of cases.This step is summarized in Algorithm 4.
Algorithm 4 Evaluation on the worst case  Г ( ) ← 0 Put in descending order all the arcs of ( ) according to their maximum deviations.

Computational experiment
Since VRPTW and RVRPTW are both NP-Hard, so to provide perfect conclusions and comparative results, we considered several kinds of instances.The robust approach examined was tested on a set of small instances based on the reference of Solomon benchmark (1987) (Solomon, 1987), and large instances of Gehring & Homberger's benchmark.Since the uncertainty of RVRPTW is simulated by discrete scenarios using Monte-Carlo Simulation, the uncertain travel time of each arc and the uncertain service time at each node are generated randomly between 0 and 10, with (Г, ) is the degree of robustness which represents the number of service times and the number of travel times assumed to be uncertain.The used instances are noted as follow _Г__, where  = {1, 2, 1, 2, 1, 2} corresponds to the class name of the benchmark of Solomon and Gehring & Homberger.Г and  represent the number of travel times and service times supposed uncertain. = {100,200,400,600,800,1000} is an index that represents the size of the instance.
Table 1 shows the results obtained for small instances (100 customers) by using Cplex for the deterministic VRPTW and the results obtained by our robust approach based on ALNS that deals with the VRPTW considering that travel times and service times are both uncertain.The column "Instance" displays the notation of the instance.The column "Initial solution" presents the initial solution with which the robust approach starts.The column "best" states the best values found by the robust approach with 10 runs while the column "mean" shows the average values found by the robust approach over 10 trials.The column "Optimal" displays the optimal solution for the deterministic VRPTW calculated by Cplex.Table 2 shows the results obtained for large instances by comparing the best-known results for the deterministic VRPTW to the results found by our robust approach based on ALNS that deals with the VRPTW considering that travel times and service times are both uncertain.

Table 2
Performance of our robust approach versus best known results

Table 2
Performance of our robust approach versus best known results (Continued) In order to visualize the impact of increasing degrees of uncertainty on objective function, we set the value of Г to 25 and we adjusted the value of Λ for multiple instances of size 100.Here is the curve obtained: The graph shows clearly that the objective function increases according to degree of uncertainty.To the best of our knowledge, this contribution is the first work to be devoted to the study of VRPTW considering both the uncertainties on travel times and service times.Due to the lack of work in this direction, we compared our results with the deterministic VRPTW literature even if the two problems are different in the sense that a good solution found for the deterministic case could become worse in the presence of uncertainties or even unreachable.
In fact, the robustness of this approach has been tested on several data sets and the results showed that the robust solutions offer great protection against delays with a slight increase in travel times and service times compared to what would have been found if a deterministic solution had been applied.The developed algorithm offers decision-making tool that allows them to choose, according to their specifications, the level of protection as well as the solution to apply.It is then clear that our approach is very powerful in terms of the robustness since it included several algorithms (Robustness verification, worst-case evaluation ...) which leads to near best solutions for all the possible realizations of the uncertainties without any further considerations but only nominal values and deviations possible uncertain data are sufficient.

Conclusion
Our main goal in this paper was to consider the robust vehicle routing problem with time windows under both travel times and service times uncertainties.For this purpose, a new robust approach has been suggested to minimize the total distance of the travel time in the presence of the maximum deviations of possible uncertain data.In this contrast, we generated all possible scenarios by using Monte Carlo simulation and we opt for the adaptive large neighborhood search ALNS algorithm to solve each subproblem related to each scenario.In this context, several destroy/repair method has been combined to explore multiple neighborhoods within the same search and defined implicitly the large neighborhood.
In order to study the feasibility of the resulting solution, efficient mechanisms have been conceived; the first concerns the verification of the robustness, while the second takes into consideration the evaluation of the solution on the worst case.
This new approach lying in the introduction of an effective way of modeling and handling several uncertainty data levels defined by pairs of uncertainty (, Г), which represent respectively the number of service times and the number of travel times assumed uncertain, has been tested on several sets of problems and showed improved robustness results for benchmark instances.Furthermore, this method offers great protection against delays with a slight increase in travel times and service times compared to what would have been found if a deterministic solution had been applied.
In this work, the computational experiments were performed to examine our proposed new approach compared to the deterministic VRPTW literature on a set of small instances based on the Solomon VRPTW benchmark and large instances of Gehring & Homberger benchmark.The results have shown the robustness of our solutions against delays and offer decision-making tool that allows choosing the level of protection, as well as the deterministic solution, is applied.
Future work will focus on the extension of the robust routing vehicle problem with time windows, in which both unexpected delays in travel time and service time may occur, to the application of parallel adaptive large neighborhood search in order to develop fast optimization procedures able to react in real time to changes in the problem information.
2.1), (3.2.2), (3.2.3) and (3.2.4).Here are a few notations used in our Algorithm:  ,Г : A realization possible  : The best robust solution  : The solution found to the  realization Cost (.):A function used to calculate the total time of displacement of a solution Weval Г (. ): A function used to calculate the worst evaluation of a solution  = ( = ,  , … ,  = ): The tour of the vehicle   = ( ,  , … ,  ): A path of the tour   Г, : The whole of the arcs which have the Г more large deviations of travel time  , : The whole of the nodes which have the  more large deviations of service time ( ) = { ,  , … ,  }: All of the nodes which constitute  ( ) = { = ( ,  ),  = ( ,  ), … ,  = ( ,  )}: The whole of the arcs which constitute the path  ( ): The maximum date of arrival of the vehicle  at node Here is the detailed algorithm of our approach:

Fig. 1 .
Fig. 1.Values of the objective function when adjusting the value of Algorithm 1 Adaptive Large Neighborhood Search Construct a feasible solution ; set  * =  Repeat Choose a destroy neighborhood  and a repair neighborhood  using roulette wheel selection based on previously obtained scores  Generate a new solution  from  using the heuristics corresponding to the chosen destroy and repair neighborhoods If  can be accepted then

Table 1
Performance of our robust approach versus deterministic VRPTW (CPLEX)